I have been reading a little about neutral dynamics in food web models (posted here). So while I was ~~procrastinating~~ working hard, I thought that I would try and write up some R code to generate static food webs using neutral processes.

I am still new to this, so it is not the prettiest code but it works (sort of). The central idea of what I’ve done here is that the probability of any two species interacting is going to be a function of population density. So the inputs of this function are just * n* which is the number of iterations, and

**which is the number of species to put into the web, meaning that the only effective parameter is**

*S**. Population densities for species 1 through*

**S****are drawn from a random lognormal distribution. Then a vector of probabilities is computed by dividing by the largest density (to get probabilities between 0 and 1).**

*S*Using the vector of probabilities it populates an SxS matrix where the probability of interaction is the sum of the probabilities of each species. Then that matrix is transformed so that the interaction probabilities are between 0 and 1. The transformed matrix provides the probabilities of drawing a 1 from a binomial distribution. Thus a binary adjacency matrix is created where a 1 means that species *i *interacts with species *j*.

The function then can spit out any number of indices to describe the network. Here I have shown an output matrix with the connectance (*C*), clustering coefficient, and nestedness (measured by NODF values).

The problem, however, is that while the model produces a range of connectance values that are marginally close to actual values they tend to be larger relative to empirically determined connectances. Also, I don’t believe it is exceedingly realistic, and there is really no justification for the probability matrix or the way probability is determined. Basically I just toyed around with the math until something fairly acceptable was produced.

I put the code at the bottom of the post, but I would be happy to hear about any thoughts people may have about generating a static neutral network. I think that being able to generate networks using a neutral model would be invaluable to comparisons of networks generated by static models such as the niche or nested heirarchy model. This does seem to be something other people are interested in as well, I posted here about a paper on neutral network models based on population dynamics models.

Neutral.Network.Model<-function(n,S){

output<-matrix(nrow=n,ncol=3,byrow=T)

for(k in 1:n){

probability.matrix<-matrix(nrow=S,ncol=S)

pop.size<-rlnorm(S)

prob<-pop.size/max(pop.size)

for(i in 1:S){

for(j in 1:S){

probability.matrix[i,j]<-prob[i]+prob[j] # add prob instead of multiply

}

}

transform.matrix<-(probability.matrix)/max(probability.matrix)

neutral.int.matrix<-matrix(nrow=S,ncol=S)

for(i in 1:S){

for(j in 1:S){

neutral.int.matrix[i,j]<-rbinom(1,1,transform.matrix[i,j])

}

}

neutral.graph<-graph.adjacency(neutral.int.matrix,mode=”directed”)

Ind<-GenInd(Tij=t(neutral.int.matrix))

output[k,1]<-Ind$C

output[k,2]<-clusters(neutral.graph,mode=c(“strong”))$no

output[k,3]<-nested(neutral.int.matrix,method=”NODF”)

}

return(output)

}