################################################################################################ ### This file generates the simulated real system and corresponding observation data ### ### using a stochastic version of the Lotka Volterra Model with different initial conditions ### ### and time varying rate constants, ### ################################################################################################ ##################################################### ### Gillespie Function (with regular time output) ### ##################################################### gillespied <- function (N, T=10, dt=0.5, ...) # Core function to simulate any stochastic chemical reaction # network. Takes the list N which describes the possible # reactions, chemicals and initial conditions and simulates # a continuous time, discrete state markov process. See # Darren Wilkinson's book for more details. { tt=0 n=T%/%dt x=N$M S=t(N$Post-N$Pre) u=nrow(S) v=ncol(S) xmat=matrix(0,ncol=u,nrow=n) i=1 target=0 repeat { h=N$h(x,t=tt, ...) h0=sum(h) if (h0<1e-10) tt=1e99 else tt=tt+rexp(1,h0) while (tt>=target) { xmat[i,]=x i=i+1 target=target+dt if (i>n) return(xmat) # return(ts(xmat,start=0,deltat=dt)) } j=sample(v,1,prob=h) x=x+S[,j] } } ### Set up the Petri Net for the Lotka Volterra Model ### N=list() N$M=c(2000,800) N$Pre=matrix(c(1,0,1,1,0,1),ncol=2,byrow=TRUE) N$Post=matrix(c(2,0,0,2,0,0),ncol=2,byrow=TRUE) N$h=function(x,th=c(0.5,0.001,0.8),func_t,t) { return( c(th[1]*x[1], th[2]*x[1]*x[2], th[3]*x[2] ) * func_t(t)) # adjust rate constants as functions of time } ### setup the model time intervals etc. ### T_final <- 10 dt <- 0.01 tseq <- seq(0,T_final,dt) rep <- 4 ### set function to alter rate parameters theta over time func_t <- function(t){ 1 - ((t-5)^2/25 - 0.5) /10 } plot(tseq,func_t(tseq)) ### scale factor on pop size to control stochasticity ### sca_pop <- 4 ### run the model ### set.seed(1) arr1 <- array(0,c(length(tseq),2,rep)) for(i in 1:rep) { N$M=c(rnorm(1,sca_pop*2000,80),rnorm(1,sca_pop*800,32)) # random initial conditions if needed N$M=sca_pop*c(1910,710) # Initial conditions: overwrite with fixed ICs arr1[,,i] <- gillespied(N,T_final+dt+1e-8,dt,th=c(1,0.00044/sca_pop,1.8),func_t=func_t) cat(i,"\n") } arr1 <- arr1/sca_pop ### plot real system simulations ### plot(tseq,arr1[,1,1],ty="n",col=4,ylim=c(0,max(arr1)),xlab="time(t)",ylab="Pred/Prey Number") for(i in 1:rep) { plot(tseq,arr1[,1,1],ty="n",col=4,ylim=c(0,max(arr1))) lines(tseq,arr1[,1,i],col=4) lines(tseq,arr1[,2,i],col=2) cat(i,"\n") # readline() } # choose simulated reality number 2 # ch_sim_real <- 2 ### set up possible observed data set observed with error ### t_obs <- seq(0.5,10,0.5) al_obs <- match(t_obs,tseq) obs_sd <- 50 # observation error sd y_obs1 <- arr1[al_obs,,ch_sim_real] # real system output y set.seed(1) y_obs <- round(y_obs1 + rnorm(n=length(y_obs1),mean=0,sd=obs_sd)) # add on obs measurement error to get z = y + e colnames(y_obs) <- c("Prey","Pred") ### plot observed data for pred and prey with measurement error ### plot(t_obs,y_obs[,1],xlim=c(0,10),ylim=c(0,max(y_obs)*1.1),pch=16,col=4,cex=0.7) arrows(t_obs,y_obs[,1] + 3*obs_sd,t_obs,y_obs[,1] - 3*obs_sd,angle=90,code=3,col=4,length=0.03,lwd=2) points(t_obs,y_obs[,2],xlim=c(0,10),ylim=c(0,max(y_obs)),pch=16,col=2,cex=0.7) arrows(t_obs,y_obs[,2] + 3*obs_sd,t_obs,y_obs[,2] - 3*obs_sd,angle=90,code=3,col=2,length=0.03,lwd=2) ### save possible observed data ### save(y_obs,obs_sd,t_obs,file="Observed Data/Possible_Observed_Data_V1.Rdata") save(arr1,ch_sim_real,file="Observed Data/Possible_Observed_Data_full_array_V1.Rdata") # load("Observed Data/Possible_Observed_Data_V1.Rdata") # load("Observed Data/Possible_Observed_Data_full_array_V1.Rdata") ############################################################################################################################# ### Extract position and height of maximum peaks for second part of paper ### obs_Prey.T.1 <- tseq[which.max(arr1[1:500,1,ch_sim_real])] obs_Prey.T.2 <- tseq[501:1001][which.max(arr1[501:1001,1,ch_sim_real])] obs_Pred.T.1 <- tseq[which.max(arr1[1:500,2,ch_sim_real])] obs_Pred.T.2 <- tseq[501:1001][which.max(arr1[501:1001,2,ch_sim_real])] obs_Prey.Y.1 <- arr1[,1,ch_sim_real][which.max(arr1[1:500,1,ch_sim_real])] obs_Prey.Y.2 <- arr1[501:1001,1,ch_sim_real][which.max(arr1[501:1001,1,ch_sim_real])] obs_Pred.Y.1 <- arr1[,2,ch_sim_real][which.max(arr1[1:500,2,ch_sim_real])] obs_Pred.Y.2 <- arr1[501:1001,2,ch_sim_real][which.max(arr1[501:1001,2,ch_sim_real])] ### store 8 peak observed outputs in one object ### obs_maximums <- c(obs_Prey.T.1=obs_Prey.T.1, obs_Pred.T.1=obs_Pred.T.1, obs_Prey.Y.1=obs_Prey.Y.1, obs_Pred.Y.1=obs_Pred.Y.1, obs_Prey.T.2=obs_Prey.T.2, obs_Pred.T.2=obs_Pred.T.2, obs_Prey.Y.2=obs_Prey.Y.2, obs_Pred.Y.2=obs_Pred.Y.2) ### add peak points to plot to check they are indeed at the peaks ### points(obs_maximums[c(1:2,5:6)],obs_maximums[c(3:4,7:8)],pch=16,col="purple") lines(tseq,arr1[,1,ch_sim_real],col=4) lines(tseq,arr1[,2,ch_sim_real],col=2) ### define observed errors on peak outputs and add to plot ### obs_max_Y_sd <- 50 obs_max_T_sd <- 0.025 arrows(obs_maximums[c(1:2,5:6)]-3*obs_max_T_sd,obs_maximums[c(3:4,7:8)], obs_maximums[c(1:2,5:6)]+3*obs_max_T_sd,obs_maximums[c(3:4,7:8)],angle=90,code=3,len=0.02) arrows(obs_maximums[c(1:2,5:6)],obs_maximums[c(3:4,7:8)]-3*obs_max_Y_sd , obs_maximums[c(1:2,5:6)],obs_maximums[c(3:4,7:8)]+3*obs_max_Y_sd ,angle=90,code=3,len=0.05) ### save observed peak data ### save(obs_maximums,obs_max_Y_sd,obs_max_T_sd,file="Observed Data/Max_peaks_Obs_and_ObsErrors_V1.Rdata") # load("Observed Data/Max_peaks_Obs_and_ObsErrors_V1.Rdata") #############################################################################################################################