################################################################## # Created by Noa Pinter-Wollman, 2016 ################################################################## # Simulation code for the paper "The effect of keystone individuals on collective outcomes # can be mediated through interactions or behavioral persistence" # By Noa Pinter-Wollman, Carl N. Keiser, Roy Wollman, and Jonathan N. Pruitt ################################################################## ## Initialize: rm(list=ls()) # clean workspace graphics.off() # close all figures # load libraries library('igraph') library('timeordered') ################################################################## # set simulation parameters n=25 # number of spiders time = 150 # number of time steps iter=100 # number of simulation iterations KIboldness=c(10,300,600) # boldness of keystone individual pint=0.3 # mean of the uniform and skewed interaction probability pintcases=c('uniform','KIunif','skewed','KIskewed', 'changing') # interaction rules minM=2 # minimum number of interactions at each time step maxM=10 # maximum number of interactions at each time step prop_bold_to_give=0.2 # A - the proportion of the difference in boldness that will be acquired by the shyer individual. prop_bold_to_keep=0.95 # P - the proportion of boldness kept at each time step, i.e., behavioral persistence. # Parameters for collective outcomes: # Disease transmission: num_ts_4_net=20 # number of time steps to include in the network for disease transmission # Prey attack: time_att = 50 # number of time steps for attack iter_att= 100 # number of simulation iterations for attack # Function that determines the relationship between interaction probability and boldness for the ‘changing’ model: fbold2int = function(bold){return(bold/600)} # initialize book keeping: boldness=array(0, dim=c(n,time,iter, length(pintcases), length(KIboldness))) # boldness of all individuals, all simulations, time steps, interaction patterns and keystone boldness values. interactions=array(0, dim=c((time*maxM)-maxM,3,iter, length(pintcases), length(KIboldness))) # all interactions for all simulations, time steps, interaction patterns and keystone boldness values. num_int=array(0, dim=c(iter, length(pintcases), length(KIboldness))) # number of interactions in all simulations, interaction patterns, and keystone boldness values. ############################################################################# ## start simulation of boldness dynamics: for (g in 1:length(KIboldness)){ # start different keystone boldness loop KIbold=KIboldness[g] for(h in 1:length(pintcases)){ # loop of different interaction rules for (k in 1:iter){ # iteration loop # set the interaction rule if(pintcases[h]=='uniform'){pintall=rep(pint,n)} # all individuals are equally likely to interact if(pintcases[h]=='skewed'){pintall=rexp(n,1/pint)} # probability to interact follows an exponential (skewed) distribution if(pintcases[h]=='KIunif'){pintall=c(pint*2, rep(pint,n-1))} # all individuals are equally likely to interact but keystone is twice as likely to interact if(pintcases[h]=='KIskewed'){pintall=c(pint*2,rexp((n-1),1/pint))} # probability to interact follows an exponential (skewed) distribution, but the keystone is twice as likely to interact as the mean of the exponential distribution boldness[,1,k,h,g]=c(KIbold,rep(0, n-1)) # set starting boldness, first individual is always the keystone interactionstemp=c(NA,NA) # set up a temporary table for book keeping of interactions for (z in 2:time){ # time loop if(pintcases[h]=='changing'){pintall=(rep(pint,n)*(fbold2int(boldness[,z-1,k,h,g]))) # probability to interact is updated every time step in proportion to the boldness of the individual (based on function determined above) pintall=pintall+runif(n, max=0.0002) # add stochastic noise so that individuals with 0 boldness can interact too } M=round(runif(1, min=minM, max=maxM)) # determine the number of interactions in this time step spider1 = sample(n,M, replace = FALSE,prob = pintall) # select the first M spiders spider2 = sample(setdiff(1:n,spider1),M, replace = FALSE,prob = pintall[-spider1]) # select the second set of M spiders (excluding individuals that were selected as spider1 to avoid interactions with self) # record interactions and time step when they happened in a temp table: interactionstemp=rbind(interactionstemp,cbind(spider1,spider2,rep(z,M))) boldness[,z,k,h,g]=boldness[,z-1,k,h,g] # make sure boldness from previous step carries over for everyone # update spider boldness based on interaction: # find who has larger boldness and update the smaller boldness level - only the individual with the smaller boldness acquires more boldness, the individual with the greater boldness does not change its boldness. # the amount of boldness acquired by the shy individual is a proportion (prop_bold_to_give) of the difference in boldness between the two individuals. boldness[spider1,z,k,h,g]=boldness[spider1,z-1,k,h,g]+ (prop_bold_to_give*apply(cbind(rep(0,M),boldness[spider2,z-1,k,h,g]-boldness[spider1,z-1,k,h,g]),1,max)) boldness[spider2,z,k,h,g]=boldness[spider2,z-1,k,h,g]+ (prop_bold_to_give*apply(cbind(rep(0,M),boldness[spider1,z-1,k,h,g]-boldness[spider2,z-1,k,h,g]),1,max)) # forget some boldness: # everyone forgets some of their boldness (1-prop_bold_to_keep) except for the keystone (always the first individual) who never forgets. boldness[-1,z,k,h,g]=prop_bold_to_keep*boldness[-1,z,k,h,g] # -1 excludes the keystone } interactions[2:dim(interactionstemp)[1],,k,h,g]=interactionstemp[-1,] # keep interactions for later analysis num_int[k,h,g]=dim(interactionstemp)[1] # get the number of interactions that happened at this time step for later analysis } } } ##################################################################### ###################### Collective outcomes ########################## ##################################################################### ################################## ##### Disease transmission dynamics (using the Timeordered package): ################################## # initialize book keeping: T50_TO=array(0, dim=c(2,iter, length(pintcases), length(KIboldness))) # first dimension contains start and end networks Vmax_TO=array(0, dim=c(2,iter, length(pintcases), length(KIboldness))) # first dimension contains start and end networks T50_KI_TO=array(0, dim=c(2,iter, length(pintcases), length(KIboldness))) # first dimension contains start and end networks Vmax_KI_TO=array(0, dim=c(2,iter, length(pintcases), length(KIboldness))) # first dimension contains start and end networks for(i in 1:length(pintcases)){ # loop of different interaction rules for(j in 1: length(KIboldness)){ # loop of different keystone boldness values for (h in 1:iter){ # loop over simulation iterations ## work on following interactions: int_from_one_run=interactions[interactions[,1,h,i,j]>0,,h,i,j] ## get the interactions for the beginning of the simulation run: int_for_net_start=int_from_one_run[int_from_one_run[,3]<(num_ts_4_net+1),] ## get the interactions for the end of the simulation run int_for_net_end=int_from_one_run[int_from_one_run[,3]>(time-num_ts_4_net),] ## time ordered analysis: # create network for the beginning of the simulation df_for_to_start=as.data.frame(cbind(int_for_net_start,int_for_net_start[,3]+1)) names(df_for_to_start)=c('VertexFrom', 'VertexTo', 'TimeStart', 'TimeStop') # add interactions again but in reverse order to create undirected network to_add_start=cbind(df_for_to_start[,2],df_for_to_start[,c(1,3,4)]) names(to_add_start)=c('VertexFrom', 'VertexTo', 'TimeStart', 'TimeStop') df_for_to_start=rbind(df_for_to_start,to_add_start) IDS_for_to=union(df_for_to_start$VertexFrom, df_for_to_start$VertexTo) start_to_net=generatetonetwork(df_for_to_start, IDS_for_to) # create network for the end of the simulation: df_for_to_end=as.data.frame(cbind(int_for_net_end,int_for_net_end[,3]+1)) names(df_for_to_end)=c('VertexFrom', 'VertexTo', 'TimeStart', 'TimeStop') # add interactions again but in reverse order to create undirected network to_add_end=cbind(df_for_to_end[,2],df_for_to_end[,c(1,3,4)]) names(to_add_end)=c('VertexFrom', 'VertexTo', 'TimeStart', 'TimeStop') df_for_to_end=rbind(df_for_to_end,to_add_end) IDS_for_to_end=union(df_for_to_end$VertexFrom, df_for_to_end$VertexTo) end_to_net=generatetonetwork(df_for_to_end, IDS_for_to_end) ### Conduct spread analysis: sa_start <- spreadanalysis(start_to_net, seq(0,20,by=1), length(V(start_to_net)), normalizebyname=TRUE) sa_end <- spreadanalysis(end_to_net, seq(0,20,by=1), length(V(end_to_net)), normalizebyname=TRUE) b_start <- transformspreadbyindividual(sa_start) b_end <- transformspreadbyindividual(sa_end) # get proportion of individuals infected for each seed vmax_start=apply(b_start, 2,max) vmax_end=apply(b_end, 2,max) # get T50 for start and end networks t50_start=rep(NA, n) t50_end=rep(NA, n) for (k in 1:dim(b_start)[2]){ t50_start[k]=which(b_start[,k]>=0.5)[1] } for (k in 1:dim(b_end)[2]){ t50_end[k]=which(b_end[,k]>=0.5)[1] } # store variables: # start network: # when generic individuals are the seed of infection: Vmax_TO[1,h,i,j]=mean(vmax_start[2:n], na.rm=TRUE) T50_TO[1,h,i,j]=mean(t50_start[2:n], na.rm=TRUE) # when the keystone is the seed of infection: Vmax_KI_TO[1,h,i,j]=vmax_start[1] T50_KI_TO[1,h,i,j]=t50_start[1] # end network: # when generic individuals are the seed of infection: Vmax_TO[2,h,i,j]=mean(vmax_end[2:n], na.rm=TRUE) T50_TO[2,h,i,j]=mean(t50_end[2:n], na.rm=TRUE) # when the keystone is the seed of infection: Vmax_KI_TO[2,h,i,j]=vmax_end[1] T50_KI_TO[2,h,i,j]=t50_end[1] } } } ############################################## ########### Collective prey attack ########### ############################################## # Get the average boldness at the end of the simulation across all simulation iterations boldness_end_to_sort=boldness[,time,,,] boldness_end_sorted=apply(boldness_end_to_sort,c(2,3,4), sort) boldness_generic_at_end=apply(boldness_end_sorted, c(1,3,4), mean) attakProbRule=c("linear","betaDist") # set rules for a probability of an individual to attack # initialize book keeping: attackers=array(0, dim=c(length(attakProbRule),n,time_att,iter_att, length(pintcases), length(KIboldness))) # define a function that alters the probability of a spider to attack as a function of number of spiders that are already out betaDist=function(numAtt){ # use numAtt from the previous time step so is just one number vect=seq(0,1, by=0.01) y=((vect)^(2-1)*((1-vect)^(7-1)))/beta(2,7) vect=vect*n # scale so that spans the number of spiders y=y*n y=y/mean(y) # scale so mean(y)=1 y[1]=y[2]/2 # make sure first value is not 0 tmp=approx(vect,y,numAtt) # interpolation return(tmp$y) } for (g in 1:length(KIboldness)){ # keystone boldness loop for(h in 1:length(pintcases)){ # interaction rule loop pAtt=boldness_generic_at_end[,h,g] # use boldness output to determine probability to attack pAtt=(pAtt/600) # scale by largest possible boldness to get probability to attack between 0 and 1 for (k in 1:iter_att){ # iteration loop for (z in 2:time_att){ # time loop for(j in 1:length(attakProbRule)){ # attack rule loop for (i in 1:n){ # individual spider loop if(attackers[j,i,z-1,k,h,g]==1){ attackers[j,i,z,k,h,g]=1 }else{ if(attakProbRule[j]=="linear"){ if (runif(1)