CATS: A high‐performance software framework for simulating plant migration in changing environments

Considering local population dynamics and dispersal is crucial to project species' range adaptations in changing environments. Dynamic models including these processes are highly computer intensive, with consequent restrictions on spatial extent and/or resolution. We present CATS, an open‐source, extensible modelling framework for simulating spatially and temporarily explicit population dynamics of plants. It can be used in conjunction with species distribution models, or via direct parametrisation of vital rates and allows for fine‐grained control over the demographic and dispersal processes' models. The performance and flexibility of CATS is exemplified (i) by modelling the range shift of four plant species under three future climate scenarios across Europe at a spatial resolution of 100 m., and (ii) by exploring consequences of demographic compensation for range expansion on artificial landscapes. The presented software attempts to leverage the availability of computational resources and lower the barrier of entry for large‐extent, fine‐resolution simulations of plant range shifts in changing environments.


| Demographic model
The demographic model is implemented as a mechanistic formulation of a stage-structured matrix model distinguishing seeds, juveniles and adults (see Figures S1 and S2). Juveniles are defined as nonreproductive plants from the seedling stage up to the year before the first reproduction, structured into annual age classes. Adults are reproductive individuals. Each ramet of a clonally propagating species is regarded as an individual (Figure 1).
The processes driving local population dynamics are quantified by environment-dependent and (facultatively) density-dependent vital rates (Table S1), which can be summarised into germination, maturation, seed production, clonal growth and year-to-year survival.
Overall, population growth is limited by the carrying capacity with surplus individuals being discarded. All demographic processes facultatively contain demographic stochasticity (drawing the number of individuals transitioning among stages from a Poisson distribution) resulting in different spatial distribution patterns at replicated runs.
CATS supports two modes to specify the spatiotemporal variation of vital rates across the study area in function of the variation in environmental conditions. In hybrid parametrisation mode, occurrence probabilities projected by SDMs serve as a proxy for a cell's suitability for a species and are translated by link functions into vital rate values which range between zero and a user-specified maximum.
The link function between suitability and each vital rate (as well as the carrying capacity) is sigmoid ( Figure 2) by default (but can be customised as linear or constant): where x is the suitability, f max is the vital rate's maximum value, OT is the occurrence threshold, a is a slope parameter, defaulting to 15, and s is the scale factor (Supporting Information 1.3.1). At a suitability value of x = OT, the vital rate value will be equal to f max × s. The user can specify each vital rate as density dependent requiring an adaptation of the vital rate function: where K is the maximum carrying capacity and N is the number of adults in the cell. This ensures λ ≥ 1 for suitability values x > OT and λ = 1 for x = OT, independent of abundance. Link functions allow for further modifications (Supporting Information 1.3).
These link functions are calibrated to guarantee stable population size (i.e. a constant number of adults in a cell) at a user-defined value of suitability. The rationale is that probabilistic SDM projections are usually translated to presences/absences using an occurrence probability threshold (e.g. Liu et al., 2005). This threshold conceptually represents the conditions of marginal suitability for a species (i.e. population growth rate λ = 1, Chase & Leibold, 2003). Sensible choices of such thresholds, such as values that maximise some SDM evaluation statistics, are provided by most SDM software packages and can directly be used for specifying the respective parameter OT (occurrence threshold). Population stability at OT is guaranteed by selecting the value of s in Equation (1) such that the dominant eigenvalue of the Leslie-matrix formulation of population dynamics (in an isolated cell) is equal to unity (see Supporting Information 1.3.1). This value of s is common to all vital rates, and represents the fraction of their maximum realised at OT. (1)

| Dispersal model
In CATS, seeds (or other dispersal units) are distributed using one or more two-dimensional dispersal kernels, each representing one dispersal vector. Dispersal kernels are user-provided raster files comprising the probability that a seed originating from the raster's centre cell will be dispersed into that particular cell. The number of seeds dispersed from a source to a target cell is the product of the number of seeds pro-

| Performance
To decrease computing time, CATS is multi-threaded. This parallelisation is realised by dividing the study area in spatially nonoverlapping regions. Thereby the demographic model, which simulates each cell in isolation, can be run in parallel for each region. For dispersal, which connects cells within and also among regions, parallel computation is restricted to regions spatially separated more than twice the maximum dispersal distance.
CATS is also available as Message Parsing Interface version, in which a single simulation is run across multiple computers. For more details see Supporting Information 5.1.

| E X AMPLE S
We showcase the functionality of CATS with one example for each demographic parametrisation mode. In the first example, using hybrid parametrisation, we individually model four species of the genus Astragalus, which represent the host plant pool of the southern European butterfly Polyommatus escheri. Their spatiotemporal range dynamics within this century were modelled based on three climate change scenarios to evaluate whether their supposedly low mobility might hinder range adaptation of the butterfly. The study-area grid covered more than 0.5 billion terrestrial cells.
The second example uses the direct parametrisation mode to explore how negative correlations between individual vital rates (Villellas et al., 2015) affect spatiotemporal population dynamics. We therefore modelled a species in a virtual landscape representing a temperature gradient. Simulations assuming all vital rates to peak at intermediate temperatures were compared to simulations assuming either adult survival, sexual reproduction or asexual reproduction to have a low at intermediate temperatures.
More detailed information on the parametrisation of the example simulations is provided in Table S4.

| RE SULTS
Under climate warming, all four Astragalus species tend to shift to the north ( Figure 3). However, the fine resolution simultaneously allows to represent fine-grain patterns of species range shifts as they emerge in mountainous terrain in particular. Overall, range margins are predicted to shift no further than 25 km. This will likely put a strong constraint on the potentially faster range shifts of the more mobile butterfly species that depends on these plants. The simulations of 80 years (30 years burn-in) for all four species (up to 38 million populated cells) took 4 h F I G U R E 3 Aggregated results of simulations of range dynamics of four Astragalus species across Europe (grid resolution: 100 m) between 2020 and 2100 under moderate climate change (RCP45) using the hybrid parametrisation mode. The number of species projected to occur in a cell in the year 2020 is subtracted from the number of species occurring in 2100. Blue, red and dark grey indicate number of species to increase, decrease and being stable respectively. Sharp borders emerge from the coarse definition of the species' current range (spatial distribution of initial populations). The inset illustrates the elevational pattern in an area of the French Massif Central with higher elevations gaining and lower elevations losing species. Results assuming weak (RCP26) and strong (RCP85) climate change are shown in Figure S1. The second example shows that all investigated, negative correlations between vital rates slow down population spread, but to a very different degree, and with different side effects (Figure 4).
Negative correlation of asexual reproduction with the other rates slightly reduces final range size, but does not affect abundance patterns. In contrast, an inverse pattern of sexual reproduction strongly constrains range size and tends to make abundance distributions more homogeneous across the environmental gradient. An inverse pattern of adult survival, finally, has an intermediate effect on range expansion, but the most extreme effect on abundance distribution, leading to consistently low abundances across the entire gradient.
These results confirm the high importance of adult survival for the local demography of perennial herbs, in particular their local abundance (e.g. Weppler et al., 2006). In addition, we show a considerable effect of the rate of sexual reproduction on the spatial spread. CATS requires the user to provide a priori fitted dispersal kernels in a specified format. Although this comes at the cost of extra effort necessary to produce the kernel files, it offers a high degree F I G U R E 4 Results of simulations in direct parametrisation mode used to exemplify the effects of demographic compensation. The spread of a virtual species was simulated across a landscape (500 × 2500 cells) representing an environmental gradient. Vital rate values were constant along the vertical axis, with noise drawn from a normal distribution, but varied on the horizontal axis as illustrated in the top-left panel, that is, either with a peak (Layer A) or a low (Layer B) at intermediate gradient values. Coloured panels represent simulation results after 50, 100 and 200 years under four scenarios assuming different subsets of vital rates to vary according to Layer B as indicated by text labels (all other vital rates followed Layer A). Colours represent the abundance of adults. Panels in the rightmost column show population growth rate's horizontal profile as a property of the environment (i.e. along the horizontal axis), that is, without density effects or seed influx.

| DISCUSS ION
Vital rate values were calculated by weighing their maximum according to the values in Layer A or B. All simulations started with the same initial distribution, with adult abundance at carrying capacity at the bottom row of cells and zero elsewhere.