Simulate stochastic bipartite network based on trait values of species - in R

Question

I would like to create bipartite networks in R. For example, if you have a data.frame of two types of species (that can only interact across species, not within species), and each species has a trait value (e.g., size of mouth in the predator allows who gets to eat which prey species), how do we simulate a network based on the traits of the species (that is, two species can only interact if their traits overlap in values for instance)?

UPDATE: Here is a minimal example of what I am trying to do. 1) create phylogenetic tree; 2) simulate traits on the phylogeny; 3) create networks based on species trait values.

# packages
install.packages(c("ape","phytools"))
library(ape); library(phytools)

# Make phylogenetic trees
tree_predator <- rcoal(10)
tree_prey <- rcoal(10)

# Simulate traits on each tree
trait_predator <- fastBM(tree_predator)
trait_prey <- fastBM(tree_prey)

# Create network of predator and prey
## This is the part I can't do yet. I want to create bipartite networks, where 
## predator and prey interact based on certain crriteria. For example, predator
## species A and prey species B only interact if their body size ratio is
## greater than X.

Answer 1

The format of the answer really depends on what you want to do next, but here's an attempt:

set.seed(101)
npred <- nprey <- 10
tree_predator <- rcoal(npred)
tree_prey <- rcoal(nprey)

## Simulate traits on each tree
trait_predator <- fastBM(tree_predator)
trait_prey <- fastBM(tree_prey)

(I used set.seed(101) for reproducibility, so these are my trait results ...

> trait_predator
         t1          t9          t4          t8          t5          t2 
-2.30933392 -3.17387148 -0.01447305 -0.01293273 -0.25483749  1.87279355 
         t6         t10          t3          t7 
 0.70646610  0.79508740  0.05293099  0.00774235 
> trait_prey
         t10           t7           t9           t6           t8           t1 
 0.849256948 -0.790261142  0.305520218 -0.182596793 -0.033589511 -0.001545289 
          t4           t5           t3           t2 
-0.312790794  0.475377720 -0.222128629 -0.095045954 

...)

The values you've generated are on an unbounded space, so it doesn't really make sense to take their ratios; we'll assume they're the logarithms of size, so exp(x-y) will be the predator/prey size ratio. Arbitrarily, I'll assume the cutoff is 1.5 ...

Use outer to compare all predator to prey traits: this creates a binary (0/1) matrix.

bmatrix <- outer(trait_predator,trait_prey,
      function(x,y) as.numeric(exp(x-y)>1.5))

One way to visualize the results ...

library(Matrix)
image(Matrix(bmatrix),xlab="Prey",ylab="Predator",sub="")

You can see for example that predator #6 (labeled t2 in the output above) is really big (log size=1.87), so it eats all the prey species ...

Using igraph (some of my approach here is a bit hacky)

library(igraph)

edges <- which(bmatrix==1,arr.ind=TRUE)  ## extract vertex numbers
## distinguish prey (columns) from pred (rows)
edges[,2] <- npred+edges[,2]             

gg <- graph.bipartite(rep(1:0,c(npred,nprey)),
             c(t(edges))) 
## c(t(edges)) collapses the two-column matrix to a vector in row order ...
## now plot ...
plot(gg,vertex.color=rep(c("cyan","pink"),c(npred,nprey)),
     edge.arrow.mode=">")

This matches the matrix format above -- predators 1 and 2 (=vertices 1 and 2) eat no-one, prey 2 (= vertex 12) is eaten by lots of different predators ... This representation is prettier but not necessarily clearer (e.g. both predator 7 and 8 eat prey 2 (vertex 12), but their arrows coincide). Having it in igraph form might be nice if you want to apply graph-theoretic approaches, though (and there are a plethora of layout options for plotting the graphs).