CC-BY 4.0
Introduction
Climate is changing at a rate rarely seen in the geological record (Burke et al., 2018; Loarie et al., 2009). Paleoecological and macroecological studies link periods of climatic instability to elevated extinction rates and to the formation of many patterns of endemism (Harrison & Noss, 2017; Song et al., 2021). Extinctions, however, are typically multicausal. Still, the rapid magnitude of contemporary and projected climate and land-use changes is expected to increase species’ vulnerability and—especially where climate change interacts with land-use change, exploitation, or biological invasions—substantially raise extinction risk (Barnosky et al., 2011; Cowie et al., 2022).
When faced with rapid environmental change, species can (i) resist new conditions, (ii) adapt, (iii) shift their distributions, or (iv) go extinct (Berg et al., 2010; Corlett & Westcott, 2013; Duputié et al., 2015). Although some adaptation and resistance are likely, range shifts are a widespread response to past and ongoing climate change, observed across many taxa including butterflies and plants (Habel et al., 2021; Graae et al., 2018; Lenoir et al., 2020). Many species are shifting poleward and upslope in response to warming (Lenoir & Svenning, 2015), but shifts may also follow other climatic drivers (e.g., water balance), anthropogenic influences, or complex multifunctional gradients, producing varied and sometimes unexpected directions of movement. (Dobrowski et al., 2013; Serra-Diaz et al., 2016b; VanDerWal et al., 2013; Wason & Dovčiak, 2017).
Empirical and modeling studies indicate that the pace of climate change often exceeds the rate at which species can naturally migrate, a problem especially acute for sessile terrestrial organisms. Slow migration produces a disequilibrium between species’ distributions and climate: trailing-edge contract but climatically suitable areas at the leading edge are unoccupied, thus increasing extinction risk (Svenning & Sandel, 2013; Corlett & Westcott, 2013; Lenoir et al., 2020).
Assisted migration (AM)–also termed assisted colonization, managed relocation, assisted range expansion and species’ translocation (Hällfors et al., 2014)– refers to the introduction of a species to an area where it is not present to support species natural migration as a consequence of climate change (Hayward, 2009; Heller & Zavaleta, 2009; Hewitt et al., 2017; Hoegh-Guldberg et al., 2008; Richardson et al., 2009). AM could be an effective strategy to compensate for the slow movement of some species in order to prevent their decline or extinction (Hällfors et al., 2018; Peterson St-Laurent et al., 2018). Nevertheless, AM is still a contentious conservation strategy (Aubin et al., 2011; Hunter, 2007; McLachlan et al., 2007; Ricciardi & Simberloff, 2021; Vitt et al., 2010), as it could lead to a number of undesired outcomes like disease transmission (Simler et al., 2019), unintended establishment of invasive species (Mueller & Hellmann, 2008) or lack of local adaptation (Tíscar et al., 2018).
Most current tools for forecasting range shifts are correlative species distribution models (SDMs) that estimate relationships between occurrences and environmental predictors and then map habitat suitability. SDMs typically project distributions for discrete climate normal periods (e.g., baseline, 2050, 2070) and thus produce static snapshots that summarize 20–30-year climatic averages (Franklin, 2010). These snapshots can identify potentially suitable future areas, but they do not capture transient processes — such as the timing of dispersal events, establishment lags, sequential colonization, or the effects of stochastic disturbance — that determine if and when those areas will become occupied.
Alternatively, there exist mechanistic or dynamic distribution modeling which explicitly account for key migration processes such as dispersal (Lehsten et al., 2019), evolution (Bocedi et al., 2021; Bush et al., 2016), disturbances (Liang et al., 2018; Serra-Diaz et al., 2015), and species interactions (Keyel et al., 2016). Some mechanistic models even incorporate detailed temporal dimension (e.g. Moving-Habitat Models (Harsch et al., 2017)). However, those models involve a high level of parametrization and may require computationally intensive simulations of range shifts. This complexity prevents us from exploring all possible combinations of introductions in terms of time and space, and therefore from optimizing an AM planification.
To date, quantifications for the need of assisted migrations have used metrics derived from species distribution models (Hällfors et al., 2016, 2017), but to our knowledge no spatial conservation planning method has dealt with the issue of assisted migration with efficiently allocate species assisted migration efforts in space and time. (Chu & Beasley, 1998)
Here, we present a new method implemented through an R package that explores the effects of introductions in space and time on species range shifts, and further optimizes the allocation of introductions. It is built to achieve a defined conservation objective — maximizing expected range size and minimizing extinction risk under spatially explicit constraints. We call this novel mathematical framework SIESTR (Introduction Effects in Space and Time on species Range dynamics). The algorithm seeks to meet a probabilistic minimum range area target, with a predefined level of confidence (e.g. 90%) by a given time horizon (e.g. 2100). This approach is computationally efficient because it relies on precomputed transition matrices to summarize dispersal and species environmental suitability dynamics, allowing rapid evaluation of introduction scenarios through metaheuristic optimization rather than full mechanistic simulations. We then show a case study where SIESTR enables the determination of an assisted migration plan for a virtual species in decline due to climate change.
The method has been implemented in a newly developed R package, designed specifically for this study, which is available on GitHub (https://github.com/ArndCallebaut/IESTR/).
Methods
The method consists of two main phases: preprocessing and optimization (Figure 1). During the preprocessing phase, transition matrices are constructed. These transition matrices are then used in the optimization phase to identify optimal species introductions in space and time.
The study area is represented as a spatial grid. Cell size may vary depending on data resolution and species’ ecology but should be chosen to reflect the establishment and dispersal dynamics of the focal species. Time is discretized into regular time steps. Two key processes are modeled: (i) species dispersal among grid cells and (ii) species survival as a function of environmental conditions within each cell at a given time.
The algorithm requires at least three input datasets (initial presences on the grid, suitability matrices over time and migration kernel) and a set of optimization parameters. Additional parameters can be specified to control the behavior of the optimization phase.
Input data and parameters
SIESTR requires three types of input data and an optional input (Figure 1, left box). The first user-input required data is the migration kernel – e.g. matrix of dispersal probability. For a given species and a given grid map, the migration kernel corresponds to the probability that the species spread into the neighboring cells during one step time. These dispersal probabilities should reflect the spatial resolution of the grid (e.g. distances) and the time step (e.g. years, decades), and it is assumed to be constant over time. This migration kernel should be a function of cell size and the dispersal capacity of the species. Knowing the dispersal range of the species, the migration kernel could be set as an isotropic exponential decay function of distance. The second required user-input is the set of suitability matrices over time. For each discrete time t, we define a suitability matrix as a grid with values between 0 and 1 corresponding to the survival chance of the species in each cell of the map between time t and time t + 1 after the dispersal process. Suitability indices may be derived from climatic conditions and/or other factors such as land use, soil properties or species interactions., and can be, for example estimated from correlative SDMs and rescaled to match demographic probabilities between t and t + 1. The third input is the occurrence map of the species at the initial time. This is a grid with a value of 0 for the absence of the species and 1 for the presence.
The parameterization of these different model inputs depends on the time and space scale used. Two examples at different scales are provided to establish each parameter in Annex S3.
The algorithm requires the user to specify the optimization objective for species introductions (Figure 1, left panel). SIESTR is currently designed to address two types of assisted migration (AM) problems, both formulated as constrained optimization tasks: (i) minimizing introduction costs while achieving a fixed target range size, or (ii) maximizing range size under a limited budget. Mathematically, both formulations resemble a knapsack-type problem, in which a set of introduction actions across space and time is selected to optimize an objective function (range size or cost) while satisfying explicit constraints (budget limits or minimum occupied area).
Several optional inputs can be provided to refine the optimization phase. The first is an introduction cost matrix. By default, introduction costs are assumed to be uniform across all cells, but users may specify spatially explicit costs to reflect economic expenses, management or land-tenure constraints (e.g. protected vs. private land), or ecological considerations such as risks to recipient ecosystems. Additional optional inputs include parameters controlling the genetic algorithm that is used to find the best set of introductions, which influence convergence speed and solution precision. These parameters are automatically derived from map size by default but can be manually adjusted if needed (see code documentation).
Figure 1 - Scheme of the SIESTR workflow. The workflow is divided into two main phases from inputs two outputs: (1) pre-processing of inputs (light blue box), and (2) optimization for introduction choices in space and time (dark blue box). Arrows detail inputs and outputs for each process.
Preprocessing
Step 1: Dispersal process - Spread matrix
The first part of the preprocessing builds the spread matrix and is based on neutral processes. This matrix defines species accessibility in space over one-time step. It is calculated using the migration kernel and it contains the probability of colonization site to site only considering the dispersal — not the survival — process. If the user has precise information about species spread rates, they can directly calculate them instead of using a migration kernel. This can be useful if the dispersal potential isn’t spatially uniform.
The migration kernel is an input of the model, in grid form, and it represents the probability of a cell colonizing the neighboring cells in a time step. The size of the dispersal matrix is at least 3 x 3 cells, dispersal is assumed isotropic and dispersal is evaluated at each time step of the simulation.
For ease of notation, we introduce the special addition operator, \(\widetilde{+}\), such as for two real numbers \(p_{1}\) and \(p_{2}\) with values between 0 and 1 (Eq1):
\(\begin{array}{r} p_{1}\widetilde{+}p_{2} = p_{1} + p_{2} - p_{1}p_{2}\# \end{array}\)
Specifically, by adding the probabilities of two independent Bernoulli events with this operator, we obtain the probability of success of at least one of them. This will enable us to add up the probabilities of colonization of the same site from different sources.
The probability of colonization of site A by site B is the probability of finding the species on site B multiplied by the probability that presence on site B at time t will lead to presence on site A at time t+1. With the special addition, we can add up the probabilities of colonization of the same site from different sources. Let t be a time-step, i and j two sites of the map, n the number of sites, \(M_{i,j}\) the probability that a presence in the site j leads to a presence in the site i at the next timestep, and \(X_{j}^{t}\) the probability of presence in the site j and at the time t of the species. Then the probability \(C_{i}^{t + 1}\) that the site i is colonized by the species at the time t+1 is, with our notation (Eq2):
\(\begin{array}{r} C_{i}^{t + 1} = \ M_{i,1}X_{1}^{t}\ \widetilde{+}M_{i,2}X_{2}^{t}\widetilde{+}\ldots\widetilde{+}M_{i,n}X_{n}^{t} = \sum_{j \in \lbrack 1,n\rbrack}^{\sim}{M_{i,j}X_{j}^{t}}\# \end{array}\)
With this notation, the probability that site B will be colonized is the ‘special sum’ operator of the probabilities of being colonized by each of the other sites. This formula can be seen as a scalar product without linear properties. Let X and Y be two vectors of size N with values between 0 and 1. We define the operator between two vectors \(\left\langle . \middle| . \right\rangle\) as (Eq3), with the special sum described in (Eq1) :
\(\begin{array}{r} \left\langle X \middle| Y \right\rangle = \ \sum_{i \in \lbrack 1,\ \ n\rbrack}^{\sim}X_{i}Y_{i}\ \# \end{array}\)
As we defined an operator between two vectors of the same size, we can associate a matrix product with it. Let A and B be two matrices of a size m*n and n*k respectively, with values between 0 and 1. Then we define the operator between two matrices \(\widetilde{*}\) as (Eq4):
\(\begin{array}{r} \ A\widetilde{*}B\ = \ \begin{pmatrix} \left\langle A_{1,\ \bullet} \middle| B_{\bullet ,1} \right\rangle & \cdots & \left\langle A_{1,\ \bullet} \middle| B_{\bullet ,k} \right\rangle \\ \vdots & \ddots & \vdots \\ \left\langle A_{m, \bullet} \middle| B_{\bullet ,1} \right\rangle & \cdots & \left\langle A_{m,\ \bullet} \middle| B_{\bullet ,k} \right\rangle \end{pmatrix}\ \# \end{array}\)
With this notation, we can express the vector \(C^{t + 1}\) with the spread matrix M and the vector of presence at the time t (Eq5):
\(\begin{array}{r} C^{t + 1} = \begin{pmatrix} C_{1}^{t + 1} \\ \vdots \\ C_{n}^{t + 1} \end{pmatrix} = \ \begin{pmatrix} \left\langle M_{1,\ \bullet} \middle| X^{t} \right\rangle \\ \vdots \\ \left\langle M_{n, \bullet} \middle| X^{t} \right\rangle \end{pmatrix}\ = M\widetilde{*}X^{t}\ \# \end{array}\)
Step 2: Survival process - Building the transition matrices
The transition matrices merge the information of the spread matrix (step 1) and the species’ suitability of each site over time (Figure 1). Transition matrices thus allow us to calculate the probability of presence in a future time step, while considering the initial distribution of the species. Let \(s_{i}^{t}\) be the suitability value of the species in the site i, at time t. The probability of species occurrence in the site i at the time t +1, \(X_{i}^{t + 1}\), is calculated as the probability that the site is colonized at t+1 (\(C_{i}^{t + 1}\)) multiplied by the probability of the species establishment and survival (suitability) \(s_{i}^{t}\) in this cell (Eq6):
\(\begin{array}{r} X_{i}^{t + 1} = \ \ s_{i}^{t}C_{i}^{t + 1}\ \# \end{array}\)
The probability of a species being present at time t+1 in a site will be equal to the probability of colonization of this site, multiplied by the probability of survival of the species over the time step [t, t+1] in this this site.
By doing this for the whole vector \(X^{t + 1}\), we obtain a formulation for the transition matrix \(T_{t \rightarrow t + 1}\)that allows us to calculate presence probabilities from \(X^{t}\) to \(X^{t + 1}\) (Eq7):
\(\begin{array}{r} X^{t + 1} = \begin{pmatrix} X_{1}^{t + 1} \\ \vdots \\ X_{n}^{t + 1} \end{pmatrix} = \ \begin{pmatrix} s_{1}^{t}\left\langle M_{1,\ \cdot} \middle| X^{t} \right\rangle \\ \vdots \\ s_{n}^{t}\left\langle M_{n, \cdot} \middle| X^{t} \right\rangle \end{pmatrix} = \ \begin{pmatrix} \left\langle s_{1}^{t}M_{1,\ \cdot} \middle| X^{t} \right\rangle \\ \vdots \\ \left\langle {s_{n}^{t}M}_{n, \cdot} \middle| X^{t} \right\rangle \end{pmatrix}\ = \ T_{t \rightarrow t + 1}\widetilde{*}X^{t}\ \\ with\ :T_{t \rightarrow t + 1} = \ \begin{pmatrix} s_{1}^{t}M_{1,\ 1} & \cdots & s_{1}^{t}M_{1,\ n} \\ \vdots & \ddots & \vdots \\ s_{n}^{t}M_{n,\ 1} & \cdots & s_{n}^{t}M_{n,n} \end{pmatrix}\# \end{array}\)
The formula is recursive. Let t1 and t2 be two times such as t1 < t2, then we can obtain a formulation of \(X^{t_{2}}\ with\ X^{t_{1}}\), by combining the different transition matrices we already have (Eq8) (Eq9):
\(\begin{array}{r} X^{t_{2}} = T_{t_{2} - 1 \rightarrow t_{2}}\widetilde{*}X^{t_{2} - 1} = \ T_{t_{2} - 1 \rightarrow t_{2}}\widetilde{*}\left( T_{t_{2} - 2 \rightarrow t_{2} - 1}\widetilde{*}X^{t_{2} - 2} \right) = \ T_{t_{2} - 1 \rightarrow t_{2}}\left( \ldots\ \widetilde{*}\left( T_{t_{1} \rightarrow t_{1} + 1}\widetilde{*}X^{t_{1}} \right) \right)\ \# \end{array}\)
\(\begin{array}{r} X^{t_{2}} = \left( T_{t_{2} - 1 \rightarrow t_{2}}\widetilde{*}(T_{t_{2} - 1 \rightarrow t_{2}}\widetilde{*}(\ldots\widetilde{*}T_{X^{t_{1}} \rightarrow X^{t_{1} + 1}} \right)\ldots))\widetilde{*}X^{t_{1}} = \ T_{X^{t_{1}} \rightarrow X^{t_{2}}}\widetilde{*}X^{t_{1}} \\ with\ :\ T_{t_{1} \rightarrow t_{2}} = \left( T_{X^{t_{2} - 1} \rightarrow X^{t_{2}}}\widetilde{*}(T_{X^{t_{2} - 2} \rightarrow X^{t_{2} - 1}}\widetilde{*}(\ldots\widetilde{*}T_{X^{t_{1}} \rightarrow X^{t_{1} + 1}} \right)\ldots))\# \end{array}\)
The calculation of \(T_{t_{1} \rightarrow t_{2}}\) in (Eq9) must be made in the right order for the equation to be true, as the \(\widetilde{*}\) operator is not associative. With this equation, we can calculate the probability of presence of the species at any point in time, knowing the probability of presence at a particular time in the past and these transition matrices.
To be able to evaluate and compare species introduction(s) effect, we set the last time-step of study, \(t_{\max},\) as a target time-step, and focus on all the transition matrix of the form \(T_{\cdot \ \rightarrow t_{\max}}\). These matrices allow for calculating the effect of introduction(s) on the probability of the presence of a species in any cell, at any time. An introduction in a cell \(i\) and a time \(t\) leads to a probability of presence in each of the cells at the end equal to the vector \({(T_{t\ \rightarrow t_{\max}})}_{\cdot ,i}\).
The effect of the multiple species presence can be summed up, as the probability of presence in a cell in the final system state is the probability that at least one of the present species (initially present or introduced) establishes in this cell. Let \(t_{0}\) be the initial time, K the number of introductions, \(\lbrack(i_{1},t_{1})\ldots(i_{K},t_{K})\rbrack\) The different choices of introduction defined as pairs of a cell (i) and a time-steps (t). Then, the probability of presence in all the cells in the final state of the system is (Eq10):
\(\begin{array}{r} X^{t_{\max}} = \ T_{{t_{0}}^{\ } \rightarrow t_{\max}}\widetilde{*}X^{t_{0}}\widetilde{+}\sum_{k \in \lbrack 1,\ K\rbrack}^{\sim}\left( T_{t_{k} \rightarrow t_{\max}} \right)_{\cdot ,\ i_{k}}\ \# \end{array}\)
The construction of these transition matrices makes it possible to obtain the results of a simulation of the introduction of species over time by a direct calculation. A quantitative worked example of its use can be found in Appendix S1.
Optimization
We calculate the effect of introducing the species in different sites and times based on transition matrices and implement optimization algorithms to evaluate alternative introduction schemes.
Optimization problem formulation
The aim is to find the best set of introductions (site, time) to satisfy the constraint while optimizing the objective function. Two kinds of problems are currently implemented in SIESTR: (1) the target range area (TRA) and (2) the best allocation of a budget (BAB). In TRA we optimize the minimum cost needed to maintain the species at a targeted range size, whereas in BAB we maximize species range area given a maximum budget.
To implement these optimizations, we need to define a choice evaluation function (CEF), that allows us to summarize the efficiency of a set of introductions. This function takes as parameters a set of introductions choices, the probability of presence in each of the sites at the end of the period without introductions (e.g. suitability and dispersal), the transition matrices, and a precision parameter set by the user α (typically 90%). The function returns the minimal number of sites colonized at the end of the period considering introduction choices, with an α%probability. It means that given the probabilities of presence at the end of the whole time period in every cell (obtained with the transition matrices), in α% of the cases, the number of cells with a presence will be equal or higher than the function output number.
In the context of the TRA optimization problem, the value of this function should be equal or higher that a set minimal range area. In the BAB, the value of this function should be maximized while keeping the total cost under a set threshold.
The two problem formulations are currently expressed in Eq. 11 for the TRA problem, and Eq. 12 for the BAB problem, with VC being the set of viable choices of introductions to choose from:
\(\begin{array}{r} \left\{ \begin{array}{r} \min_{A \subset VC}{\sum_{(i,t) \in A}^{}{cost(i)}} \\ CEF\left( T_{t_{0} \rightarrow t_{\max}}\widetilde{*}X^{t_{0}}\widetilde{+}\sum_{(i,t) \in A}^{}\left( T_{t_{\begin{array}{r} \ \\ \ \end{array}} \rightarrow t_{\max}} \right)_{\cdot ,\ i},\alpha \right) \geq minimal\ final\ range\ area \end{array} \right.\ \# \end{array}\)
\(\begin{array}{r} \left\{ \begin{array}{r} \max_{A \subset VC}{CEF\left( T_{t_{0} \rightarrow t_{\max}}\widetilde{*}X^{t_{0}}\widetilde{+}\sum_{(i,t) \in A}^{}\left( T_{t_{\ } \rightarrow t_{\max}} \right)_{\cdot ,\ i},\alpha \right)} \\ \sum_{(i,t) \in A}^{}{cost(i)} \leq maximum\ total\ cost \end{array} \right.\ \# \end{array}\)
In case the costs are not specified, the sum of costs becomes equivalent to the number of introductions.
Step 3: Filtering of introductions and occurrence probability predictions without introductions.
We reduce all potential combinations of introductions in space and time to a smaller set by using two filters. The result of this pre-evaluation will lead us to assess the optimization only using this subset of viable solutions among the panel of all introduction sets.
The first filter assesses the effectiveness of the introduction at a site at a given time. Each cell and time of introduction (i,t) is assigned in the transition matrices a probability vector of presence at the end of the study period. Summing the values of this vector gives the mean number of sites where the species will be present at the end of the period by introducing the species at that site at the corresponding time (Eq13):
\(\begin{array}{r} eval\left( \left( T_{t_{\ } \rightarrow t_{\max}} \right)_{\cdot ,\ i} \right) = \ \sum_{j \in \lbrack 1,n\rbrack}^{}\left( T_{t_{\ } \rightarrow t_{\max}} \right)_{j,\ i}\ \# \end{array}\)
A threshold parameter controls, for each potential choice of introduction, if this choice offers a mean number of colonized sites by the end of the study period sufficient. We have set a default value of 1 for this threshold, which filters introductions that would bring on average less than one presence of the species at the end of the whole time period. This threshold can nevertheless be modified. Many possibilities for a small number of introductions could lead to a more drastic elimination of the inefficient sites, and therefore a need to increase the threshold.
The second filter discards sub-optimal times of introduction. For a set of times of introduction for a given site, the filter eliminates all the times that are giving lower probabilities of presence at the end of the period for all sites, compared to the other sites. In other words, if it’s strictly better to introduce a species at a time compared to another time, the second time will be discarded from potential solutions, as there would be no reason to choose it over the other.
Genetic algorithm
We implemented a genetic algorithm to optimize the introductions in space and time for a given goal (see input parameters) on a set of viable solutions (see Pre-evaluation section). A genetic algorithm is an evolutionary algorithm whose goal is to find the most optimal solution to an optimization problem, potentially under constraints, with a mechanism inspired by biological evolutionary phenomenon (Chu & Beasley, 1998). The objective here for the algorithm is to find the set of pairs of cells and time solutions to the problems (2.3.1). Details on the algorithm workflow and parametrization are specified in Appendix S2.
A case study on assisted migration
SIESTR is designed as a general framework that can be adapted to different spatial extents and resolutions. Input data and spatial scale depend on the ecological characteristics of the focal species and must therefore be defined on a case-by-case basis. We provide guidance on appropriate data products and spatial scales, with a particular focus on tree species, for which assisted migration programs are already underway (Appendix S3). We used a virtual case study rather than a direct empirical application. Unlike real-world case studies, whose results are inherently specific to the chosen species and landscape, a virtual system allows us to isolate and better understand the behavior of the model under different constraints and assumptions.
We applied SIESTR to a simulated assisted migration scenario for a plant species. The study area is a circular island characterized by a climatic gradient generated by a diagonally oriented mountain range (Figure 2a). The island is represented as a 50 × 50 grid, with each cell representing a potential site of species presence or absence. Climate change is simulated as an average temperature increase of 1.5 °C across the island, with warming occurring more rapidly at higher elevations (Figure 2b,c). The target species is initially present in 50 cells located at the margin of its ecological niche (Figure 2d), forming two disjunct populations (Figure 2e). The objective is to identify optimal assisted migration locations and timing that prevent range contraction under climate change between 2030 and 2060, while maintaining the species’ occupied area at 50 cells. To simplify the system, habitat suitability is based solely on temperature. Temperature in each cell is determined by three factors: altitude (cooler at higher elevations), latitude (warmer in the south than in the north), and time (progressive warming). The climate change simulation spans 30 annual time steps, corresponding to the 30-year study period.
The equation used is the following, with altitude, latitude and time factors between 0 and 1, and with \(T_{\min,alt = 0}\), \({\mathrm{\Delta}T}_{latitude}\), \({\mathrm{\Delta}T}_{altitude}\), \({\mathrm{\Delta}T}_{climatechange}\) as external factors:
\(\begin{array}{r} Temperature = T_{\min,alt = 0} + \ {\mathrm{\Delta}T}_{latitude}*latitude + \\ {\mathrm{\Delta}T}_{altitude}*altitude + \ {\mathrm{\Delta}T}_{climatechange}*time\ \# \end{array}\)
The species input parameters to the algorithm consist of a migration core matrix, a set of matrices describing species probability of surviving from t to t+1.over space and time, and a matrix identifying the current species distribution. We set the migration kernel as a 3x3 matrix with a value of 1 in the center representing the current presence of the species and 0.01 in all other cells. That is, a 0.01 probability of colonizing each year each of the eight cells around the center cell. The suitability matrices are based on species climatic suitability on the island through time (50 timesteps). The suitability is calculated from the temperature as follows, with the ordered external variables \(T_{min,opt}\) \(T_{max,opt}\ \), \(T_{\max}\)and \(T_{min,opt}\), where\(T_{\max}\)and \(T_{\min}\ \)represent species niche margins while \(T_{min,opt}\) and \(T_{max,opt}\) limit the optimal part of the range where the survival rate is the highest (Figure 2d).
We applied SIESTR to three assisted migration situations: Experiment 1 – range dynamics without assisted migration operations (Experiment 1, Figure 3); Experiment 2 – with assisted migration using equal costs among all land (Figure 4), and Experiment 3 – using introduction costs decreasing with the altitude (Experiment 3, Figure 5). Because the optimization stems from a pseudo-stochastic process, we evaluated part of the uncertainty of the algorithm by running 100 repetitions of the algorithm and comparing the stability of the results for Experiment 2 and 3. Experiment 1 does not optimize any introduction, thus repetitions are not necessary.
Figure 2 - Topographic and climatic conditions of the simulation, and species suitability characteristics. The surface of the island is represented by a 50-column by 50-row grid, each cell being a potential site of presence/absence of the species. (a) Virtual island with a Southwest-Northeast mountain range, (b) initial climate conditions with a decreasing temperature gradient both with altitude and toward NW direction and (c) final temperature conditions at the end of the simulation period (30 time steps) showing a uniform temperature increase of 1.5 °C (d) Species niche parameters for the study species for temperature and (e) initial presence of the species at 50 sites (cells), and suitable area dynamics under simulated climatic change.
Results
The SIESTR method was implemented using R version 4.3.1, using the seed 123 (for R and C++ part). The algorithm ran during 0’12’’per run, with a use of 512 Mb memory on RAM on an Intel® Core™ i7-10875H CPU @ 2.30GHz for each of the two cost scenarios. Experiment 1, where no optimization is performed, was performed in 0’05’’. No parallelization was used.
In addition to the work presented in this section, Appendix S5 presents a sensitivity analysis of the parameters on the number of introduction sites (Figure S5.1), and a brief presentation of the calculation time as a function of the size of the map of the study problem (Table S5.1).
Experiment 1: Presence predicted without introduction(s)
Species disappeared where the climate becomes unfavorable – the southern section of the island – and colonize its new range, according to the shifts of the climate suitability through time (Figure 3). At the end of the period, with a 90% certainty, the species is at least present in 18 sites: 10 initial sites (Figure 3, blue cells) and at least 8 colonized sites (Figure 3 in green cells), while it was present in 50 sites at the start. This first output of the method thus shows species range change without introductions, and it is calculated with the transition matrix \(T_{t_{0} \rightarrow t_{\max}}\).
Figure 3 - Species distribution shifts obtained through SIESTR algorithm after 30 timesteps of climate change. The output shows occurrences that maintained (blue, probability of presence estimated higher than 90%) or lost (red, probability of presence estimated lower than 10%) as a result of shifts in climate suitability. The probability of species colonization (green gradient) along the 30 years of climate change accounts for both dispersal and climate (via species suitability) at each site. Species occurrence shifts are overlaid on climate suitability shifts in gray.
Experiment 2: Introductions with uniform costs to maintain range area
At the end of 30 years (time steps) of simulated climate change, the species survived in 10 sites (Figure 4a, blue points), and it was necessary to introduce it in 7 sites (yellow to brown points) to maintain species range area to 50 sites. Note that the algorithm does not directly predict where colonization occurs. It identifies at the end of the period what are the probabilities of presence in each cell.
Figure 4 - Introduction scheme according to SIESTR output in order to maintain the initial range area of 50 cells: (a) Introduction sites in space and time and their effect on the probability of presence of the species at the end of the period (30-time steps), and (b) temporal distribution of introduction events, with N the number of introductions at each time-step. The final range area comprises 10 initial sites, 7 introduced sites and at least 33 sites colonized via species dispersal (e.g. colonization probability >0.90).
For this experiment, introductions are found in the north and south of the island, in areas with more favorable suitability gradients (Figure 2b-c).
The results of the model suggested that planting as early as possible (that is to say, in our experience, as close as possible to year 2030) is an optimal solution to enhance species migration (Figure 4b). Introductions in areas that become unsuitable were never selected, nor introductions in areas that become suitable during the studied period.
Experiment 3: Introduction with spatially heterogeneous costs to maintain range area.
If cost surfaces are provided to the algorithm (Figure 5a), the introduction choices achieved by SIESTR differed from those obtained with uniform cost (Experiment 2). At the end of 30 years (time-steps) of simulated climate change, the species persisted in 10 sites (Figure 5b, blue points). Thus, introductions in 7 sites were necessary to maintain species range area to 50 sites (Figure 5b, yellow to brown points). This represents the same introduction number as in the case of uniform costs.
Optimal introductions are located in the Southern section of the island (Figure 5b), which is in stark contrast with the Northern introductions suggested in Experiment 2 (Figure 4b). For these sites, it should also be noted that the optimal introduction time is the same as in Experiment 2 (Figure 4c).
With a decreasing introduction cost with altitude, the algorithm seeks to minimize costs by favoring areas that have a lower cost concentrated on the north west.
Figure 5 - Introduction scheme according to SIESTR output in order to maintain the initial range area of 50 sites under non-uniform cost surfaces: (a) Cost map used for the case study, modeled as higher costs on lowlands and lower costs in the mountain tops (see Figure 2), (b) Introduction sites in space and time and their effect on the probability of presence of the species at the end of the period, and (c) temporal distribution of introduction events. The final range area comprises 10 initial sites, 7 introduced sites and at least 33 sites colonized via species dispersal (colonization probability >0.90).
Output uncertainty
Since the optimization algorithm is metaheuristic, it returns a single solution according to a stochastic algorithm. However, the algorithm may find different sets of introduction sites as solutions under multiple iterations with the same input parameters. We explored the uncertainty around AM costs using 100 iterations with the same inputs.
We observed that none of the sites were consistently selected more than 10% of the times for uniform costs (Figure 6a grey cells). For the non-uniform costs, 3 sites were selected more than 50% of the times and 9 were selected between 10% and 49% of the times (Figure 6b, grey & black cells).
We observed an overall larger area of potential introduction sites under uniform than non-uniform cost surfaces (Fig 6, light grey area), albeit with low selection rates (i.e., less than 10% of the iterations selecting the site).
Figure 6 - Spatial representation of introduction choices for 100 runs of the SIESTR algorithm. For each site, the percentage of times it was chosen as an introduction location is shown, in the (a) uniform costs scenario and (b) heterogenous costs scenario.
Discussion
In this study, we developed a new framework designed to address a central challenge in dynamic assisted migration planning: identifying when and where species introductions are most effective under climate change while accounting for uncertainty, costs, and dispersal dynamics. SIESTR, the new algorithm and R package presented here, has the potential to help decision making in spatial conservation planning. Using a virtual species case study, we demonstrated that the method is capable of meeting the target of maintaining species range size under increasing adverse climate change conditions while explicitly quantifying uncertainty in introduction decisions.
We found that under smooth climate gradients and spatial uniform costs of introduction, the algorithm suggests planting early in sites that are suitable during the whole period to facilitate species range shifts. This is because early introductions benefit most from dispersal, especially if climatic conditions are favorable for the development of the species. This is still the case when simulations introduce spatially variable costs that may counteract that effect. However, the incentive to introduce species early should be interpreted cautiously, as early introductions may also entail ecological risks, including increased exposure to future unfavorable conditions or biotic interactions not captured by the model. These results suggest also that as risks or costs vary in space, and as temporal variation increases – typically the conditions expected in real case studies, e.g. extreme events— it will become crucial to identify the precise time and space where to help species migrate in the landscape. In addition, other objectives or constraints in SIESTR could be used for different purposes, leading to other recommendations.
SIESTR proposes a single solution via a genetic algorithm, but we show that there is not a single solution for the optimization. In fact, multiple runs of the same models produced results that show that some optimal spatial locations for introductions may be equally likely (Figure 6), indicating that in practical applications many candidate sites may be functionally interchangeable for assisted migration. Nevertheless, some sites were consistently selected by the algorithm, suggesting that they were systematically selected, and therefore should be prioritized in planning assisted migration.
The virtual case study presented in this study represents a strong test of the framework, as both climate and cost surfaces were modeled as smooth gradients, a configuration that increases solution uncertainty. In more heterogeneous landscapes—either in species suitability or in costs—introduction choices are expected to be more constrained, thereby reducing variability in SIESTR outputs. This is likely in mountainous regions with strong spatial variability in climate and microclimates (Davis et al., 2019; He et al., 2017), in highly seasonal systems or regions affected by strong climatic perturbations such as El Niño, which generate pronounced temporal variation in colonization and extinction dynamics (Davis et al., 2016; Serra-Diaz et al., 2016a), and in landscapes with complex agroforestry mosaics or restrictive conservation policies that limit feasible introduction sites (Fahrig, 2003; Laforge et al., 2022). Although increased heterogeneity may reduce uncertainty in site selection, it can also impose ecological costs. Higher landscape heterogeneity has been shown to limit species dispersal (Serra-Diaz et al., 2015), potentially requiring a greater number of introductions to achieve a target range size.
SIESTR is based on data and assumptions that need to be carefully considered in each case. Our model inherits the same key simplifications of the KISSmig model (Nobis & Normand, 2014): the entire colonization process is described by the probability of direct spread from a presence area, thus integrating the set of complex processes of reproduction, dispersal and establishment into a single value over time. Indeed, our dispersal kernel does not vary in space and time and does not reflect variation in dispersal effort via reproduction (Ronce et al., 2000), or variation in dispersal vectors such as wind or animals (González-Varo et al., 2021; Kling & Ackerly, 2021; Pires et al., 2018). This information (albeit static) in SIESTR is considered when adjusting the migration kernel that, for instance, could be based on traits (Tamme et al., 2014). Importantly, in our example, the dispersal kernel does not consider low probability of long-distance dispersal. Long distance dispersal has been key in plants shifting their distribution in past climatic changes (Nathan, 2006), and it is still important under present conditions (Alsos et al., 2007; Baltzinger et al., 2019). We argue that in its current version, SIESTR is designed to obviate such low probability events to guide (assisted) migration. Nevertheless, migration kernels could be modified to reflect long-distance dispersal events. Beyond our example, SIESTR does not currently consider long-distance propagation beyond what is modelled by the migration kernel. Connectivity between patches is included in the calculations via the migration kernel, so with in general low distance patches. Such developments have not been made but could be the subject of further research.
Similarly, species fitness is represented using a single parameter derived from an ecological niche model or species distribution model. There is ongoing debate regarding whether outputs from such models reflect biomass growth, establishment probability, or other population-level processes (Serra-Diaz et al., 2013; Thuiller et al., 2014). In this case study, we rely on the widely used concept of climatic suitability as an indicator of overall exposure to climate change (Dawson et al., 2011), and interpret it as a proxy for population survival. Alternative metrics could, however, be incorporated into the framework. For example, population growth rates derived from integral projection models provide demographically grounded measures of fitness (Merow et al., 2014a; Merow et al., 2014b). Suitability could also be decomposed by life stage by combining multiple models, allowing different aspects of vulnerability between juveniles and adults to be represented explicitly (Serra-Diaz et al., 2016a). Ultimately, the choice of fitness or survival metric will depend on the ecological knowledge available, and the data and modeling approach applicable to the focal species.
Biotic interactions strongly influence population dynamics, but they are difficult to integrate directly into SIESTR because the model first calculates species presence probabilities at the end of the study period before choosing introductions. As a result, the environment cannot change in response to the presence or absence of our species, which can be a limiting assumption. However, these interactions can still be considered indirectly by incorporating them into the suitability matrices. As these matrices represent the species’ survival probability in each cell at each time step, users can include the effects of biotic interactions in these inputs directly, provided they know how these interactions affect survival and can predict how biotic conditions will change over time.
The method introduced here builds on heuristic approaches from spatial conservation planning, such as Zonation (Lehtomäki & Moilanen, 2013), while explicitly incorporating high-resolution spatio-temporal dynamics. This extension represents a refined assessment of the need for assisted migration beyond static suitability outputs from species distribution models (Hällfors et al., 2016). Importantly, the explicit inclusion of time enables the framework to account for temporal stochasticity and the effects of extreme events on conservation planning. This is particularly relevant given the projected increase in disturbance frequency and intensity under climate change and their influence on species’ ability to track suitable conditions (Early & Sax, 2011; Fischer et al., 2021; Liang et al., 2018, 2023; Seidl et al., 2017; Serra-Diaz et al., 2015).
In SIESTR, disturbances are implicitly represented through the suitability matrices, for example via the effects of droughts or heatwaves on establishment and survival. However, the framework does not explicitly model biotic disturbances or their spatial propagation. Drawing further from spatial conservation planning, we introduced a cost matrix to represent biological and/or anthropogenic constraints. Incorporating costs produced marked shifts in both the location and timing of introductions (Figure 4 vs. Figure 5), consistent with results from spatial prioritization studies in which alternative cost surfaces alter conservation allocations (Jung et al., 2021; Moilanen et al., 2011).
Costs in SIESTR should not be interpreted solely as economic or implementation expenses. They may also represent ecological values, such as risks to recipient ecosystems, or a combination of ecological and socioeconomic constraints. Framing costs in this way enables constraining areas such as exclusion zones, restricted areas, or prioritize regions for assisted migration. A current limitation is that costs are static through time, despite expected changes in land-use intensity and management under future socioeconomic pathways (Schrammeijer & Verburg, 2019). Such temporal shifts are already observed in conservation areas for taxa such as trees, where assisted migration initiatives are being implemented (Guo et al., 2022).
The optimization process of SIESTR has been specifically designed to explore different scenarios of species assisted migration. We provide some guidance about products and scales that could be used to build such applications with a special focus on tree species as programs of assisted migration are already underway (Appendix S3). However, the outputs from the model could be used in other applications. For instance, transition matrices can provide site-level information on the effects of planting at a given site and time. Indeed, the use of transition matrix could inform invasion biology, by finding probable initial presence sites for a current observed distribution, knowing their dispersal capabilities. That is, the optimization process and its parameters could be modified to identify sites that maximize the spread of an invasive, thus identifying high risk invasion sites beyond the use of climate similarity or raw outputs of species distribution models (Cardador et al., 2022; Redding et al., 2019). Another potential application could be developed in paleoecology, to hindcast the most likely past refugia from which species could have been spread to result in current species distributions. These widely used applications were formerly used either through static models, or through time-costly simulations.
The large effects of climate change on species redistributions require flexible methods that aid decision-making. Given the importance and yet contentious debate on introductions (either purposely or not), SIESTR will help explore the risks and potential schemes for assisted migration based on quantitative analysis. It will provide insights into effects of species introduction based on data available for many species and support new research perspectives on conservation ecology.
Acknowledgements
Preprint version 3 of this article has been peer-reviewed and recommended by Peer Community In Ecology (https://doi.org/10.24072/pci.ecology.100782, Simova, 2026)
Funding
AC acknowledges funding from pôle A2F. J.M.S.D. acknowledges the support of ANR-21-CE32-0003 and the Ramon y Cajal Program RYC2022-035668-I, funded by the Spanish Ministry of Science and Universities MCIU/AEI/10.13039/501100011033 and FSE+, and the project DISECO PID2024-162433NB-I00
Conflict of interest disclosure
The authors declare that they comply with the PCI rule of having no financial conflicts of interest in relation to the content of the article.
Data, scripts, code, and supplementary information availability
The appendices are available through Zenodo (Callebaut, 2026; https://doi.org/10.5281/zenodo.22678205).
The R package is available through Github at https://github.com/ArndCallebaut/IESTR and in Zenodo (Callebaut, 2025; https://doi.org/10.5281/zenodo.18037082). The data of the case study are built with a script available in the R package.