####Caruso T and Rillig MC 2022
###Simulation codes for paper submitted to:
###'Plant and Soil'; special issue on 'Plant-soil feedback: the next generation'. 



##########################
###Bever J 2002. Soil community feedback and the coexistence of competitors: conceptual frameworks and empirical tests
### New Phytologist (2003) 157: 465-473

############The following function correponds to Bever's deterministic model
############equations eq.1 (for plants) and eq. 3 (for the soil community)

###Plant A is pA
###Plant B is pB
###the soil community axis is sA
###sA means that when the soil axis = 1, the soil community is the one that plant A alone would be associated to



#############################################
########Solution of deterministic model using the function "ode" in "deSolve"

#############Model definition
BeverPS<-function(t, state, parameters) {
  with(as.list(c(state, parameters)),{
    # rate of change
    
    dpA <- (ra*pA)*((1+(alpha_A*sA))+(beta_A*(1-sA))-((pA+cB*pB)/KA))
    dpB <- (rb*pB)*((1+(alpha_B*sA))+(beta_B*(1-sA))-((pB+cA*pA)/KB))
    dsA <- sA*(1-sA)*( (pA/(pA+pB)) - v*(pB/(pA+pB)) ) 
    
    # return the rate of change
    list(c(dpA,dpB,dsA))
  }
  ) # end with(as.list ...
}


# time variable for integration
times <- seq(0, 1000, by = 0.01)

###the parameters, set  identical to Bever 2002, with the slight variation of having equal KA and KB
parameters<-c(ra=0.7,rb=0.5,
              
              alpha_A=-0.03,
              alpha_B=0.1,
              
              beta_A=0.1,
              beta_B=-0.2,
              cB=0.98,
              cA=0.885,
              KA=100,
              KB=100,
              v=0.8)

# initial states
state <- c(pA=50,pB=50,sA=0.5)

####solution
library(deSolve)
out<-ode(y = state, times = times, func = BeverPS, parms = parameters)

###extract solution
simBeverPS<-as.data.frame(out)


###Plot results

##Replot results into one panel
par(mfrow=c(1,1),mar=c(4.5,5.75,1.5,1),tck=-.03,mgp=c(3,.5,0),cex.axis=0.7)

plot(simBeverPS$time+1,simBeverPS$pA,type="n",cex.lab=1.5,cex.axis=1.5,xlab="time",ylab="Population size",main="Temporal Trajectory",ylim=c(30,80))

#######Plant A in default black
lines(simBeverPS$time+1,simBeverPS$pA)

#######Plant B in blue
lines(simBeverPS$time+1,simBeverPS$pB,col="blue")

####Soil axis in pink, in "percentage" on the y-axis
lines(simBeverPS$time+1,100*simBeverPS$sA,col="pink")








############################
######################
#############Stochastic version

#########################
##################
####Definition of key function 
####with iterative scheme based on the Euler-Maruyama method 

Bever_Stoch_Paths<-function(para, x_ini, Time, steps){
  n<-steps
  pA<-numeric(n) ###plant A
  pB<-numeric(n) ###plant B
  sA<-numeric(n) ###soil
  z1<-numeric(n) ###noise 1
  z2<-numeric(n) ###noise 2
  z3<-numeric(n) ###noise 3
  pA[1] <- x_ini[1]
  pB[1] <- x_ini[2]
  sA[1] <- x_ini[3]
  
  for (i in 1:n) {
    ####using normal distribution for white noise
    z1[i]<-rnorm(1)
    z2[i]<-rnorm(1)
    z3[i]<-rnorm(1)
###the equations with drift equals to Bever's model and diffusion
###based on Cholesky factorization of the variance-covariance matrix 
    pA[i+1]<- pA[i]+((para["ra"]*pA[i])*((1+(para["alpha_A"]*sA[i]))+(para["beta_A"]*(1-sA[i]))-((pA[i]+para["cB"]*pB[i])/para["KA"])))*(Time/steps)+para["sigmaA"]*sqrt(Time/steps)*z1[i]
    pB[i+1]<- pB[i]+((para["rb"]*pB[i])*((1+(para["alpha_B"]*sA[i]))+(para["beta_B"]*(1-sA[i]))-((pB[i]+para["cA"]*pA[i])/para["KB"])))*(Time/steps)+para["sigmaB"]*sqrt(Time/steps)*(para["rho12"]*z1[i] + para["sigmaB"]*z2[i]*sqrt(Time/steps)*(sqrt(1-(para["rho12"]*para["rho12"]))))
    sA[i+1]<- sA[i]+((sA[i]*(1-sA[i]))*( (pA[i]/(pA[i]+pB[i])) - para["v"]*(pB[i]/(pA[i]+pB[i])) ))*(Time/steps)+z1[i]*para["sigmaS"]*sqrt(Time/steps)*para["rho13"]+z2[i]*para["sigmaS"]*sqrt(Time/steps)*(((para["rho23"]-para["rho13"]))/sqrt(1-(para["rho12"]*para["rho12"])))+para["sigmaS"]*sqrt(Time/steps)*z3[i]*sqrt(1-para["rho13"]^2-(((para["rho23"]-para["rho13"])^2)/(1-para["rho12"]^2)))
    
    
    }
  
  output<-data.frame(pA,pB,sA)
  output
  
}


############################################
###################
#Simulation

####Scenario 1: deterministic, all variances and correlation null
###set parameters
para<-c(ra=0.7,rb=0.5,
        
        alpha_A=-0.03,
        alpha_B=0.1,
        
        beta_A=0.1,
        beta_B=-0.2,
        
        cB=0.98,
        cA=0.885,
        KA=100,
        KB=100,
        v=0.8,
        sigmaA=0,
        sigmaB=0,
        sigmaS=0,
        rho12=0,
        rho23=0,
        rho13=0
)
###note, sigmaS should be < 0.1 for stable systems
###note the correlation and variance terms set to zero, to collapse the diffusion term
###set initial conditions
x_ini<-c("pA"=50,"pB"=50,"sA"=0.5)


Sim_stoch<-Bever_Stoch_Paths(para, x_ini, Time=1000, steps=10000)
######
par(mfrow=c(1,1),mar=c(4.5,5.75,1.5,1),tck=-.03,mgp=c(3,.5,0),cex.axis=0.7)

Time<-1000
steps<-10000
plot(seq(0,Time,Time/steps),Sim_stoch$pA,type="n",cex.lab=1.5,cex.axis=1.5,xlab="time",ylab="population size",ylim=c(30,80))
lines(seq(0,Time,Time/steps),Sim_stoch$pA)

#plot(seq(0,Time,Time/steps),Sim_stoch$pB,type="n",cex.lab=1.5,cex.axis=1.5,xlab="time",ylab="pB",main="Temporal Trajectory")
lines(seq(0,Time,Time/steps),Sim_stoch$pB,col="blue")

#plot(seq(0,Time,Time/steps),Sim_stoch$sA,type="n",cex.lab=1.5,cex.axis=1.5,xlab="time",ylab="sA",ylim=c(-0.5,1))
lines(seq(0,Time,Time/steps),100*Sim_stoch$sA,col="brown")
abline(h=50)




legend("bottomright", title="PSF",
       legend=c("Soil","Plant B","Plant A"),
       lwd=2, cex=0.7,y.intersp=0.8,col=c("pink","blue","black","lightgreen","pink"), lty=c(1,1,1),bty = "n")

#############################
########same ouput as for the deterministic version solved with "ode"


####Scenario 2: stochastic, all equal variances but no correlation
###set parameters
para<-c(ra=0.7,rb=0.5,
        
        alpha_A=-0.03,
        alpha_B=0.1,
        
        beta_A=0.1,
        beta_B=-0.2,
        
        cB=0.98,
        cA=0.885,
        KA=100,
        KB=100,
        v=0.8,
        sigmaA=1,
        sigmaB=1,
        sigmaS=0.01,
        rho12=0,
        rho23=0,
        rho13=0
)
###note, sigmaS should be < 0.1 for stable systems
###set initial conditions
x_ini<-c("pA"=50,"pB"=50,"sA"=0.5)


Sim_stoch<-Bever_Stoch_Paths(para, x_ini, Time=1000, steps=10000)
######

#####Ten paths per species
sims<- replicate(10,Bever_Stoch_Paths(para, x_ini, Time=100, steps=10000), simplify=FALSE)
str(sims)

###############
#######collect multiple trajectories to plot them and calculate stats
sp_trajs<-function (simsobj,sp_n) {
  n<-length(simsobj)
  traj_list<-vector('list',n)
  for(k in 1:n) traj_list[[k]]<-as.matrix(simsobj[[k]])[,sp_n]
  traj_list
  
  }  

pA_trajs<-sp_trajs(simsobj=sims,sp_n=1)
str(pA_trajs)
pA_trajsm<-do.call(cbind, pA_trajs)

par(mfrow=c(3,1),mar=c(4.5,5.75,1.5,1),tck=-.03,mgp=c(3,.5,0),cex.axis=0.7)


Time<-1000
steps<-10000
pA_trajs<-sp_trajs(simsobj=sims,sp_n=1)
str(pA_trajs)
pA_trajsm<-do.call(cbind, pA_trajs)
matplot(seq(0,Time,Time/steps),pA_trajsm,type="l",xlab="Time",ylab="Pop Size",main="multiple paths plant A") 

Time<-1000
steps<-10000
pB_trajs<-sp_trajs(simsobj=sims,sp_n=2)
str(pB_trajs)
str(pA_trajs)

pB_trajsm<-do.call(cbind, pB_trajs)
matplot(seq(0,Time,Time/steps),pB_trajsm,type="l",ylab="Pop Size",xlab="Time",main="multiple paths plant B") 



Time<-1000
steps<-10000
sA_trajs<-sp_trajs(simsobj=sims,sp_n=3)
str(sA_trajs)

sA_trajs<-do.call(cbind, sA_trajs)
matplot(seq(0,Time,Time/steps),sA_trajs,type="l",xlab="Time",ylab="Soil Community Axis",main="multiple paths sA") 





####Scenario 3: stochastic, all equal variances and correlation
###set parameters
para<-c(ra=0.7,rb=0.5,
        
        alpha_A=-0.03,
        alpha_B=0.1,
        
        beta_A=0.1,
        beta_B=-0.2,
        
        cB=0.98,
        cA=0.885,
        KA=100,
        KB=100,
        v=0.8,
        sigmaA=1,
        sigmaB=1,
        sigmaS=0.01,
        rho12=0.5,
        rho23=0.5,
        rho13=0.5
)

###set initial conditions
x_ini<-c("pA"=50,"pB"=50,"sA"=0.5)


Sim_stoch<-Bever_Stoch_Paths(para, x_ini, Time=1000, steps=10000)
######

#####Ten paths per species
sims<- replicate(10,Bever_Stoch_Paths(para, x_ini, Time=100, steps=10000), simplify=FALSE)
str(sims)

###############
#######collect multiple trajectories to plot them and calculate stats
sp_trajs<-function (simsobj,sp_n) {
  n<-length(simsobj)
  traj_list<-vector('list',n)
  for(k in 1:n) traj_list[[k]]<-as.matrix(simsobj[[k]])[,sp_n]
  traj_list
  
}  

pA_trajs<-sp_trajs(simsobj=sims,sp_n=1)
str(pA_trajs)
pA_trajsm<-do.call(cbind, pA_trajs)

par(mfrow=c(3,1),mar=c(4.5,5.75,1.5,1),tck=-.03,mgp=c(3,.5,0),cex.axis=0.7)


Time<-1000
steps<-10000
pA_trajs<-sp_trajs(simsobj=sims,sp_n=1)
str(pA_trajs)
pA_trajsm<-do.call(cbind, pA_trajs)
matplot(seq(0,Time,Time/steps),pA_trajsm,type="l",xlab="Time",ylab="Pop Size",main="multiple paths plant A") 

Time<-1000
steps<-10000
pB_trajs<-sp_trajs(simsobj=sims,sp_n=2)
str(pB_trajs)
str(pA_trajs)

pB_trajsm<-do.call(cbind, pB_trajs)
matplot(seq(0,Time,Time/steps),pB_trajsm,type="l",ylab="Pop Size",xlab="Time",main="multiple paths plant B") 



Time<-1000
steps<-10000
sA_trajs<-sp_trajs(simsobj=sims,sp_n=3)
str(sA_trajs)

sA_trajs<-do.call(cbind, sA_trajs)
matplot(seq(0,Time,Time/steps),sA_trajs,type="l",xlab="Time",ylab="Soil Community Axis",main="multiple paths sA")