# 10 Models of evolution in networks

**Fernando Pedraza**

**session 31/03/2023**

## 10.1 About this session

In this session you will be implementing the coevolution model proposed in Guimaraes et al. 2017^{1}. This tutorial is split into 4 sections.

- A brief description of the coevolution model, its equations and some code to implement it in R.
- A demonstration of how to simulate the trait evolution of a pairwise interaction.
- An exercise where you will simulate the trait evolution of a small network.
- An exercise where you will simulate the trait evolution across networks with different structures.

The goal of this tutorial is for you to get an intuitive understanding of how the model works and what its predictions entail. You should gain an understanding of how changing the strength of mutualistic selection and the structure of a network shapes the coevolution of traits.

You will have to write some code to complete the exercises, but it will be mostly similar to the code I provide you in the demonstrations. If you find the code overwhelming I can walk you through it, feel free to ask questions during the exercise session or email me at fernando.pedraza@ieu.uzh.ch. If you’re up for a challenge, you can rewrite my code to make it more efficient.

## 10.2 Section 1: The coevolution model

We will be using the coevolutionary model proposed in Guimaraes et al 2017. As described during the lecture, the model incorporates coevolutionary dynamics and network structure to model trait evolution. The equation is the following:

\[ Z_i^{(t+1)} = Z_i^{(t)} + \phi_{i}\left[\sum_{j,j\neq{i}}q_{ij}^{(t)}(Z_{j}^{(t)}-Z_{i}^{(t)}) + \left(1 - \sum_{j,j\neq{i}}^Nq_{ij}^{(t)}\right)(\theta_{i}-Z_{i}^{(t)})\right]\]

where \(Z_i^{(t)}\) is the mean trait value of the local population of speices *i* at time *t*. \(\phi\) is a compound parameter that affects the slope of the selection gradient. *N* is the number of species in the network. \(q_{ij}^{(t)}\) is the interaction weight that describes the evolutionary effect of species *j* on the selection gradient of species *i* and evolves over time because of trait matching. \(\theta\) is the phenotype favoured by environmental selection. For more information on model parameters, refer to Guimaraes et al. 2017.

Note that:

\[ q_{ij}^{(t)} = m_{i}\frac{a_{ij}e^{-\alpha(Z_{j}^{(t)}-Z_{i}^{(t)})^{2}}}{\sum_{k, i\neq{k}}^{N}a_{ik}e^{-\alpha(Z_{k}^{(t)}-Z_{i}^{(t)})^{2}}}\]

where \(m_i\) is the level of mutualistic selection and \(\alpha\) is a scaling constant that controls the sensitivity of the evolutionary effect to trait matching (this is set to 0.2 as in Guimaraes et al. 2017).

Given the model, we need to build a function in R that takes the required parameters and computes trait evolution. This the function that we will be using to run the model. Notice that the function `coevolution_model`

also incorporates a parameter (`epsilon`

) that determines whether the model has reached an equilibrium. The function also has a parameter (`t_max`

) that sets the maximum number of time steps to simulate. Look over the function and try to link each line of code to operations in the model.

```
# Code for coevolutionary model
= function(f, m_value){
coevolution_model # Function to model trait evolution in a mutualistic network
# Arguments:
# f: square adjacency matrix representing the mutualistic network
# m: strength of mutualistic selection
# Returns:
# A data frame with the timesteps, trait values, species names, and theta values
#
# Set parameters that model will use
#
# Define the number of species
= dim(f)[1]
n_sp # Define alpha
= 0.2
alpha # Define the distribution for phi
= 0.7
phi_mean = 0.01
phi_sd # Draw phi value for species
set.seed(1)
= rnorm(n_sp,phi_mean,phi_sd)
phi # Define mutualism selection distribution
= m_value
m_mean = 0.01
m_sd # Draw mutualism selection value for species
set.seed(2)
= rnorm(n_sp,m_mean,m_sd)
m # Define range of theta values
= 0
theta_min = 10
theta_max # Draw theta value for species
set.seed(3)
= runif(n_sp,min = theta_min,max = theta_max)
theta # Draw initial trait value for species
set.seed(4)
= runif(n_sp, min = theta_min, max = theta_max)
init # Define threshold to stop simulations
= 0.000001
epsilon # Define the max number of generations to simulate
= 10000
t_max
# Initialize an empty matrix to store trait values.
# The number of columns of the matrix is equal to the number of
# species in the initial matrix
# The number of rows of the matrix is equal to the maximum number of
#simulation steps
= matrix(NA,nrow=t_max,ncol=n_sp)
trait_sim_mat colnames(trait_sim_mat) <- paste0('Species ', 1:n_sp)
# Store initial trait values in first row of result matrix
1,] = init
trait_sim_mat[
# Perform the simulation for a maximum of t_max-1 time steps
for (time_step in 1:(t_max - 1)){
# Define the trait values for the current simulation
= trait_sim_mat[time_step,]
current_trait_values # Compute the differences in trait values for all possible species
# combinations
# Results are stored as matrix
= t(f*current_trait_values) - f*current_trait_values
trait_differences_mat # Compute the q matrix
# Values in this matrix represent the evolutionary effect that
# species have on each other
= f*(exp(-alpha*(trait_differences_mat^2)))
q # Standardize the matrix to a maximum value of 1
= q / apply(q,1,sum)
q_n # Compute selection by mutualism
= q_n * m
q_m # Calculate selection differentials
= q_m * trait_differences_mat
sel_dif # Calculate trait change as a result of mutualism
= phi * apply(sel_dif, 1, sum)
r_mut # Calculate trait chagne as a result of environment
= phi * (1 - m)* (theta - current_trait_values)
r_env # Calculate trait change given mutualism and environment and update
# matrix with results from simulation
+ 1,] = current_trait_values + r_mut + r_env # updating z values
trait_sim_mat[time_step # Calculate how much trait changed from previous time step to
# current one
= mean(abs(current_trait_values - trait_sim_mat[time_step + 1,]))
dif # If the the observed change in trait value is smaller than
# epsilon, then exit loop
if(dif < epsilon)
break
}
# get the output
<- trait_sim_mat[1:(time_step+1),]
matrix_output # generate time steps
<- rep(1:nrow(matrix_output), times = n_sp)
time_step_vector # extract theta values
<- rep(theta, each = nrow(matrix_output))
theta_vector # convert to dataframe
<- as.data.frame(matrix_output)
dataframe_output # add column
<- cbind(time_step_vector, melt(dataframe_output))
dataframe_output # add the theta column
$theta <- theta_vector
dataframe_output# rename columns
colnames(dataframe_output) <- c('time', 'species', 'trait_value', 'theta')
# Return a matrix containing the evolution of traits
return(dataframe_output)
}
```

After defining our function, lets run it under different scenarios.

## 10.3 Section 2: Running the coevolution on model on a pairwise interaction

Now we run `coevolution_model`

with two interacting species. In this section we will run the simulation with progressively higher \(m\) values. We want to know how the strength of mutualistic selection alters the coevolution of traits.

Before we start, we load the packages we will be using for this exercise session:

```
# Load required packages
library(ggplot2)
library(ggthemes)
library(igraph)
library(reshape2)
```

### 10.3.2 Run simulations

First we define the adjacency matrix we will be using for the simulations.

```
# Define the adjacency matrix of the community
<- matrix(c(0,1,1,0),nrow=2,ncol=2) f
```

Next we run the simulations with weak mutualistic selection

```
# Define the strength of mutualistic selection
<- 0.3
weak_m # Run simulations for weak m
= coevolution_model(f,weak_m)
coevolution_output_weak_m # Add the m value to our data frame
$m <- weak_m coevolution_output_weak_m
```

Now we move on to the scenario with intermediate mutualistic selection

```
# Define the strength of mutualistic selection
<- 0.6
mid_m # Run simulations for mid m
= coevolution_model(f,mid_m)
coevolution_output_mid_m # Add the m value to our data frame
$m <- mid_m coevolution_output_mid_m
```

Finally, we run the scenario with strong mutualistic selection

```
# Define the strength of mutualistic selection
<- 0.9
strong_m # Run simulations for strong m
= coevolution_model(f,strong_m)
coevolution_output_strong_m # Add the m value to our data frame
$m <- strong_m coevolution_output_strong_m
```

Next, we put all our results together and visualize them. In the plot I show how the trait value of each species changed through time. I show the environmental optimum of each species as dashed lines on the plot.

```
# Now we put all the results together
<- rbind(coevolution_output_weak_m,
coevolution_output
coevolution_output_mid_m,
coevolution_output_strong_m)
# Plot result
ggplot(data = coevolution_output,aes(x=as.numeric(time),y=trait_value,col=species)) +
geom_line() + geom_point(size=0.8) + theme_minimal() +
geom_hline(aes(yintercept=theta, color = species),linetype="dashed") +
xlab("Time") + ylab("Trait Value") + theme(legend.position = 'bottom') +
facet_grid(~m)
```

## 10.4 Section 3: Running the coevolution on model on a small network

Now you will run `coevolution_model`

with four interacting species. Again, we will run the model with progressively stronger mutualistic selection. We want to know how the strength of mutualistic selection alters the coevolution of traits.

### 10.4.2 Your task

You will now simulate coevolution on the network shown above under three scenarios of mutualistic selection: 0.3, 0.6, and 0.9. You should produce a plot similar to the one we produced when simulating coevolution for two species.

```
# Work on you Rscript file to complete the task. You should do the following tasks:
# 1. Build the adjacency matrix
# 2. Simulate coevolution for the three different m values (remember to add the m value to your output!)
# 3. Combine the output from the three runs into a single dataframe
# 4. Plot the results
```

Once you’ve produced the plot, answer (write the answers in your Rscript) the following questions:

**How does the strength of mutualistic selection affect the dynamics of trait evolution?**

**What differences can you observe in the trait evolution of species when going from a pairwise interactions to a (small) network?**

## 10.5 Section 4: Simulating coevolution across networks with different structures

Now we are interested in exploring how network structure influences the outcome of coevolution of entire communities. To do this, we will simulate coevolution on a set of networks of identical size but with different connectance. I am providing you with these networks. Note that I have generated the networks artificially, by means of a model. Even though these are not empirical networks, they will allow us to explore our question of interest in this toy example.

First we load the networks.

```
# The networks are stored inside a list saved as an Rdata file.
load('~/ecological_networks_2023/downloads/Data/03-31_models_of_evolution_in_networks/list_of_networks.Rdata')
```

Next we define two new functions. The first is a slightly modified version of the coevolution model. This alternative version simulates coevolution on a given function and calculates the average trait matching between all species in the community. Larger values indicate higher trait similarity between species. This function also calculates the connectance of a network.

```
# Code for coevolutionary model
= function(f, m_value){
coevolution_model_across_networks # Function to model trait evolution in a mutualistic network
# Arguments:
# f: square adjacency matrix representing the mutualistic network
# m: strength of mutualistic selection
# Returns:
# A dataframe with the timesteps, trait values, species names, and theta values
#
# Set parameters that model will use
#
# convert incidence matrix to adjacency
<- as.matrix(as_adjacency_matrix(graph_from_incidence_matrix(as.matrix(f))))
f # Define the number of species
= dim(f)[1]
n_sp # Define alpha
= 0.2
alpha # Define the distribution for phi
= 0.7
phi_mean = 0.01
phi_sd # Draw phi value for species
set.seed(1)
= rnorm(n_sp,phi_mean,phi_sd)
phi # Define mutualism selection distribution
= m_value
m_mean = 0.01
m_sd # Draw mutualism selection value for species
set.seed(2)
= rnorm(n_sp,m_mean,m_sd)
m # Define range of theta values
= 0
theta_min = 10
theta_max # Draw theta value for species
set.seed(10)
= runif(n_sp,min = theta_min,max = theta_max)
theta # Draw initial trait value for species
set.seed(11)
= runif(n_sp, min = theta_min, max = theta_max)
init # Define threshold to stop simulations
= 0.000001
epsilon # Define the max number of generations to simulate
= 10000
t_max
# calculate the connectance of the current matrix
<- graph.density(graph_from_adjacency_matrix(f))
network_connectance
# Initialize an empty matrix to store trait values.
# The number of columns of the matrix is equal to the number of
# species in the initial matrix
# The number of rows of the matrix is equal to the maximum number of
#simulation steps
= matrix(NA,nrow=t_max,ncol=n_sp)
trait_sim_mat colnames(trait_sim_mat) <- paste0('Species ', 1:n_sp)
# Store initial trait values in first row of result matrix
1,] = init
trait_sim_mat[
# Perform the simulation for a maximum of t_max-1 time steps
for (time_step in 1:(t_max - 1)){
# Define the trait values for the current simulation
= trait_sim_mat[time_step,]
current_trait_values # Compute the differences in trait values for all possible species
# combinations
# Results are stored as matrix
= t(f*current_trait_values) - f*current_trait_values
trait_differences_mat # Compute the q matrix
# Values in this matrix represent the evolutionary effect that
# species have on each other
= f*(exp(-alpha*(trait_differences_mat^2)))
q # Standardize the matrix to a maximum value of 1
= q / apply(q,1,sum)
q_n # Compute selection by mutualism
= q_n * m
q_m # Calculate selection differentials
= q_m * trait_differences_mat
sel_dif # Calculate trait change as a result of mutualism
= phi * apply(sel_dif, 1, sum)
r_mut # Calculate trait chagne as a result of environment
= phi * (1 - m)* (theta - current_trait_values)
r_env # Calculate trait change given mutualism and environment and update
# matrix with results from simulation
+ 1,] = current_trait_values + r_mut + r_env # updating z values
trait_sim_mat[time_step # Calculate how much trait changed from previous time step to
# current one
= mean(abs(current_trait_values - trait_sim_mat[time_step + 1,]))
dif # If the the observed change in trait value is smaller than
# epsilon, then exit loop
if(dif < epsilon)
break
}
# extract the trait values at equilibrium
<- trait_sim_mat[time_step + 1,]
traits_at_equilibrium # compute the degree of trait matching in the community
<- as.matrix(dist(traits_at_equilibrium,upper=TRUE,diag=TRUE))
matching_diff <- (exp(-alpha*((matching_diff)^2)))
matching_at_equilibrium_all <- mean(rowMeans(matching_at_equilibrium_all,na.rm = TRUE))
matching_at_equilibrium # return a vector with two elements:
# The first is the connectance of the network.
# The second is the euclidean distance of the traits in the community
return(c(network_connectance,matching_at_equilibrium))
}
```

The second function allows us to run the coevolution model across all the networks in our list. As parameters, the function requires the list of networks we want to use and the strength of mutualistic selection we want the model to run. The function returns a vector with two elements. The first column contains the connectance values of each network. The second column contains the average trait matching in the community.

```
<- function(list_of_networks, m_value){
simulate_coevolution_across_networks
# create empty dataframe
<- NULL
df # loop through the list of networks
for (i in 1:length(list_of_networks)) {
# run the coevolution model for the current network with the specified mutualistic selection value
<- coevolution_model_across_networks(list_of_networks[[i]], m_value = m_value)
current_coevolution_result # append the results from the simulation to the dataframe
<- rbind(df, current_coevolution_result)
df
}# make sure that the dataframe is in the correct format
<- as.data.frame(df)
df # rename the dataframe columns
colnames(df) <- c('connectance', 'trait_matching')
# return the dataframe
return(df)
}
```

### 10.5.1 Your task

We are interested in seeing how connectance is related to the trait matching in a community after coevolution. Essentially, we want to know how network structure affects the similarity of traits in a community after they have coevolved. We also want to know if the relationship between network structure and trait matching changes as we increase the strength of mutualistic selection. To do this you will run `simulate_coevolution_across_networks`

with `list_of_networks`

twice, one with \(m = 0.1\) and the second with \(m = 0.9\).

After running the models you should produce two plots, one for \(m = 0.1\) and another for \(m = 0.9\). In both you should show the relationship between network connectance and trait matching. To help you to analyse the trend, fit a simple linear model to the plot (you can add a linear model to a ggplot using the following command: `ggplot(...) + geom_smooth(method = 'lm')`

).

Once you have produced the plots, answer the following questions:

**How does the connectance of a network relate to the similarity of traits in a community when coevolution is weak? What happens when coevolution is strong?**

**Can you come up with a hypothesis as to why we observe these patterns? What more simulations could you run to see if your hypothesis is true? (There is no wrong answer here, it’s perfectly fine to speculate!) **

Guimarães, P., Pires, M., Jordano, P. et al. Indirect effects drive coevolution in mutualistic networks.

**Nature***550*, 511–514 (2017).↩︎