# 7 Simulating networks

**Alessandro Vindigni and Fernando Pedraza**

**session 28/03/2023**

## 7.1 Models for simulating networks

Theoretical models used for generating networks can be useful for at least two reasons. First, they can help us understand the underlying processes driving the structure of ecological networks. Second, if the models perform well, they can help us expand the number of networks we work with.

Several theoretical models exist for generating artificial networks. Most were originally proposed for food webs and later were adapted for mutualistic interactions. Here we will explore how two different models (the cascade and the niche models) proposed by Pires et al. 2011 perform in reproducing the structural properties of a set of mutualistic networks. Next we provide a brief description of the assumptions of each model. When describing the models, for simplicity, we will refer to one guild as the “animals” and the other as “plants”. Both models establish interactions based on the traits of species. These traits are numbers that are randomly assigned to each animal and plant species in the community.

### 7.1.1 The cascade model

The cascade model assumes that interactions are based on hierarchical feeding. In other words, a species can interact with species whose trait values are smaller than their own. The algorithm for the model is the following:

Assign a random trait value drawn from a uniform distribution on the interval [0,1] to each animal and plant in the community. For mutualistic networks, these trait values are drawn independently for each guild.

Compare the traits of all animals and plants to determine whether or not they interact. Animals can potentially interact with plants whose trait values are smaller than their own but can never interact with other animals.

Determine the probability of interactions occurring: \(q = E/T\), where \(E\) is the number of recorded interactions and \(T\) is the number of possible interactions after species positions are defined.

### 7.1.2 The niche model

The niche model also establishes interactions based on species’ traits. Yet, it assumes that animals can only interact with plants whose traits fall within a range. The algorithm for the model is the following:

Assign a random trait value drawn from a uniform distribution on the interval [0,1] to each animal and plant in the community. For mutualistic networks, these trait values are drawn independently for each guild.

For each consumer (\(i\)) calculate it’s niche range(\(r_i\)). \(r_i=xn_i\), where \(0 ≤ x ≤ 1\) is a random variable with a beta-distributed probability density function \(p(x) = \beta(1-x)(\beta-1)\) with \(\beta=(1/2C) - 1\).

Determine the range center (\(c_i\)) for each species. \(c_i\) is a uniformly random number between \(r_i/2\) and \(min(n_i, 1 - r_i/2)\).

Determine the interaction interval (\(I(D_i)\)) for each species. \(I(D_i) = [c_i - r_i/2, c_i + r_i/2]\).

Establish interactions between plants and animals. An animal \(i\) will interact with all plants \(j\) whose \(n_j\) falls within its interaction interval.

## 7.2 Running the models for one community

We are providing you the code to implement the models. However, because
simulating a network can be a rather slow process and R is already
plenty slow(!), we coded the functions in the Julia language (Julia is a
relatively new language that merges readability with speed). The code is
stored on the R server, in `ecological_networks_2023/downloads/Functions/build_networks.jl`

.
However, to keep things manageable, we will be running the Julia code
from R. To do this, we will use the `JuliaCall`

package.

To initialise a Julia session from R, run the following commands:

```
# load the package
library(JuliaCall)
# if the package is not installed, run the following command:
# install.packages(JuliaCall)
# Initialise a Julia session from R
julia_setup("/home/ubuntu/julia-1.8.1/bin/")
# Then install some required packages (note you only have to do this once!)
julia_install_package("Distributions")
julia_install_package("Random")
# Source the functions that we will use to run the models (change test_student with your user name)
julia_source("/home/test_student/ecological_networks_2023/downloads/Functions/build_networks.jl")
```

Now we can run the model to generate networks. The function to run the
cascade model is `build_cascade_networks`

while `build_niche_networks`

runs the niche model. Both functions require three arguments:

- the number of animal species in the community.
- the number of plant species in the community.
- the target the connectance (i.e. the connectance that the artificial network should have).

The function returns an incidence matrix that corresponds to the requested artificial network. For example, to generate an artificial incidence matrix composed of 20 animals and 30 plants with a target connectance of 0.2 using the niche model, you would run the following command:

```
# run the julia code with the specified parameters and store the incidence matrix.
<- julia_eval("build_niche_networks(20, 30, 0.2)") niche_artificial
```

You can then analyse the artificial incidence matrix like you have done so in the past days. For example, you can plot it or calculate some structural properties (e.g. nestedness or modularity).

## 7.3 Compare artificial networks against empirical ones

Using the niche and the cascade models, we have created a set of artificial networks and stored them in a file. To load that data run the following commands:

```
# this block has to occur in sequence
# clean up
rm(list=ls())
# load simulated networks
load("~/ecological_networks_2023/downloads/Data/03-28_simulating_networks/artificial_matrices.Rda")
<- mget(ls()) artificial_networks
```

As you can check, the variable `artificial_networks`

is a named list with elements

```
<- names(artificial_networks)
simulated_nw_names simulated_nw_names
```

```
## [1] "M_PL_008_cascade_inc" "M_PL_008_niche_inc" "M_PL_025_cascade_inc"
## [4] "M_PL_025_niche_inc" "M_PL_037_cascade_inc" "M_PL_037_niche_inc"
## [7] "M_PL_038_cascade_inc" "M_PL_038_niche_inc" "M_PL_059_cascade_inc"
## [10] "M_PL_059_niche_inc" "M_SD_014_cascade_inc" "M_SD_014_niche_inc"
## [13] "M_SD_021_cascade_inc" "M_SD_021_niche_inc" "M_SD_024_cascade_inc"
## [16] "M_SD_024_niche_inc" "M_SD_025_cascade_inc" "M_SD_025_niche_inc"
## [19] "M_SD_029_cascade_inc" "M_SD_029_niche_inc"
```

each element containing the incidence matrix of a simulated network. As you can read from the list above, the model used in the simulation is defined in the name of each element. The initial part of the name refers, instead, to the corresponding empirical network of web of life.

Before proceeding, it is useful to load the following packages:

```
library(igraph)
library(rweboflife)
library(formattable)
```

We have also pre-computed the connectance, nestedness, and modularity of the 10 considered empirical networks. You can load that data as well from file:

```
# load indices computed for empirical networks
load("~/ecological_networks_2023/downloads/Data/03-28_simulating_networks/empirical_nw_df.Rda")
%>% formattable() empirical_nw_df
```

interaction_type | network_name | nestedness | connectance | modularity |
---|---|---|---|---|

Pollination | M_PL_008 | 0.3978981 | 0.2535885 | 0.3428711 |

Pollination | M_PL_025 | 0.5615351 | 0.2500000 | 0.2825811 |

Pollination | M_PL_037 | 0.2866498 | 0.1800000 | 0.4926698 |

Pollination | M_PL_038 | 0.3549508 | 0.2351190 | 0.4081077 |

Pollination | M_PL_059 | 0.9346739 | 0.4201183 | 0.1748661 |

SeedDispersal | M_SD_014 | 0.8245156 | 0.4448529 | 0.1875213 |

SeedDispersal | M_SD_021 | 0.6623913 | 0.2559524 | 0.2577670 |

SeedDispersal | M_SD_024 | 0.6318829 | 0.4761905 | 0.2321875 |

SeedDispersal | M_SD_029 | 1.0000000 | 0.5000000 | 0.2650000 |

SeedDispersal | M_SD_025 | 0.5714092 | 0.3307692 | 0.2985398 |

Your goal is to produce a similar dataset using the simulated networks. To help you in this task, we will guide you through the computation of nestedness, modularity, and connectance for just one artificial network.

Assume you want to focus on the 4th artificial network of the list named

```
<- 4 # 4th simulated network
nw_label simulated_nw_names[nw_label]
```

`## [1] "M_PL_025_niche_inc"`

You can extract the incidence matrix as follows

```
# extract incidence matrix form the list
<- artificial_networks[[nw_label]] inc_matrix
```

and compute the relative nestedness using the `rweboflife`

package

```
# compute netsedness
<- rweboflife::nestedness(inc_matrix)
nest nest
```

`## [1] 0.4732627`

and the relative modularity (following the procedure discussed in a previous exercise session)

```
# compute modularity
<- graph_from_incidence_matrix(inc_matrix)
nw_graph <- cluster_fast_greedy(nw_graph)
modules <- modularity(modules)
mm mm
```

`## [1] 0.3818828`

If you wish to compute the connectance from the incidence matrix (instead of using the edge list as in the first exercise session), we provide you a helping function:

```
<- function(M){
connectance_from_inc_matrix # Make sure we are working with a matrix
<- as.matrix(M)
M # Binarize the matrix
<- as.matrix((M>0))
mat class(mat) <- "numeric"
# extract dimensions
<- nrow(mat)
num_rows <- ncol(mat)
num_cols
# obtain number of links
<- sum(mat)
links
# compute connectance
<- links/(num_rows * num_cols)
con
# return
return(con)
}
```

Then you can apply this function directly to the incidence matrix of simulated networks:

```
<- connectance_from_inc_matrix(inc_matrix)
conn conn
```

`## [1] 0.2727273`

### Exercise 2

Using the data and the commands provided above, produce a pair of plots (box plot, scatter plot, etc.) that highlight the differences between the observed and the predicted values of modularity and nestedness.

Comment on the differences between the values predicted using either the cascade or the niche model.