Schrödinger’s Range-Shifting Cat: How Skewed Temperature Dependence Impacts Persistence with Climate Change

The majority of species display strongly asymmetric responses to climatic variables, yet most analytic models used to investigate how species will respond to climate change assume symmetric responses, with largely unknown consequences. Applying a known mapping of population dynamical equations onto corresponding well-studied problems from quantum mechanics, we extend analytical results to incorporate this asymmetry. We derive expressions in terms of parameters representing climate velocity, dispersal rate, maximum growth rate, niche width, high-frequency climate variability, and environmental performance curve skew for three key responses: (1) population persistence, (2) lag between range displacement and climate displacement, and (3) location of maximum population sensitivity. We find that asymmetry impacts these climate change responses, but surprisingly, under our model assumptions, the direction (i.e., warm skewed or cool skewed) of performance curve asymmetry does not strongly contribute to either persistence or lags. Conservation measures to support range-shifting populations may have most benefit near their environmental optimum or where the environmental dependence is shallow, irrespective of whether this is the leading or trailing edge. A metapopulation simulation corroborates our results. Our results shed fresh light on how key features of a species’ environmental performance curve can impact its response to climate change.


Introduction
Climate change is driving range shifts of species across the world (Parmesan and Yohe 2003;Lenoir et al. 2020), making understanding the underlying ecological dynamics and drivers a key component of conservation efforts (Urban et al. 2016).For species with narrow environmen-tal niches, their likelihood of survival will be determined by their relative rates of extirpation from sites at trailing edges and colonization at leading edges (Kerr 2020).Range shifts are not instantaneous, creating communities not at equilibrium with their climatic niche (Svenning and Sandel 2013;Alexander et al. 2018;Rumpf et al. 2019;Lenoir et al. 2020).While species distribution models (Elith and Leathwick 2009) can indicate the potential end state after climate change, dynamic models are required to examine how a species' range may shift through time (Zurell et al. 2009;Alexander et al. 2018).To meet this need, many modeling approaches have been applied, spanning the full breadth of possible trade-offs among model precision, realism, and specificity (Levins 1966).
Large generalized spatially explicit simulation models (Brooker et al. 2007; Urban et al. 2012;Lurgi et al. 2015;Thompson and Gonzalez 2017;Thompson and Fronhofer 2019) are able to flexibly capture many processes relevant to range shift dynamics and have been highly informative.However, the complexity of these models makes systematic interrogation of conclusions drawn from particular parameter choices a challenge.Results can depend strongly on the underlying assumptions of the model (Zurell et al. 2018).Moving-habitat integrodifferential equation models (Zhou and Kot 2011;Kot and Phillips 2015;Harsch et al. 2017;Hurford et al. 2019), which include a continuous spatial element and discrete time, can capture many nuances of dispersal processes and are somewhat analytically tractable (Kot and Phillips 2015).However, conclusions are still largely restricted to inspection of simulation results.
To complement these more specific models, there remains a strong need for fundamental theory to identify critical determinants of climate change responses.Purely analytic models provide another perspective to the problem, offering an "all-else-equal" baseline for consideration.
To address this need, partial differential equation models building on well-established reaction-diffusion equation models (Cantrell and Cosner 2003) for population and gene spread (Fisher 1937;Skellam 1951;Nelson and Shnerb 1998;Hastings et al. 2005) have been applied to climate change scenarios (Pease et al. 1989;Potapov and Lewis 2004;Berestycki et al. 2009;Li et al. 2014).This work has shown the critical rate of climate change that a species can survive to be a function of dispersal rate and population growth rate from rare, allowing direct comparison with data (Leroux et al. 2013).
However, to date analytic models have been restricted to simple representations of species performance across environments.In particular, the functional dependence of the species' intrinsic growth rate on its environment-that is, the environmental performance curve (EPC)-is assumed to be symmetric.It is widely appreciated that most species show highly asymmetric environmental responses (Savage et al. 2004), for example, a gradual increase up to an optimum followed by a sharp decline.Disparities in environmental sensitivity between trailing and leading range edges may be expected to influence range-shift dynamics.For example, using a numerical approach, Hurford et al. (2019) demonstrated a case where the EPC asymmetry and dispersal interact to cause divergent effects depending on the direction of the asymmetry.Furthermore, these asymmetries may be particularly relevant when underlying climatic variability is considered (Nadeau et al. 2017).
Here, we extend previous analytic theory to incorporate asymmetric environmental dependence and directly investigate the impact of asymmetry on three key responses: likelihood of persistence, range-shift lag, and the location of peak sensitivity to conservation interventions.We then corroborate these results in a simulation exhibiting complex dynamics.

Setting and Specification of Core Model
Our approach follows the work of Nelson and Shnerb (1998) and Dahmen et al. (2000) by reformulating the species-movement problem as a Schrödinger equation, unlocking mathematical methods and results from the quantum mechanics literature.For simplicity, we outline only the mathematical derivations here.See table 1 for a summary of parameters and variables used.Full details are given in supporting information (SI) 1 (SI 1-4 are available online in the supplemental PDF), and SI 2 is an evaluated Mathematica notebook to replicate our analysis.Throughout, we assume a single spatial dimension (x) that spans a linear gradient in an environmental variable E, where initially E p hx, where the steepness of the environment gradient is h p 1 per unit length.Our E could represent any environmental variable, but it is most easily conceptualized as an average temperature variable.Climate change is introduced by a linear time dependence of E: E x,t p h(x 2 vt), where v is the rate of change in E, a reasonable assumption over medium timescales.Hence, the location where E p 0 is displaced with constant velocity v.
In a warming scenario, x could be distance from a pole or to a mountain peak and v would be negative.The correspondence between space and E means that v can also be interpreted as a "climate velocity" of dimensions of length#time 21 (Brito-Morales et al. 2018).A species' local population growth rate in a particular environment is given by an EPC function g(E).Following, for example, Fisher (1937), Pease et al. (1989), Nelson andShnerb (1998), andHastings et al. (2005), we model changes in the local population density b p b(x, t) of a species through time t at location x, including dispersal at a rate D, as The parameter c represents the sum of the strengths of all direct and indirect intraspecific competition (density dependencies).The effect of migration is incorporated via the final term of equation (1).The rate of net migration depends on the curvature of the population density to either side of the focal point.Populations near the peak of the biomass distribution (where the curvature is negative) lose population density to net migration, while the edges of the distribution (where the curvature is positive) gain.
The dispersal rate D is 0.5 times the mean squared displacement along x of a lineage per unit time and can be related directly to the mean movement of individuals per generation (Kareiva and Shigesada 1983).
Central to the approach is a change of reference frame to track the region of suitable environment across space.We define y p x 2 vt and use y as our principal spatial variable (fig.1).This allows the distribution u(y, t) p u(x 2 vt, t) p b(x, t) of the population in a reference frame comoving with the environment to be examined.Evaluation of the derivatives results in the following partial differential equation for u(y, t), equivalent to the equation for b(x, t) above: Approximate General Model Solution Using a Schrödinger Equation Building on a general approach also deployed by Nelson and Shnerb (1998), we now construct an approximate solution of equations ( 1) and (2) for the case of vulnerable populations-that is, populations that are close to extirpation.Focusing on vulnerable species simplifies the problem because when populations are low the quadratic terms in equation ( 2), describing density dependencies, will generally be small compared with both the dispersal terms and the term containing the EPC g(⋅), with the latter two mostly balancing each other (SI 1, sec.S1.3).The effect of density dependence can then be taken into account as a small perturbation.As a result, the equilibrium solution u(y, t) will, up to a scaling factor, be very similar to the form that u(y, t) would attain while the population was growing from initial low abundance (Dahmen et al. 2000;Kenkre and Kumar 2008).Hence, while the approach breaks down if the species reaches its local carrying capacity at each specific location, if dispersal losses are a major contributor and the population remains far from local carrying capacity it can be a reasonable assumption.
The approximation requires that the spatial distribution of the population prior to climate change (v p 0) is either known or has been computed from either equation (1) or (2) (SI 1).We denote this spatial distribution by the function w(y).It describes the shape of the population distribution across space without climate change and at low population density.
Mathematically, w(y) is obtained not directly as the solution of equation (2) for v p 0 but as the solution of the corresponding equation without density dependence: Here, l 0 is the maximum feasible growth rate of the population growing from low abundance without climate change, also referred to as population fitness or the invasion growth rate (Grainger et al. 2019).This equation defines w(y) up to a constant factor that does not matter for the following.It is possible to obtain from w(y) and l 0 an approximation of the species' response to climate change at any given velocity v without solving another differential equation (SI 1, secs.S1.2 and S1.3).Instead, the population's approximate distribution in the presence of climate change where U is a scaling parameter given by and Here, l is the invasion growth rate under climate change.
Under our assumptions, equation ( 4) holds up to a relative error of the magnitude of l=max y g(y).The population's actual distribution prior to climate change is a good first approximation of w(y) (pw(x)) and can be used in place of w above.Alternatively, equation ( 3) can be solved numerically.
Another possibility is to draw on the formal equivalence between equation ( 3) and the time-independent Schrödinger equation (Schrödinger 1926) for the wave function w of a particle moving in a one-dimensional energy potential V(y) (Nelson and Shnerb 1998): In this equation, E denotes the total energy, m is the particle's mass, and ħ is a constant.Comparing terms, one sees that the potential energy V(y) corresponds to the environmental performance function g(⋅) (with a sign flip).The final term describes the kinetic energy of the particle and corresponds to the dispersal term in equation ( 3).The value of these correspondences come from the fact that equation ( 7) has been studied extensively in quantum physics.This has led to a strong body of intuition about the nature of the solutions of the eigenvalue problem given by equation ( 7) and the discovery of many functional forms for V(y) for which closed-form solutions can be derived (Mattis 1993).Importantly, this approach is not adopting an interpretation of population density as a quantum mechanical wave function.Results generated from this approach are distinct from the fundamental uncertainty aspects of quantum mechanics that may be familiar to some readers, although previous authors have suggested that those could have ecological applications too (Bull 2015;Real et al. 2017).Rather, this approach allows the transfer of existing mathematical results from the quantum mechanics literature.We consider only real-valued solutions, and so the problem is closely related to classical reaction-diffusion models (Nagasawa 1993).

Representing EPCs
In our model, the function g(E) describes how the local intrinsic growth rate of a population (i.e., without dispersal effects) depends on environmental conditions.It could be any function that is negative as E becomes very large or very small.However, by selecting one of the many functions for which equation ( 7) has been solved analytically (for an accessible list, see http://wikipedia.org/wiki/List_of_quantum-mechanical_systems_with_analytical _solutions), considerable progress can be made.
An asymmetric EPC (fig.2b, 2c) can be defined in analogy to the Morse potential (Morse 1929) function where R max is the local intrinsic growth rate of the species at the environmental optimum E p f, w is a niche width parameter, and a is a nonzero asymmetry parameter.In the limit a → 0, the function simplifies to a parabolic or "harmonic" potential function (fig.2a).For simplicity, we will assume the environmental optimum to be at f p 0 and assume that jawj ! 1 throughout, which assures that both very high and very low values of E lead to negative growth rates g Mor (E).EPC g(E) provided in other functional forms are best approximated by g Mor (E) by matching the first three derivatives at the point f of maximum performance (g 0 (f) p 0).This is achieved by setting w p (2R max =jg 00 (f)j) 1=2 and a p g 000 (f)=3g 00 (f).
Care is needed when referring to skew direction-a positive value of a leads to a "negative" or left-tailed skew (fig.2b, 2c) in terms of E. With a negative v, the E variable increases with climate change-if E is a temperature variable, higher values are therefore warmer, and negative a values would be described as "warm skewed" because the heavy tail is to the warmer (more equatorial, lower elevation) side of the optimum (Hurford et al. 2019).
On top of these underlying functions, the effect of highfrequency (frequency ≫ l) temporal environmental variation can be modeled by convolving these curves with a probability density function that represents the environmental variation.When temporal environmental variability is represented by a zero-centered Gaussian distribution with standard deviation j, the convolution with the Morse potential function equation ( 8) maps onto a transformed Morse potential (SI 2).In figure 2d, we show the effect of introducing variability with j p 1 on the overall effective Range Shifts with Asymmetric EPCs 165 EPC: it softens the edges and shifts and flattens the peak of the performance curve (Ruel and Ayers 1999).

Analytic Results
We use our model to derive analytic predictions for how three response properties depend on key parameters: (1) the capacity for species to sustain themselves under climate change, (2) the lag between a species' distribution and where it "should" be if it kept pace with climate change, and (3) the location where conservation interventions would be most efficacious.

Response 1: Critical Speed of Climate Change
The sign of the invasion growth rate under climate change l is crucial for a moving population.If positive, a population at low abundance will grow (despite climate change) to the equilibrium state described by equation ( 4).If negative, the population will decline to extinction.While actual extinction events are a stochastic, not deterministic, process (Kessler and Shnerb 2009), l values can still be highly informative and provide an analytic target without the need for simulations.By equation ( 6), climate change is always detrimental to invasion growth rate.Remarkably, this invasion growth rate decline is entirely independent of the form of the EPC g(E).Also noteworthy is that the impact of climate change rate on invasion growth rate is quadratic, implying that linear extrapolations from observations of slow change will not correctly predict the impact of faster changes.
Asymmetry does, however, affect the total population size given by integrating equation (4) over y.Mathematically, this follows from the fact that w(y) itself is asymmetric and the direction of asymmetry affects the value of the denominator integral in equation ( 5).Intuitively this stems from the extent to which the skewed population distribution aligns with the skewed population growth rate.We work through an example in section 2.3 of SI 2, and this can also be seen in figure 3.
Equation ( 6) can be rearranged to determine the critical speed of climate change (v * ) at which a species can no longer keep up and will become extinct (i.e., l p 0): Alternatively, one can solve for D to identify the critical rate of dispersal that a species must exceed to maintain its population (Hastings et al. 2005;Leroux et al. 2013).
Again, neither of these results depend on the shape of g().
The shape of the performance curve does, however, influence the pre-climate-change population growth rate (l 0 ).Using the asymmetric EPC defined above to determine l 0 within equation ( 6), we can obtain (SI 2, sec.2.3) a relatively simple expression for the population fitness under climate change, This is valid as long as D ! 4R max a 24 w 22 (SI 1, sec.S1.6).The first three terms of equation ( 10) correspond to previously found conclusions about the impact of climate change in the case of a symmetric (quadratic) EPC (Pease et al. 1989).Consistent with intuition, the model predicts that greater maximum population growth rate (R max ) increases invasion growth rate and that climate change is detrimental to population fitness.The effect of dispersal rate D is multifaceted-the second term shows how greater dispersal can offset the impact of faster climate change, but the third term shows a negative effect of greater dispersal due to losses from the central part of the population.This effect is mitigated by larger niche widths.
The fourth term predicts that asymmetry (a) acts quadratically and is therefore independent of skew directionwhether the long tail in performance is on the leading or trailing range edge is not relevant to population fitness.This follows naturally for any EPC g(y) from the fact that asymmetry only influences the baseline population fitness l 0 , which does not have any sense of directionality.Most of the parameters in equation ( 10) contribute only to the l 0 component, not the response to climate change as such.
Through Jensen's inequality, temporal variance in the environment will impact the overall relationship between mean environmental values and average performance (Ruel and Ayers 1999;Vasseur et al. 2014).Specifically, high-frequency temporal variation (where average effects are representative) does not change the form of g() but modifies parameters (SI 2, secs.2.3.2 and 2.3.3).The general effect of this on l is complex (SI 2, sec.2.3.4).However, the marginal effect on l of increasing the environmental variance j 2 from a low value, is easier to interpret.Taking into account the condition D ! 4R max a 24 w 22 for equation (10), increasing variability will always lead to a reduction in growth rate.However, it shows that the marginal effect of temporal variation on invasion growth rate is dependent on a large number of model terms.The range width denominators confirm the intuition that temporal variation is most influential with narrow-ranged species.Both dispersal and asymmetry in equation ( 8) effectively widen the niche and correspondingly reduce the negative impacts of environmental variation.Again, the direction (sign) of the asymmetry is not relevant, since a enters quadratically.
Response 2: Range Shift Lag Behind Climate Change Since both dispersal and population growth take time, there is expected to be a measurable lag (in space and time) between the suitable climatic range and the distri-bution of a population (Alexander et al. 2018).With a constant rate of climate change and assuming a species' population is able to sustain itself, it will eventually reach an equilibrium in coordinates comoving with the changing climate.Even in an idealized model system like ours, there are multiple metrics that can be used to describe a species range in space (Yalcin and Leroux 2017).Here, we measure lag in the moving reference frame between the pre-climate-change peak of population density and the peak with climate change once it reaches the traveling wave stage (fig.1).
Using the asymmetric EPC including the effect of temporal variation, we can derive expressions for these lags (SI 2).Focusing on the case where the rate of climate change v is relatively small (v → 0), the steady-state lag (D) in time becomes This lag in time can be converted to lag in space D space i (fig. 1) by multiplying by the velocity of climate change.
If "lag" is instead measured in terms of the center of mass of the population's distribution, the corresponding result differs by a multiplication of the w 2 =R max term by 2 (SI 2, sec.2.4).Where there is no asymmetry (a p 0), only the first term is relevant.By the first term, lags are greater with wider niche widths-the population is not forced to respond as rapidly when a larger part of the environment is suitable.The denominator, 2 ffiffiffiffiffiffiffiffiffiffiffi ffi DR max p is sometimes referred to as the Fisher velocity because it describes the rate of propagation of a population front in a homogeneous environment with constant local invasion growth rate R max (Fisher 1937;Hastings et al. 2005).The lag in time is therefore proportional to the time it takes a population to propagate by a distance w.
From the second term of equation ( 12), it can be seen that the exact impact of a depends in a complex way on other parameters but responds only to the magnitude of asymmetry, not the sign.Overall, greater maximum growth rate R max and greater dispersal D always decrease that lag, and greater niche width w, asymmetry a, and climate variability always enhance it.The effect of climate variation j is tied to the asymmetry of the EPC, with variation combining with asymmetry increasing effective niche width, while in the symmetric case j has no effect.

Response 3: Sensitivity to Conservation Actions
Our model can be analyzed to examine where interventions are most consequential for the overall population growth rate and might therefore be of particular conservation priority.We do this by assessing the sensitivity Range Shifts with Asymmetric EPCs 167 of the overall population growth rate l to local perturbations in growth rate at particular points in space.We seek to find y Smax , the most sensitive location to perturbations in local growth rate.Note that we use the moving spatial reference frame y and optimum f p 0, so that y Smax is relative to the environmental optimum at each point in time as climate change progresses.This method is analogous to determining the most sensitive life stage (e.g., Caswell 2012Caswell , 2019)), but substituting age by space.The sensitivity of l is proportional to the local product of the reproductive value and the abundance in the stable distribution.We show in section S1.5 of SI 1 that, as a result, sensitivity is maximal at the maximum of w(y).In particular, it is independent of the rate of environmental change.For the Morse EPC, this evaluates to (SI 2, sec.2.5) For symmetric EPC, there is no shift relative to the environmental optimum at all.For our asymmetric EPC, maximum sensitivity always occurs on the long-tail side of the environmental optimum (fig.3).

Comparison with Simulations
To derive analytic results, our framework makes a number of simplifying assumptions.To test the translatability of our predictions to more complex systems, we built a spatially explicit population model.Expanding on the continuous time and continuous space analytic results, we built this model to include randomness and discreteness in both space and time.Because of the inherent challenges in determining effective parameters (Nelson and Shnerb 1998), we focus on qualitative comparisons between our analytic predictions and the simulations but show in SI 4 that quantitative predictions are of the correct order of magnitude.Where possible we use the same symbols for simulation and analytic parameters where they are similar but note that they are not directly comparable because of differences in overall model structure.

Simulation Model Specification
Full details of our simulation model are given in SI 3. In brief, we generated rectangular arenas with dimensions 40#10, each with 200 randomly located nodes connected to their nearest neighbors (fig.4a).A single environmental variable E linearly increases along the x-axis of each arena from 0 to 40.Population dynamics of each species were modeled in discrete time and followed logistic growth with dispersal between nodes, where the population growth rate depended on the current environment at that node through performance functions based on equation ( 8), with R max p 10, w p 1, and a p 50:9.To represent temporal environmental variability, for each set of 15 time steps (which we refer to as a "year") we added a value drawn from a Gaussian distribution (mean p 0, j p 0:5) to all E values during that year.Immigration into neighboring sites was determined by an exponential distance decay function, with dispersal rate chosen such that extinctions are possible once climate change is introduced.
We assembled 100 sets of species with either direction of environmental performance skew (a p 10:9 or 20.9).We did not test a symmetric EPC case, as it is not possible to standardize all aspects of the performance curve for a fair comparison.For each set, 100 species were generated with environmental optima (f) each drawn from a uniform distribution between 20 and 30, maintaining a region of suitability throughout the simulation to mitigate possible edge effects.The arenas were seeded with initial colonists and the model integrated for 200 years to fill their initial range.Species that at any point fell below a threshold biomass (10 26 ) across all nodes were considered extinct and removed.
We examined three response metrics that could be derived from both the simulation model and the analytic model, to assess whether the pattern of parameter dependencies holds.

Critical Rate of Climate Change
We sequentially simulated each assembled set of species under varying rates of climate change (v p 0 to 0.5 in steps of 0.05) and identified the fraction of species that had fallen below 1% of their starting total population size at the end of 50 years of climate change.In line with our expectations, there was a relatively precipitous decline in survival chance.Both a p 10:9 and 20.9 showed very similar responses (fig.4bi).We confirmed that the asymmetry parameter was indeed impacting v * with trials of a p 20:7, 0.3, 0.3, and 0.7 (SI 3).We fit a generalized linear mixed effects model that included as main predictors v, a, their interactions, and assemblage as a random effect to estimate v 50 , the v at which the 50% of species become extinct.When left skewed (a p 10:9), v 50 p 0:292, while when right skewed (a p 20:9), v 50 p 0:305.

Movement Lags
We tested the lag in movement of the center of mass of the population distribution for each of the assemblies during climate change.We first assessed the starting location of the "center of mass" of each species before any climate change (m i,start ) as the average x-coordinate of each node weighted by the biomass of species i at each node, averaged over 20 years.We then ran 70 years of climate change at v p 0:1 and measured the average lag in space ( D space i ) for each species over the final 20 years: Overall, a p 10:9 led to a marginally greater average spatial lag (0.618) than a p 20:9 (0.467), but this was small compared with the overall amount of variation in the simulations (fig.4bii).

Location of Peak Sensitivity
Localized conservation interventions were represented by increasing the intrinsic growth rate by two at all sites within an intervention window one spatial unit wide.This band was centered L spatial units in the x-dimension from the optimum of each species.Climate change was introduced at a rate of v p 0:28, which in the absence of intervention would cause around half the species to become extinct.When climate change was introduced, the intervention window moved with it.The simulation of 50 years of climate change was repeated with values of L ranging from 13 to 23 (in steps of 0.5), and the percentage increase in species surviving the climate change period with the conservation intervention, compared with simulations without the intervention, was recorded.
In line with the analytic expectations, we found that the location where the sensitivity had the most impact depended on the direction of the skew of the asymmetry of the EPC (fig.4biii).Interventions were most efficacious at preventing extinctions when they were located away from the optimum on the shallow environmental sensitivity (long tail of the EPC) side of the moving range, regardless of whether this was the leading or trailing edge of the range.

Discussion
Our results shed new light on how the shape of a species' environmental niche and other key drivers may impact responses to climate change.Our finding that while the extent of asymmetry is potentially highly influential, the direction of skew is not relevant to either the likelihood of population persistence or the range shift lag in our model is particularly surprising.Meta-analyses of the rates of species range shifts have not identified strong trait or taxonomic signals, instead identifying range size and habitat breadth as the most informative predictors (Maclean and Beissinger 2017; Lenoir et al. 2020).If it is assumed that the direction of environmental performance asymmetry is linked to taxonomy, this would align with our results; however, data specifying EPCs are not presently available to conduct strong tests.
In physics, beyond the simplest case of the movement of a hydrogen atom in a vacuum, the complexities of multibody problems limit the resolution with which predictions can be made.The same is true in ecology, and there are a large number of additional processes, including species interactions, evolution, demography, and environmental heterogeneity, that will influence the observed dynamics (Urban et al. 2016).While the mathematical framework of using partial differential equations in a moving environment has been validated in a laboratory bacterial system (Lin et al. 2004), there is an inevitable substantial gap between the modeled processes and known real-world ecological complexity.Analytic calculations will struggle to concurrently include the full diversity of processes that simulation models can achieve, and certain features pose particular challenges-for example, heavy-tailed dispersal models that do not have moment-generating functions (Liu and Kot 2019).Nonetheless, the expan-sion of the analytic theory is fundamental to building a foundation of expectations for how the natural world will react.
Using our model, we demonstrated that interventions would have the greatest benefit for a rare range-shifting species when located somewhat near the center of the range, where the environmental conditions are those that lead to maximal population density before environmental change set in.This region corresponds to the long tail of an asymmetric performance response.The result is valid regardless of whether this region is toward the leading or trailing edge of the at-risk species' range.This conclusion reframes discussions about whether it is most helpful to support a species at its trailing edge (where it is at risk of imminent disappearance) or toward its leading edge (where assisted translocation is possible) and could form a possible rule of thumb when information is sparse.

Applicability and Scope
In principle, the parameters in the analytic model could be directly measured for natural populations (e.g., Leroux et al. 2013).However, measuring dispersal rates remains challenging, especially given the importance of rare longdistance events (Kerr 2020).The predictability of dispersal and colonization rate is limited even under tightly controlled experimental conditions (Melbourne and Hastings 2009), leading to hard limits to the accuracy of any direct model prediction.
In our analysis, we interpreted b() as denoting a population density (individuals per unit area).However, one can instead interpret b() as the density of patches occupied (occupied patches per unit area), in the spirit of metapopulation modeling introduced by Levins (1969).This second interpretation may lead to a closer alignment between the assumptions of the model and empirical realities, as well as being more directly amenable to empirical verification or parameterization using grid-cell occupancy data.Advances in mathematical methods might further strengthen the link between the metapopulation picture and coarse-grained population distributions (Nelson and Shnerb 1998).
An important step in the derivation of our analytic results is the focus on species that are close to their extinction threshold.This allows the impact of density dependence to be assumed minimal and the application of perturbation theory.However, it is also a state where extinction risk from demographic stochasticity is significant, which can significantly impact the long-term prospects of the population (Lande 1993;Kessler and Shnerb 2009).A steady state of a system such as equation (2) can broadly be achieved in two ways-either through high dispersal and low density dependence or through higher density dependence and much lower influence of dispersal on growth rate.We examined species with minimal density dependence, and therefore dispersal terms will strongly influence local growth rates, including a not-insignificant reduction in population density in central areas due to emigration.Our model is therefore likely most directly applicable at local scales, where mass effects are more significant than occasional rare long-distance dispersal events.This context is helpful to understand why our findings contrast with a previous modeling result.Hurford et al. (2019) found through simulations of an integrodifferential equation model that positively skewed EPCs were associated with reduced lags compared with negatively skewed EPCs.Hurford et al. attributed this result to increased immigration into newly suitable areas from high-density populations near the leading edge.A key difference between our assumptions that may explain the discordance is that we focus on species close to an extinction threshold, while in Hurford et al.'s model extinctions do not occur.

Limitations and Extensions
We base our results around the properties of a stable traveling wave and use the toolbox of quantum physics to precisely describe its motion.Our analytic framework is therefore orientated around long-term equilibrium solutions.Yet there is considerable scope for the analysis of transient dynamics following the start of the perturbation (in this case, climate change), which may be more representative of available ecological observations (Hastings 2016).While our work significantly extends the scope of analytic theory in this field, there remain many further processes that are not directly considered.Our approach to modeling temporal variability effectively only modifies the underlying growth curve.It therefore only indirectly captures the impact of discrete stochastic events, which can be highly influential over short time frames relevant to contemporary climate change responses.
Furthermore, neither our analytic model nor simulations explicitly include interactions between species precluding potential "box-car effects," where competitive interactions slow the rate of advance (Urban et al. 2012;Legault et al. 2020), or the potential for extirpation from areas due to the climate-driven arrival of new species.Although the shape of the environmental performance could be attributed to indirect biotic as well as direct influences of the environment, this would not necessarily capture the complexities of interaction with other species, particularly when climate variability is considered (Terry et al. 2022).Models for interacting species within the reactiondiffusion framework have been developed (Cantrell and Cosner 2003;Potapov and Lewis 2004), and incorporating interactions within a moving-environment framework is an interesting future avenue of research.Last, our model assumes fixed species traits, but species have the potential to adapt their traits to changing climates through plasticity or evolution (Hoffmann and Sgrò 2011).Trait adaptation can be built directly into the partial differential equation framework (Pease et al. 1989;Chevin et al. 2010) and represents a further promising area for future research, particularly as nonsymmetric trait-fitness functions have been identified as an area of concern (Osmond and Klausmeier 2017;Klausmeier et al. 2020).

Conclusion
Although spatial partial differential equation models have a long pedigree within ecology (Fisher 1937;Hastings et al. 2005), our results show how rich seams of results remain to be harnessed to generate fresh ecological perspectives and more detailed baseline expectations.Our focus on the special-but critically important-case of species close to extirpation allows a simplification to an essentially linear problem and the incorporation of knowledge from other disciplines that can bring new and surprising analytical insight.

Figure 1 :
Figure 1: a, Example of movement of species range with climate change under a fixed spatial reference frame (x).Climate change rate (v) is set to 10.45 starting at t p 0. Brightness of color signifies population density; peak density is indicated by the central blue line.The red boundary delineates where the population density is above 0.01.b, Species range through time under climate change in the moving reference frame (y)-parameters are otherwise identical to a.With ongoing climate change, after a period of adjustment the population reaches a steady traveling wave state, with a lower overall population density and a constant lag behind the moving climate.Parameter values are listed in SI 2.

Figure 2 :
Figure 2: Illustration of environmental performance curves under different models that permit analytic solutions.In each case, the species has an optimum f p 0 and R max p 10. a, Harmonic potential function, w p 2 (effectively the Morse function with a p 0). b, Morse potential function, where a p 10:9, w p 1. c, Morse potential function, where a p 20:9, w p 1. d, Morse potential function but incorporating the effects of a variable, E (j p 1, a p 10:9, w p 1).Note how the peak of the curve has shifted with the introduction of this environmental variability.

Figure 3 :
Figure3: Distribution of sensitivity to intervention and population density at equilibrium with climate change.In both panels, the horizontal axis represents y, the environmental optimum is y p 0, the species leading edge is to the right, and the species trailing edge is to the left.Parameters are identical for both panels, except for the direction of asymmetry.Three quantities are plotted: the environmental performance curve (EPC) function, showing environmental optima at y p 0 (blue); the population density (orange, rescaled); and the sensitivity to intervention (green, rescaled).

Figure 4 :
Figure 4: Metapopulation simulations setting and results.a, Illustration of example virtual region of connected habitats the species inhabit in the simulation.Initial mean environmental values (E) are given by the blue (low) to red (high) coloration.The size of the yellow diamonds illustrates the population density at each site.The environmental performance curve (EPC) of the species (right skewed, a p 20:9) is drawn above.With climate change and an increasing underlying environmental variable, the species will have to shift its range leftward.b, Numerical results from the simulations, partitioned and colored by the direction of the skew of the EPC.i, Species survival as a function of the speed of climate change (v).The horizontal line shows 50% survival.Colored lines are fitted binomial generalized linear models.ii, Density plots of the lag by which species fall behind the movement of climate at a moderate rate of climate change.iii, Observed increase in survival during climate change, at different locations of conservation intervention relative to each species' optimum.Solid colored lines are generalized additive models fitted through simulation results with 95% confidence intervals.Black dashed lines illustrate the EPCs for reference.Peak efficacy aligns with the long tail of the EPC in both cases.For all three responses, the qualitative predictions of our analytic model are supported.

Table 1 :
Summary of parameters and variables used in the analytic theory * Critical rate of climate change at which the population cannot persist (i.e., l p 0) y Smax Location of peak sensitivity of l, relative to optimum D space Lag in space between total climate displacement and displacement of population at steady state D time Lag in time between changes in the environmental values and the population density Note: EPC p environmental performance curve.