# Metapopulation projections

#### David Garcia-Callejas, Maria Paniw, and cxr team

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.

Conceptual overview of metapopulation functions. We model stage-specific demographic transitions from t to t+1 (here, t = one year) using a matrix population model for populations of n species in h sites. Stages are juveniles, J, non-reproductive, N, and reproductive, R. Stage transitions depend on demographic rates of survival (S), transition probability to the reproductive stage (T), and recruitment (O). Non-reproductive individuals can also disperse (DN) and have a probabaility of SDN to survive dispersal. The example logistic regression for species 1 shows that rates are modeled as functions of an environmental variable (E) and intra- (here NS1) and interspecific (here NS2, NS3, NS4) density.The different colors highlight that the mean demographic rates (the α in the models) differ among species and habitats.

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:

    initial.densities[[1]]

    ##      [,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.

    example_param[["s1"]][["sa"]]

    ##    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:

    data("glm_example_coefs",package = "cxr")

    glm_example_coefs

    ##                  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:

    example_param$s2$sa["Sr",5] <- 0.001

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)

1. 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)

2. 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.

   The top-predator (S2)

3. The top predator (S3) is specialized on the prey and therefore responds strongly (and positively) to prey abundance. It is negatively affected by the meso-predator, but with caveats as explained above.

4. All three species disperse between 2 sites. Dispersal is highest for the top predator, followed by the meso-predator, and prey. Dispersal mortality takes the opposite pattern, being highest for the prey.

Generating coefficient randomly

And, finally, coefficients can also be generated by sampling from a normal distribution with a given mean and standard deviation. Below we enter the coefficient for the intercept (`vr.coef = "alpha"`) of the survival of reproductive adults (`vital.rate = c("Sr")`), for the first species (`sp = "s1"`) and all sites (by not specifying the `sites` argument), setting the mean and sd values from which random samples will be drawn.

    example_param <- generate_vital_rate_coefs(param = example_param,
                                       sp = "s1",
                                       # sites = "sa",
                                       vital.rate = c("Sj"),#,"Sn","Sr","Rn","Rr","D","Ds","O"),
                                       vr.coef = "alpha",
                                       mean.coef = .1,sd.coef = 0)