#R code for the paper "Dynamic economic thresholds for insecticide applications 
#against agricultural pests: importance of pest and natural enemy migration

#Authored by Tamar Keasar, Eric Wajnberg, George Heimpel, Ian Hardy, 
#Liora Shaltiel, Daphna Gottlieb, Saskya van Nouhuys

#This code was written by Eric Wajnberg and Tamar Keasar

#It models the number of insecticide applications needed to keep a pest population below 
#a fixed economic threshold or a dynamic threshold adjusted to the abundance of a natural enemy
#It also calculates the farmer's seasonal revenue.

#This version is from September 2022


library(plotly)

setwd("d:/Projects/IIAS 2022/Pesticides team")
pesticide=function (r,beta,m,alpha,P0,B0,Ip, Ib,Im,Cp,C_NE,K, Cost.spray, Yield.plant, Crop.value,
                    Yield.reduction, Value.kg, Plants.hectare)
{
  # procedure for the pesticide model for Tamar Keasar
  # Eric Wajnberg - March 2022  #
  # Here is the meaning of the parameters:
  #
  # r: Pest reproductive rate/week (default: 1)
  # beta: NE birth rate/week (default: 4.9)
  # m:    NE mortality/week (default: 4.2)
  # alpha: NE's predation rate (default: 2.2)
  # P0: Initial pest population size (default: 10)
  # B0: Initial NE population size (default: 2)
  # Ip: Proportion of pests killed by pesticide (default: 0.75)
  # Ib: Proportion of NEs killed by pesticide (default: 0.2)
  # Im: Proportion of NEs migrating after pesticide (default: 0)
  # Cp: Out-field pests colonization parameter(default:0.001)
  # C_NE: Out-field NE colonization parameter (default: 0.001)
  # EIL: Economic injury level
  # K: Field carrying capacity for pests (default: 50)
  # Cost.spray: Pesticide application cost per hectare (default: USD 120)
  # Crop.value: Market price of 1 kg crop (default: USD 0.12 for tomato)
  # Yield.reduction: damage to yield per pest individual (in percent of yield of a single plant, default: 0.0121)
  # Plants.hectare: Planting density (default: 30000 plants/hectare for tomato)
  #
  # scenario 1: When pest number > Economic Threshold (EIL/(1+r))
  # scenario 2: When pest number > Economic Injury Level (>EIL/((1+r)-alpha*B))
  #
  
  # week numbers
  weeks=0:24
  
  #Calculate EIL
  EIL=Cost.spray/(Yield.plant*Crop.value*Plants.hectare*Yield.reduction*Ip)
  
  #Calculate ET
  Fixed.threshold=EIL/(1+r)
  
  #Make array for DTs
  Dynamic.threshold=rep(NA,length(weeks))
  
  # computing both scenarios
  Pt.total=matrix(NA,length(weeks),2)
  P.born=matrix(NA,length(weeks),2)
  P.killed.by.NE=matrix(NA,length(weeks),2)
  P.killed.by.pesticide=matrix(NA,length(weeks),2)
  P.moved.in=matrix(NA,length(weeks),2)
  Bt.total=matrix(NA,length(weeks),2)
  B.born=matrix(NA,length(weeks),2)
  B.died=matrix(NA,length(weeks),2)
  B.killed.by.pesticide=matrix(NA,length(weeks),2)
  B.moved.to.out.field=matrix(NA,length(weeks),2)
  B.moved.to.in.field=matrix(NA,length(weeks),2)
  Seasonal.pest=matrix(NA,2)
  Seasonal.revenue=matrix(NA,2)
  Seasonal.spray=matrix(NA,2)
  
  Spray=matrix(NA,length(weeks),2)
  
  #Set initial densities of pest and NE
  Pt.total[1,]=P0
  Bt.total[1,]=B0
  
  #No spraying on week 1
  Spray[,1]=0
  Spray[1,]=0
 
  for (i in 2:length(weeks))
  {
    Spray[i,1]=ifelse(Pt.total[i-1,1]>Fixed.threshold,1,0)
    Dynamic.threshold[i]=(EIL+alpha*Bt.total[i-1,2])/(1+r)
    Spray[i,2]=ifelse(Pt.total[i-1,2]>Dynamic.threshold[i],1,0)
    
    for (j in 1:2)
    {
      B.moved.to.in.field[i,j]=Pt.total[i-1,j]*C_NE
      B.moved.to.out.field[i,j]=ifelse(Spray[i,j]==0,0,Bt.total[i-1,j]*Im)
      B.killed.by.pesticide[i,j]=ifelse(Spray[i,j]==0,0,Bt.total[i-1,j]*Ib)
      B.died[i,j]=Bt.total[i-1,j]*m
      B.born[i,j]=Bt.total[i-1,j]*Pt.total[i-1,j]*beta
      Bt.total[i,j]=B.born[i,j]-B.died[i,j]-B.killed.by.pesticide[i,j]-
      B.moved.to.out.field[i,j]+B.moved.to.in.field[i,j]
      if(Bt.total[i,j]<0) Bt.total[i,j]=0
      P.moved.in[i,j]=(K-Pt.total[i-1,j])*Cp
      if(P.moved.in[i,j]<0) P.moved.in[i,j]=0
      P.killed.by.pesticide[i,j]=ifelse(Spray[i,j]==0,0,Pt.total[i-1,j]*Ip)
      P.killed.by.NE[i,j]=Bt.total[i-1,j]*alpha
      P.born[i,j]=Pt.total[i-1,j]*(1+r)
      Pt.total[i,j]=P.born[i,j]-P.killed.by.NE[i,j]-P.killed.by.pesticide[i,j]+P.moved.in[i,j]
      if(Pt.total[i,j]<0) Pt.total[i,j]=0
            }
  }
  for (j in 1:2)
    {
    Seasonal.pest[j]=sum(Pt.total[,j])
    Seasonal.spray[j]=sum(Spray[,j])
    Seasonal.revenue[j]=Yield.plant*(1-Seasonal.pest[j]*Yield.reduction)*Crop.value*Plants.hectare
    }
    
  Outputs <<-list(Seasonal.spray, Seasonal.revenue) 
}


Make.plot=function(For.plot, xtitle, ytitle, ztitle, Surf.color, Fig.name, zmin, zmax, step, meshcolor)
{
  axx= list(title=xtitle)
  axy=list(title=ytitle)
  axz=list(title=ztitle, nticks = step, range = c(zmin, zmax))
  colnames(For.plot)=c("xvalue", "yvalue", "zvalue")
    plot_ly(For.plot, x = ~xvalue, y = ~yvalue, z = ~zvalue, 
          type = 'mesh3d',
          facecolor = rep(meshcolor, 1000),
          opacity=0.4)%>%
   layout(scene=list(xaxis=axx, yaxis=axy, zaxis=axz))%>%
   layout(font=list(size = 15))
  
}

#Fig. 1, 2
#Manipulate pest and NE colonization rates

Figs12=data.frame(matrix(ncol = 4, nrow = 0))
for (x in c(0, 0.003, 0.006, 0.009, 0.012, 0.015))
{
  for (y in c(0, 0.003, 0.006, 0.009, 0.012, 0.015)) 
#These colonization values were selected because they result in positive revenues
    {
    pesticide(r=1,beta=4.9, m=4.2,alpha=2.2,P0=10,B0=2,Ip=0.75, Ib=0.2, Im=0, Cp=x, C_NE=y, 
    K=50, Cost.spray=120, Yield.plant=4, Crop.value=0.12, Yield.reduction=0.0121, Plants.hectare=30000)
    z=c(x,y,Outputs[[1]], Outputs[[2]])
    Figs12=rbind(Figs12, z)
    }
}
names(Figs12)[1]='Pest.colonization'
names(Figs12)[2]='NE.colonization'
names(Figs12)[3]='Sprays.ET'
names(Figs12)[4]='Sprays.DT'
names(Figs12)[5]='Revenue.ET'
names(Figs12)[6]='Revenue.DT'
Figs12[7]=Figs12[3]-Figs12[4]
names(Figs12)[7]='Sprays.diff'
Figs12[8]=Figs12[1]*1000
names(Figs12)[8]='Pest.colonization.1000'
Figs12[9]=Figs12[2]*1000
names(Figs12)[9]='NE.colonization.1000'


#Fig. 1

Fig1a=data.frame(Figs12[,c(8,9,3)])
Make.plot(For.plot=Fig1a, xtitle="Pest immigration rate", ytitle="Natural enemy immigration rate", 
     ztitle="Number of insecticide applications", zmin=0, zmax=24, step=6, meshcolor="black")
Fig1b=data.frame(Figs12[,c(8,9,4)])
Make.plot(For.plot=Fig1b, xtitle="Pest immigration rate", ytitle="Natural enemy immigration rate", 
    ztitle="Number of insecticide applications", zmin=0, zmax=24, step=6, meshcolor="black")
Fig1c=data.frame(Figs12[,c(8,9,7)])
Make.plot(For.plot=Fig1c, xtitle="Pest immigration rate", ytitle="Natural enemy immigration rate", 
    ztitle="Difference in number of insecticide applications", zmin=0, zmax=5, step=5, meshcolor="black")

#Fig. 2

Fig1black=data.frame(Figs12[,c(8,9,5)])
Make.plot(For.plot=Fig1black, xtitle="Pest immigration rate", ytitle="Natural enemy immigration rate", 
          ztitle="Seasonal revenue (US$)", zmin=0, zmax=10000, step=10, meshcolor="grey")


#Fig. 3, 4
#Manipulate insecticide toxicity to pest and NE
Figs34=data.frame(matrix(ncol = 4, nrow = 0))
for (x in seq(0.4, 1, 0.06))
{
  for (y in seq(0, 1, 0.1))
  {
#Revenue is positive at toxicity levels >0.4  
    pesticide(r=1,beta=4.9, m=4.2,alpha=2.2,P0=10,B0=2,Ip=x, Ib=y, Im=0, Cp=0.001, C_NE=0.001, 
   K=50, Cost.spray=120, Yield.plant=4, Crop.value=0.12, Yield.reduction=0.0121, Plants.hectare=30000)
    #Selected Cp and C_NE values that would produce positive revenues
    z=c(x,y,Outputs[[1]], Outputs[[2]])
    Figs34=rbind(Figs34, z)
  }
}
names(Figs34)[1]='Toxicity.pest'
names(Figs34)[2]='Toxicity.NE'
names(Figs34)[3]='Sprays.ET'
names(Figs34)[4]='Sprays.DT'
names(Figs34)[5]='Revenue.ET'
names(Figs34)[6]='Revenue.DT'
Figs34[7]=Figs34[3]-Figs34[4]
names(Figs34)[7]='Sprays.diff'


#Fig. 3

Fig3a=data.frame(Figs34[,c(1,2,3)])
Make.plot(For.plot=Fig3a, xtitle="Toxicity to pest", ytitle="Toxicity to natural enemy", 
    ztitle="Number of insecticide applications", zmin=0, zmax=24, step=6, meshcolor="black")
Fig3b=data.frame(Figs34[,c(1,2,4)])
Make.plot(For.plot=Fig3b, xtitle="Toxicity to pest", ytitle="Toxicity to natural enemy", 
    ztitle="Number of insecticide applications", zmin=0, zmax=24, step=6, meshcolor="black")
Fig3c=data.frame(Figs34[,c(1,2,7)])
Make.plot(For.plot=Fig3c, xtitle="Toxicity to pest", ytitle="Toxicity to natural enemy", 
  ztitle="Difference in number of insecticide applications", zmin=0, zmax=5, step=5, meshcolor="black")


#Fig. 4

Fig1black=data.frame(Figs34[,c(1,2,5)])
Make.plot(For.plot=Fig1black, xtitle="Toxicity to pest", ytitle="Toxicity to natural enemy", 
          ztitle="Seasonal revenue (US$)", zmin=0, zmax=10000, step=10, meshcolor="grey")


#Fig. 5
#Manipulate pest reproductive rate and NE feeding rate
Fig5=data.frame(matrix(ncol = 4, nrow = 0))
for (x in seq(1, 10, 1))
{
  for (y in seq(1, 10, 1))
  {
    pesticide(r=x,beta=4.9, m=4.2,alpha=y,P0=10,B0=2,Ip=0.75, Ib=0.2, Im=0, Cp=0.001, C_NE=0.001, 
    K=50, Cost.spray=120, Yield.plant=4, Crop.value=0.12, Yield.reduction=0.0121, Plants.hectare=30000)
    
    z=c(x,y,Outputs[[1]], Outputs[[2]])
    Fig5=rbind(Fig5, z)
  }
}
names(Fig5)[1]='Pest reproduction'
names(Fig5)[2]='NE feeding rate'
names(Fig5)[3]='Sprays.ET'
names(Fig5)[4]='Sprays.DT'
names(Fig5)[5]='Revenue.ET'
names(Fig5)[6]='Revenue.DT'
Fig5[7]=Fig5[3]-Fig5[4]
names(Fig5)[7]='Sprays.diff'

Fig5a=data.frame(Fig5[,c(1,2,3)])
Make.plot(For.plot=Fig5a, xtitle="Pest reproduction (Offspring/individual/week)", ytitle="NE feeding rate (individuals/week)", 
    ztitle="Number of insecticide applications", zmin=0, zmax=24, step=6, meshcolor="black")
Fig5b=data.frame(Fig5[,c(1,2,4)])
Make.plot(For.plot=Fig5b, xtitle="Pest reproduction (Offspring/individual/week)", ytitle="NE feeding rate (individuals/week)", 
    ztitle="Number of insecticide applications", zmin=0, zmax=24, step=6, meshcolor="black")
Fig5c=data.frame(Fig5[,c(1,2,7)])
Make.plot(For.plot=Fig5c, xtitle="Pest reproduction (Offspring/individual/week)", ytitle="NE feeding rate (individuals/week)", 
  ztitle="Difference in number of insecticide applications", zmin=0, zmax=5, step=5, meshcolor="black")


#Fig. 6
#Manipulate pest damage to the plants and the crop value
Fig6=data.frame(matrix(ncol = 4, nrow = 0))
for (x in seq(0.001, 0.05, 0.005))
{
  for (y in seq(0.01, 1, 0.1))
  {
    pesticide(r=1,beta=4.9, m=4.2,alpha=2.2,P0=10,B0=2,Ip=0.75, Ib=0.2, Im=0, Cp=0.001, C_NE=0.001, 
              K=50, Cost.spray=120, Yield.plant=4, Crop.value=y, Yield.reduction=x, Plants.hectare=30000)
    
    z=c(x,y,Outputs[[1]], Outputs[[2]])
    Fig6=rbind(Fig6, z)
  }
}
names(Fig6)[1]='Pest damage to yield'
names(Fig6)[2]='Crop value'
names(Fig6)[3]='Sprays.ET'
names(Fig6)[4]='Sprays.DT'
names(Fig6)[5]='Revenue.ET'
names(Fig6)[6]='Revenue.DT'
Fig6[7]=Fig6[3]-Fig6[4]
names(Fig6)[7]='Sprays.diff'

Fig6a=data.frame(Fig6[,c(1,2,3)])
Make.plot(For.plot=Fig6a, xtitle="Pest damage to yield", ytitle="Crop value", 
    ztitle="Number of insecticide applications", zmin=0, zmax=24, step=6, meshcolor="black")
Fig6b=data.frame(Fig6[,c(1,2,4)])
Make.plot(For.plot=Fig6b, xtitle="Pest damage to yield", ytitle="Crop value", 
    ztitle="Number of insecticide applications", zmin=0, zmax=24, step=6, meshcolor="black") 
Fig6c=data.frame(Fig6[,c(1,2,7)])
Make.plot(For.plot=Fig6c, xtitle="Pest damage to yield", ytitle="Crop value", 
  ztitle="Difference in number of insecticide applications", zmin=0, zmax=5, step=5, meshcolor="black")