**Introduction**

The dynamics and composition of natural communities are ultimately determined by how individuals interact with their biotic and abiotic environment via demographic rates (survival and reproduction) in order to persist in space and time. Here, we develop hierarchical functions to scale from stage-specific demographic rates to community dynamics. The functions first describe stage-specific demographic rates for various species that are affected simultaneously by environmental factors and biotic interactions across different habitats. These rates then determine transitions among different life-cycle stages in a local populations and dispersal among populations, which, in turn, characterizes the metapopulation dynamics of all species in a community.

The metapopulation framework implemented in these functions can model the dynamics of any number of species present in any number of sites, with the only constraint that we assume three distinct life stages, which generally correspond to juvenile individuals, non-reproductive adults, and reproductive adults. The methodology follows the vec-permutation approach introduced in Hunter & Caswell (2005) for projecting metapopulation abundances.

The vital rates considered in this framework include the survival probability of each life stage, the transitions between stages, reproductive outputs, and dispersal. Specifically:

- Sj: Survival probability of juveniles
- Sn: Survival probability of non-reproductive adults
- Sr: Survival probability of reproductive adults
- Rn: Probability of a non-reproductive adult to transition to the reproductive stage
- Rr: Probability of a reproductive adult to remain in the reproductive stage
- D: Probability of dispersal (emigration)
- Ds: Probability of surviving all stages of dispersal and establishing in a new site (immigration)
- O: Recrutiment, i.e., per capita number of offspring that survive to the juvenile stage produced (≥ 0)

These vital rates can be parameterized to differ weakly or strongly between sites and species. They arise from the combined effect of species-specific values (given by the density of the focal species, \(N_{s1}\) below), the effect of individuals of other species present in the same site (\(N_{s2,...,sn}\)), and potentially the effect of environmental variables (\(E\)). In particular, for a given vital rate \(v\) of a species \(s1\), they are computed as:

\[Logit(v_{s1}) = \alpha + \beta_1*N_{s1} + \sum_{i \in {s2,...,sn}} \beta_i*Ni + \beta_E*E\]

In this formulation, interactions between species densities and the environment are not shown for readability, but these can be included as well.

**Data structures**

The vec-permutation approach relies on a block-diagonal formulation of the demographic and dispersal processes. For each species, demographic changes are described by B, a block-diagonal matrix where the \(i^{th}\) diagonal block \(K_i\) is a 3x3 matrix population model at site \(i\). Dispersal is described by \(M\), another block-diagonal matrix where the \(j^{th}\) diagonal block \(L_j\) is a \(i*i\) matrix of disperal probabilities for each stage \(j\). In our case, non-reproductive adults disperse to become reproductive adults in a different site, and we therefore populate the corresponding off-diagonal element of \(M\). Our metapopulation model assumes that local demography happens first and then individuals disperse.

To run our approach of metapopulation dynamics of interacting species, we need to first define a series of constant values and data structures. Note that we explicitly consider environmental forcing: currently, our framework accepts a single environmental factor that varies throughout the number of projected timesteps.

```
# load the package
library(cxr)
# define species names
sp <- c("s1","s2","s3")
# number of species
num.sp <- length(sp)
# define site names
sites <- c("sa","sb")
# number of sites
num.sites <- length(sites)
# number of demographic stages - this should be always fixed
num.stages <- 3
# vital rates to account for - these names should be fixed
rates <- c("Sj","Sn","Sr","Rn","Rr","D","Ds","O")
# years (or time steps) of simulations
years <- 100
# simulate environment from normal distribution
set.seed(123)
env <- rnorm(years, mean=0, sd=1)
```

At this stage, we need to generate or import the coefficients for
obtaining vital rates for every species in every site. This data
structure thus holds the \(\alpha\) and
\(\beta\) coefficients from Eq. 1. In
this example, we import these coefficients in an already created list
named `metapopulation_example_param`

. We suggest to separate
in different scripts the generation of the data structure and the
projection of the dynamics. Nevertheless, at the end of this vignette we
explain in more detail the different options to generate or import these
coefficients.

```
data("metapopulation_example_param", package = "cxr")
# coefficients for species "s1" in site "sa"
metapopulation_example_param[["s1"]][["sa"]]
```

```
## alpha beta1 beta2 beta3 beta4 beta5
## Sj 1.0388308 1.3457376 0.069977335 -0.078802049 -0.0599899328 0.002809132
## Sn 3.0471897 1.5928148 -0.050277967 -0.044215025 -0.0214356158 -0.012578936
## Sr 3.8510474 2.3649396 -0.037753281 -0.081006020 -0.0408717604 -0.001894939
## Rn 1.2175032 0.0000000 -0.089419233 0.000000000 0.0000000000 0.000000000
## Rr 0.9076192 0.1785048 -0.022440715 -0.029090383 0.0039708918 0.000000000
## D -0.3153730 0.1924553 0.007495063 -0.007089005 0.0001276958 0.000000000
## Ds 0.7000000 0.1200000 0.000000000 0.000000000 0.0000000000 0.000000000
## O 2.6156340 0.7026241 -0.081252277 -0.061627015 -0.0826407224 0.001873481
## beta6 beta7
## Sj 0.004493186 0.0006890698
## Sn -0.020976595 0.0187594251
## Sr -0.052637814 -0.0333495434
## Rn 0.000000000 0.0000000000
## Rr 0.000000000 0.0000000000
## D 0.000000000 0.0000000000
## Ds 0.000000000 0.0000000000
## O -0.026285279 0.0000000000
```

Next, we define the structure that will hold the vec-permutation matrices, with the given number of species, sites, and stages. There are three types of matrices in the vec-permutation methodology: permutation matrices, dispersal, and demography. There is one matrix of each type for each species, and these matrices all have dimensions of number of stages multiplied by number of sites in both rows and columns.

```
# build the templates for the vec-permutation matrices
# returns a nested list: [[matrix.type]][[species number]]
# where matrix.type is "demography", "dispersal", or "permutation"
vpm <- vec_permutation_matrices(num.sp,num.sites,num.stages)
# for example
vpm[["demography"]][[1]]
```

```
## [,1] [,2] [,3] [,4] [,5] [,6]
## [1,] 0 0 0 0 0 0
## [2,] 0 0 0 0 0 0
## [3,] 0 0 0 0 0 0
## [4,] 0 0 0 0 0 0
## [5,] 0 0 0 0 0 0
## [6,] 0 0 0 0 0 0
```

We also need to define manually the initial abundances of each species at each stage in each site. These are stored in a list in which each element holds the abundance of a given species in a matrix of dimensions number of sites in rows and life stages in columns.

```
# initial.densities: list: [[species]][sites*stages]
# this needs to be filled manually
initial.densities <- list()
# sp1
initial.densities[[1]] <- matrix(c(10,10,15,10,10,13),
nrow = num.sites,
ncol = num.stages,
byrow = TRUE)
# sp2
initial.densities[[2]] <- matrix(c(15,5,3,15,5,3),
nrow = num.sites,
ncol = num.stages,
byrow = TRUE)
# sp2
initial.densities[[3]] <- matrix(c(5,5,2,3,4,2),
nrow = num.sites,
ncol = num.stages,
byrow = TRUE)
```

For example, the initial abundances of the first species are as follows:

```
## [,1] [,2] [,3]
## [1,] 10 10 15
## [2,] 10 10 13
```

Where there are 10 juvenile individuals at site 1, 15 reproductive adults at site 1, and so on.

Lastly, we need to define the data structure that will store the whole projected dynamics, i.e. species densities in each timestep, site, and stage. This is, again, a nested list with species in the first level, year (or more generally, timestep) in the second level, and each element of the third level being a matrix of dimensions of number of sites (rows) by number of stages (columns).

```
# projected.densities list: [[species]][[years]][[sites*stages]
projected.densities <- list()
for(i.sp in 1:num.sp){
projected.densities[[i.sp]] <- list()
for(i.year in 1:years){
projected.densities[[i.sp]][[i.year]] <- matrix(0,
nrow = num.sites,
ncol = num.stages)
}
}
```

**Projecting population densities**

The densities of each species at each timestep are calculated within
a loop. In this loop, first we update the transition matrix of each
species, given their prior densities and coefficients. This is done with
the function `fill_transition_matrix`

. Second, we update the
vec-permutation demography and dispersal matrices of each species, with
the functions `fill_demography_matrix`

and
`fill_dispersal_matrix`

. With these two matrices, we can
finally obtain the expected density of each species at the next timestep
using the function `calculate_densities`

.

```
# create an auxiliary list, to keep track of the densities per timestep
current.densities <- initial.densities
for(i.year in 1:years){
# Update projected densities at this timestep
for(i.sp in 1:num.sp){
projected.densities[[i.sp]][[i.year]] <- current.densities[[i.sp]]
}
# update transition matrices ----------------------------------------------
# this is a list per species and site
transition_matrices <- list()
for(i.sp in 1:length(sp)){
transition_matrices[[i.sp]] <- list()
for(i.site in 1:length(sites)){
# store the transition matrix for this sp and site
transition_matrices[[i.sp]][[i.site]] <- fill_transition_matrix(focal.sp = i.sp,
site = i.site,
param = metapopulation_example_param,
env = env[i.year],
current.densities = current.densities)
}# for each site
names(transition_matrices[[i.sp]]) <- sites
}# for each species
names(transition_matrices) <- sp
# update demography and dispersal matrices --------------------------------
for(i.sp in 1:length(sp)){
# demography
vpm[["demography"]][[i.sp]] <- fill_demography_matrix(focal.sp = i.sp,
vpm = vpm,
transition_matrices = transition_matrices)
# dispersal
vpm[["dispersal"]][[i.sp]] <- fill_dispersal_matrix(focal.sp = i.sp,
num.sites = num.sites,
param = metapopulation_example_param,
vpm = vpm,
env = env[i.year],
current.densities = current.densities)
}# for i.sp
# update densities --------------------------------------------------------
# all stages and sites for each species
for(i.sp in 1:length(sp)){
current.densities[[i.sp]] <- calculate_densities(focal.sp = i.sp,
vpm,
current.densities)
}# for i.sp
}# for i.year
```

Finally, we can plot the resulting dynamics over time, species, and/or sites

```
# transform list of projected densities to dataframe
df <- densities_to_df(projected.densities)
# tidy
df$species <- dplyr::recode(as.character(df$species), "1" = "S1", "2" = "S2", "3" = "S3")
# plot
dynamics.plot <- ggplot2::ggplot(df,ggplot2::aes(year,density,col=species))+
ggplot2::geom_line()+
ggplot2::facet_grid(stage~site,scales = "free")+
ggplot2::scale_color_manual(name="",values=c("darkgreen","orange","darkred"))+
ggplot2::xlab("Simulation year")+ggplot2::ylab("Total density")+
ggplot2::theme_bw(base_size=20)+
ggplot2::theme(panel.grid = ggplot2::element_blank())+
ggplot2::theme(panel.grid.major = ggplot2::element_blank(),
panel.grid.minor = ggplot2::element_blank(),
strip.background = ggplot2::element_blank(),
panel.border = ggplot2::element_rect(colour = "black"))
```

**Appendix: Different options for building the coefficient data
structure**

Users can generate the list with parameter coefficients through a set
of helper functions. The data structure is created with the function
`build_param`

, and populated with the function
`generate_vital_rate_coefs`

.

The function `build_param`

accepts the name of the species
modelled, the number or name of sites, a string giving the vital rates
included, and an argument `env`

specifying whether
environmental forcing is to be accounted for. Lastly, the coefficients
from Eq. 1 may include the interaction between the density of one or
more species and the environment. The user has the possibility to limit
the number of env:species interactions by setting the argument
`num.params`

. For example, given three species and including
environmental forcing, there are a minimum of 1+1+3 coefficients (the
intercept, one coefficient for the environmental effect, and one
coefficient for the effect of each species), and a maximum of 1+1+3+3
coefficients, by including each interaction. In this example, we include
all environment:density interactions, for a total of eight coefficients.
It is assumed that the interactions included follow the order in which
the species are entered.

```
# this builds an empty data structure
example_param <- build_param(sp = sp,
sites = sites,
rates = rates,
env = env,
num.params = 8)
```

The `example_param`

structure is a nested list, sorted by
species in the first level and site in the second level. Each element is
a matrix of dimensions number of vital rates by number of
coefficients.

```
## alpha beta1 beta2 beta3 beta4 beta5 beta6 beta7
## Sj NA NA NA NA NA NA NA NA
## Sn NA NA NA NA NA NA NA NA
## Sr NA NA NA NA NA NA NA NA
## Rn NA NA NA NA NA NA NA NA
## Rr NA NA NA NA NA NA NA NA
## D NA NA NA NA NA NA NA NA
## Ds NA NA NA NA NA NA NA NA
## O NA NA NA NA NA NA NA NA
```

There are a number of ways in which these matrices can be populated.
The elements, as shown by Eq. 1, are equivalent to the coefficients (on
the link scale) of a Generalized Linear Model (GLM). Therefore, if the
user has a tidy table with the same number of coefficients, this table
can be input to the `generate_vital_rate_coefs`

function via
the `glm.object`

argument. For example, using an example
coefficient table from a GLM:

```
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.0388308163 0.380762792 2.7282887 6.366387e-03
## dens 0.0699773352 0.010531993 6.6442633 3.047375e-11
## env 1.3457375997 0.353523654 3.8066409 1.408671e-04
## densS2 -0.0788020486 0.014226983 -5.5389150 3.043514e-08
## densS3 -0.0599899328 0.027802360 -2.1577281 3.094898e-02
## dens:env 0.0028091316 0.009296323 0.3021766 7.625174e-01
## env:densS2 0.0044931860 0.012676525 0.3544494 7.230022e-01
## env:densS3 0.0006890698 0.023931772 0.0287931 9.770296e-01
```

Notice that the coefficient names are different from those in the
`example_param`

tables. In order to import the GLM table into
the structure, we also need to specify the name equivalences:

```
glm.coef.equivalence <- list("alpha" = "(Intercept)",
"beta1" = "env",
"beta2" = "dens",
"beta3" = "densS2",
"beta4" = "densS3",
"beta5" = "dens:env",
"beta6" = "env:densS2",
"beta7" = "env:densS3")
```

And now, we can import these values. Assuming that these coefficients
refer to the survival rate of juveniles for species `s1`

at
site `sa`

:

```
example_param <- generate_vital_rate_coefs(example_param,
sp = "s1",
sites = "sa",
vital.rate = "Sj",
glm.object = glm_example_coefs,
glm.coef.equivalence = glm.coef.equivalence)
# the resulting table. Note that we only entered the coefficients for survival of juveniles (Sj)
example_param[["s1"]][["sa"]]
```

```
## alpha beta1 beta2 beta3 beta4 beta5 beta6
## Sj 1.038831 1.345738 0.06997734 -0.07880205 -0.05998993 0.002809132 0.004493186
## Sn NA NA NA NA NA NA NA
## Sr NA NA NA NA NA NA NA
## Rn NA NA NA NA NA NA NA
## Rr NA NA NA NA NA NA NA
## D NA NA NA NA NA NA NA
## Ds NA NA NA NA NA NA NA
## O NA NA NA NA NA NA NA
## beta7
## Sj 0.0006890698
## Sn NA
## Sr NA
## Rn NA
## Rr NA
## D NA
## Ds NA
## O NA
```

Otherwise, coefficients can be entered manually:

**Metapopulaiton example**

In `metapopulation_example_param`

, the coefficients for
each demographic rate were obtained from GLMs fitted on simulated data
from a relatively simple community with interaction across trophic
levels: one prey species (S1), one meso-predator (S2), and one top
predator (S3). You can think of a system such as rabbit, fox, and
lynx.

The underlying data were simulated using the following rules:

- The prey (S1) is relatively more sensitive to environmental fluctuations than either the predators. Survival of prey (across all stages) is negatively affected by the predators, and more so by the top predator (S3), which is assumed to be more specialized on the prey.

The meso-predator (S2) is a generalist, and its survival across all stages is positively affected by prey density, but not as strongly as the survival of the top predator; it competes with the top predator and is more strongly negatively affected by its density than the other way around.