Population Genetics Meets Ecology: A Guide to Individual‐Based Simulations in Continuous Landscapes

Abstract

ABSTRACT

Individual‐based simulation has become an increasingly crucial tool for many fields of population biology. However, continuous geography is important to many applications, and implementing realistic and stable simulations in continuous space presents a variety of difficulties, from modeling choices to computational efficiency. This paper aims to be a practical guide to spatial simulation, helping researchers to implement individual‐based simulations and avoid common pitfalls. To do this, we delve into mechanisms of mating, reproduction, density‐dependent feedback, and dispersal, all of which may vary across the landscape, discuss how these affect population dynamics, and describe how to parameterize simulations in convenient ways (for instance, to achieve a desired population density). We also demonstrate how to implement these models using the current version of the individual‐based simulator, SLiM. We additionally discuss natural selection—in particular, how genetic variation can affect demographic processes. Finally, we provide four short vignettes: simulations of pikas that shift their range up a mountain as temperatures rise; mosquitoes that live in rivers as juveniles and experience seasonally changing habitat; cane toads that expand across Australia, reaching 120 million individuals; and monarch butterflies whose populations are regulated by an explicitly modeled resource (milkweed).

Introduction

Explicit spatial models are indispensable for understanding how species live, interact, and evolve across geographic landscapes. However, formulating sensible models of demography in continuous space is fraught with pitfalls and choices unfamiliar to many researchers interested in spatial modeling. For instance, in the commonly used, nonspatial Wright–Fisher model, the population size is directly specified. However, in spatial models with locally defined dynamics, the number of individuals is a stochastic, emergent property. It takes some expertise to coerce the model to produce a desired equilibrium size. Even fundamental population-genetic concepts such as selection coefficients cease to have a single obvious interpretation in a spatial context.

Implementing an individual-based simulation requires great specificity—choices must be made regarding many mechanisms and parameters. Simulated organisms must separately give birth and die, unlike in more abstract theory where often only the net effects of birth minus death enter (e.g., Cantrell and Cosner 2004). Here we present modeling strategies for implementing individual-based simulations in explicit geographic space. These include spatial movement, such as dispersal, as well as spatial interactions, such as the feedback between local population density and net reproductive rate that is necessary to avoid unbounded growth.

Why individual-based simulations, and why in continuous space? Since real individuals are discrete, and live in continuous space, such simulations can in principle more accurately model many real-world situations. Some phenomena simply require individual-based simulations (e.g., locally adaptive genotypes). Non-individual-based models can be more computationally efficient, but at an often unknown cost to accuracy (Stillman et al. 2015). Similarly, discretized spatial landscapes are usually associated with model assumptions that do not provide a consistent approximation to continuous-space dynamics (e.g., Barton et al. 2002). For instance, Battey et al. (2020) showed that some aspects of genetic variation in fine grids of randomly mating demes of fixed size were irreconcilably different from continuous-space models: in fact, discretization error increased with finer grids. Given these effects of discretization, it is in our experience simpler to move directly to continuous space. Indeed, individual-based simulations in continuous space may even require less modeling expertise than other strategies, because it is usually relatively easy to come up with order-of-magnitude estimates of aspects of an organism's life cycle from natural history data, which can be used to parameterize an individual-based simulation. More abstracted modeling frameworks can involve analytical approximations, compound parameters, and other technicalities that require careful checking and mathematical experience. For instance, it is more difficult to translate “offspring disperse around 100 m” to a migration rate between adjacent grid cells than to directly input a mean dispersal parameter of 100 m. Although abstract models with fewer choices might feel more general for theoretical work or methods development, this can be a false generality, as more work is required after the fact to determine to which real-world organisms a given result or method is applicable.

Explicit individual-based population models are not new to ecology (DeAngelis and Yurek 2017; Grimm 1999). A great deal of ecological work has sought to quantify the effects of density-dependent demographic feedback. For a single species, negative feedback between population density and growth rate is necessary to avoid the population growing without bound, although there are a great many ways to set this up in practice (De Wit 1960; Beverton and Holt 1957; Ellner et al. 2016). Demographers have a deep understanding of how to describe and parameterize the statistical properties of birth and death, and what the emergent consequences are for population growth, lifespan, age distribution, and long-term fitness (Tuljapurkar 2013). Although temporal stochasticity is relatively well understood in demography (Tuljapurkar 1989), the consequences of spatial heterogeneity—particularly outside of metapopulation models (Hanski 1997)—have received less attention.

Geography can have strong effects on patterns of genetic variation (Wright 1943; Malécot 1969; Rousset 2000; Charlesworth et al. 2003; Battey et al. 2020; Min et al. 2022) and on evolutionary processes (Felsenstein 1976; Uecker et al. 2014; Savolainen et al. 2007). Genetic differentiation is shaped by the movement of individuals, and hence distance, geographical features, and the spatio-temporal history of the species (Hewitt 2011; Rosenberg et al. 2005; Ramachandran et al. 2005). Geography is, therefore, not only important, but also a relatively untapped source of information to inform inference (Bradburd and Ralph 2019). However, it is difficult to obtain analytical predictions from spatial population genetics models (Felsenstein 1975; Barton et al. 2002). Most spatial work in population genetics uses partial differential equations that do not represent genetic differentiation (Montgomery 1973; Barton 1979; Sedghifar et al. 2016; Etheridge et al. 2024), or specifically looks at the fronts of population expansions (e.g., Barton et al. 2013; Paulose et al. 2019; Nullmeier and Hallatschek 2013; Etheridge and Penington 2022).

Simulation has long been a useful tool in the study of populations (Grimm 1999; DeAngelis and Yurek 2017), particularly for the purposes of prediction (e.g., for population viability analysis, see Dunning Jr. et al. 1995)—even predating the common usage of digital computers (Pearson 1960). Simulations are also useful for inference, ranging from exploratory studies to training for deep learning. Their use depends on their computational cost: many modern machine-learning methods are only feasible with relatively fast simulations. Introducing geographic space increases computational complexity, making the task of producing training data more difficult (but see Smith et al. 2023; Champer et al. 2021). Today, there are several sophisticated software suites targeted at just the sort of landscape-scale simulations we discuss here, with both discretized spatial models (Landguth et al. 2017; Bocedi et al. 2021; Landguth and Cushman 2010; Schumaker and Brookes 2018; Neuenschwander et al. 2018; Rebaudo et al. 2013) and continuous spatial coordinates (Haller and Messer 2023).

How to Use This Paper

This paper is intended as a guide to the territory of spatial modeling, with a focus on individual-based simulations in discrete time in continuous geographic space.

The paper has two broad parts: The first part (Sections 2 through 8) describes modeling strategies to capture ecological and geographic processes in individual-based models. The second part, Case Studies, provides concrete examples of how to combine these modeling strategies to simulate a biological system.

First, we tackle population regulation (Section 2), as this is the topic that in our experience is the greatest barrier to new modelers. We then describe our strategy for parameterizing and analyzing spatial simulations (Section 3). Next, we cover spatial movement and mate choice (Sections 4 and 5, respectively), highlighting some interesting challenges. Subsequently, we discuss how to implement spatial heterogeneity and further concepts in population regulation, including stochasticity (Section 6), and finally natural selection (Section 8).

We hope this material will be useful to researchers who wish to produce simulations that are—at least roughly—modeled on concrete empirical systems. This includes empirical researchers wanting to explore the plausibility of hypotheses or the power of study designs, methods developers who want to test their methods on realistic spatial models, theoreticians wishing to explore spatial models, and managers wishing to explore alternative scenarios. So, we aim to make it easy to build a model starting from those quantities that we generally have good estimates of for particular organisms. This differs from the modeling philosophy in much theoretical work, which often begins from a “simplest possible” model.

To give practitioners a head start on simulating the models we describe, the main text is accompanied by Boxes containing code that implements each modeling concept. The code may be run with the program SLiM (Haller and Messer 2023), a flexible and powerful individual-based eco-evolutionary simulator that, as of the newly released version 4.2, includes a full set of tools for modeling interactions between multiple species across geographic landscapes. We recommend that readers open the minimal SLiM template (provided in Appendix A and at ) in SLiM's GUI to experiment with while reading. The online code repository also contains SLiM scripts for each of the Case Studies.

Population Dynamics and Density

We start by describing how to maintain a stable population. Consider an extremely simple spatial simulation: organisms are asexual, and do not move during their lifespan. Each time step, each adult gives birth to a Poisson($$$ f $$$) number of offspring, that each disperse to a random location, whose displacement from their parent's location is drawn from a Normal distribution with mean 0 and standard deviation $$$ {\sigma}_D $$$. Then, each individual dies with a probability $$$ \mu $$$.

Simulating this model, its flaw is immediately apparent as it runs: either all individuals rapidly disappear (if $$$ f\le \mu $$$), or the computer grinds to a halt as the number of individuals explodes (if $$$ f>\mu $$$). We need some kind of density-dependent feedback to maintain a stable population size—in other words, we need the net population growth rate to change from positive to negative as the population density grows past some point. When population density in an area is high, either birth rates need to decrease or death rates need to increase.

To provide density-dependent feedback, we first need a notion of “population density” at a point in space, a measure of the number of individuals nearby per unit area. Let $$$ K $$$ denote the desired equilibrium population density (in individuals per unit area). A general way to define $$$ n\left(\boldsymbol{x}\right) $$$, the population density around a location $$$ \boldsymbol{x} $$$, is to specify an interaction kernel, $$$ \rho \left(\boldsymbol{x}\right) $$$, which is a nonnegative function with $$$ \int \rho \left(\boldsymbol{x}\right)\mathbf{d}\boldsymbol{x}=1 $$$; an interaction scale $$$ {\sigma}_X $$$; and then if the two-dimensional locations of the individuals in the population are $$$ {\boldsymbol{x}}_{\mathbf{1}},\dots, {\boldsymbol{x}}_{\boldsymbol{n}} $$$, define $$$ n\left(\boldsymbol{x}\right) $$$ by1$$$ n\left(\boldsymbol{x}\right)=\frac{1}{\sigma_X^2}\sum \limits_i\rho \left(\frac{\boldsymbol{x}-{\boldsymbol{x}}_{\boldsymbol{i}}}{\sigma_X}\right). $$$

Since $$$ \int \rho \left(\boldsymbol{x}/{\sigma}_X\right)\mathbf{d}\boldsymbol{x}={\sigma}_X^2 $$$, the value of $$$ n\left(\boldsymbol{x}\right) $$$ is in units of individuals per unit area. A common choice for $$$ \rho \left(\boldsymbol{x}\right) $$$ is the Gaussian density function. One concrete interpretation is that if $$$ \rho \left(\left(\boldsymbol{x}-\boldsymbol{y}\right)/{\sigma}_X\right) $$$ gives the proportion of time that an individual at $$$ \boldsymbol{y} $$$ spends near $$$ \boldsymbol{x} $$$, then $$$ n\left(\boldsymbol{x}\right) $$$ is proportional to the time spent by all individuals near $$$ \boldsymbol{x} $$$. (More concretely, $$$ {\int}_An\left(\boldsymbol{x}\right)\mathbf{d}\boldsymbol{x} $$$ is the total amount of time spent by all individuals in the region $$$ A $$$).

Now suppose that in each time step of the model, each individual has a chance to produce offspring. The number of offspring depends on the individual's location. An individual at $$$ \boldsymbol{x} $$$ produces a random number of offspring with a mean of $$$ f\left(n\left(\boldsymbol{x}\right)/K\right), $$$ where $$$ f $$$ is the birth rate or fecundity. Each juvenile then disperses to a nearby location whose displacement from $$$ \boldsymbol{x} $$$ is chosen from a given probability distribution. Then, all individuals (including those just born) survive to the next time step with probability $$$ 1-\mu \left(n\left(\boldsymbol{x}\right)/K\right) $$$, where $$$ \mu (u) $$$ is the mortality rate at scaled density $$$ u $$$. This way, birth and death rates depend on an individual's location $$$ \boldsymbol{x} $$$ via the smoothed population density $$$ n\left(\boldsymbol{x}\right) $$$, scaled by a parameter $$$ K $$$ that controls the equilibrium density. We build on this general form throughout the paper to describe how life history, geography, and selection make individual fecundity and mortality rates more complex.

Equilibrium Population Density

Will this model stabilize, and if so, to what density? We expect an equilibrium when births balance deaths. If we define the local per-capita net reproductive rate (the expected increase due to birth minus the decrease due to death, per individual) to be2$$$ F(u)=f(u)\left(1-\mu (u)\right)-\mu (u). $$$then we expect an equilibrium at a density of $$$ {n}_{\ast } $$$, solving $$$ F\left({n}_{\ast }/K\right)=0 $$$. Note that newborns are subjected to the same mortality rate $$$ \mu (u) $$$ as the rest of the population, so the expected increase is $$$ f(u)\left(1-\mu (u)\right) $$$ rather than $$$ f(u) $$$ in Equation (2). For this reason, it is convenient to choose the functional forms of $$$ f $$$ and $$$ \mu $$$ so that $$$ F(1)=0 $$$, in which case we expect the equilibrium $$$ {n}_{\ast } $$$ to be roughly equal to $$$ K $$$. (However, we will see in Sections 6, 7.3, and Appendix B.4 that often $$$ K $$$ is not exactly the equilibrium density.) In addition to $$$ F(1)=0 $$$, for the population to be stable we also need $$$ {F}^{\prime }(1)<0 $$$ so that the net reproductive rate decreases with density near to the equilibrium. The argument $$$ u $$$ in $$$ f(u) $$$, $$$ \mu (u) $$$, and $$$ F(u) $$$ is the scaled population density $$$ u=n\left(\boldsymbol{x}\right)/K $$$, written this way so that it is easy to control the equilibrium population density by simply changing $$$ K $$$, independently of other factors.

A brief note on what we have done here: it might seem most natural to set up a simulation using fecundity and mortality rates based on empirical observation. If so, then mean population density would be an emergent property of the simulation—in other words, our script would not have a parameter $$$ K $$$ that we could directly adjust. However, to do this we need empirical estimates of how fecundity and mortality depend on local density, which are very difficult to obtain. A much more common situation is to have estimates of fecundity and mortality rates and population density at equilibrium. The parameterization we outline here, in which $$$ K $$$ is a directly tunable rather than emergent parameter, is designed to make this use case natural.

Functional Forms

The next question is: what forms should we use for the birth and death rate functions $$$ f(u) $$$ and $$$ \mu (u) $$$? There is surprisingly little guidance from the theoretical or empirical literature. Population models mostly come in two flavors: “phenomenological” (or “top-down”) models that only consider the net reproductive rate, $$$ F(u) $$$; and “mechanistic” (or “bottom-up”) models that explicitly consider birth and death separately (Geritz and Kisdi 2012). We need a mechanistic model, since our approach here is individual-based, but most literature on density dependence uses phenomenological models (for notable exceptions, see Coulson et al. (2008), or integral projection models Ellner et al. (2016)). Matrix population models (Caswell 2000) are (mostly) mechanistic and widely used for management, but rarely incorporate density dependence. Estimating the functional forms of density dependence from empirical data is a logistically and statistically daunting task, even without considering environmental variation. As we intend to simulate from our model, we will simply choose a mathematically convenient form, reducing the problem to estimating parameters given that functional form. Perhaps the most obvious strategy is to multiply the base fecundity or survival rate by a function of the local density that decreases when density is high. Appendix C works through examples using common functional forms for $$$ F(u) $$$.

Regulation by Mortality

Suppose we would like to use density-dependent feedback only on mortality. In this case, the average number of offspring is constant: we may say $$$ f(u)=f $$$. Rearranging Equation (2), we see the survival probability in a location experiencing a scaled population density $$$ u $$$ is3$$$ 1-\mu (u)=\frac{1+F(u)}{1+f}. $$$

The Beverton–Holt form for the net effects of density (Beverton and Holt 1957) would have that $$$ F(u) $$$ is proportional to $$$ \left(1+a\right)/\left(1+ au\right)-1 $$$, for some constant $$$ a $$$ that controls the strength of the density-dependent feedback. So, to set up the model to have “Beverton–Holt” feedback, we plug in $$$ F(u)=\left(1+a\right)/\left(1+ au\right)-1 $$$ and obtain that4$$$ 1-\mu (u)=\frac{1+a}{\left(1+f\right)\left(1+ au\right)}\kern0.48em \mathrm{if}\ a\le f. $$$

The leftmost plot in Figure 1 visualizes Equations (3) and (4); the others are discussed in Section 7. Mortality regulation is also demonstrated in Box 1, and most of the examples in this paper use this model with $$$ a=f $$$. In practice, an empirical estimate of the survival probability at low density (the limit as $$$ u\to 0 $$$) $$$ {s}_0 $$$ can be incorporated so that $$$ a={s}_0\left(1+f\right)-1 $$$.

[IMAGE OMITTED. SEE PDF]

The condition that $$$ a\le f $$$ is so that this gives us a valid probability: otherwise, the survival probability can be negative. Note that the strategy used here cannot work with the discrete logistic: since $$$ F(u)=1-u $$$ can get arbitrarily negative, eventually random density fluctuations will lead to negative probabilities, and errors in the simulation.

There are many other possible ways to set up $$$ f(u) $$$ and $$$ \mu (u) $$$ that result in the same net form of density dependence. (See Section 7.3 and Appendix C for how stochasticity alters the realized equilibrium density and alternative functional forms.) The general forms of $$$ f(u) $$$ and $$$ \mu (u) $$$ may also need to be modified depending on the organism. If there are two sexes and only one sex can bear offspring, then $$$ f $$$ should not be the total number of offspring, but rather the number of offspring of the offspring-bearing sex. More generally, if vital rates depend on any aspect of the individual (e.g., only adults can reproduce, or mortality is age-dependent), then the equivalent calculations must be done with a matrix population model (Caswell 2000) or an integral projection model (Ellner et al. 2016). Our examples mostly ignore such complications, although several models have distinct life stages (Case Studies 9.2 and 9.4).

1

BOX

Regulating mortality and fecundity in SLiM.

The Wright–Fisher model is a population model with a fixed population size. A Wright–Fisher model in continuous space, therefore, has global population regulation—individuals are affected by others arbitrarily far away, resulting in counter-intuitive, unrealistic consequences (Felsenstein 1975). SLiM's default model is the Wright–Fisher model (“WF”). In this study, our goal is to model particular species in realistic ways, so we use only SLiM's non-Wright–Fisher (“nonWF”) model, which requires explicitly choosing how births and deaths occur.

In a SLiM nonWF model, the “fitness” attribute of an individual is the probability of survival until the next time step. So, controlling mortality in SLiM is simply a matter of setting individuals' fitness. To compute local density, we use the localPopulationDensity() function, which computes density just as described in Equation (1) (with options for the choice of kernel). To do this we need to first (during setup) define an InteractionType object, here using a Normal kernel with a scale of SX (as in, $$$ {\sigma}_X $$$) and a maximum distance of (SX * 3):

1 initializeInteractionType(1, “xy”, maxDistance=SX * 3);2 i1.setInteractionFunction(“n”, 1.0, SX);

Now, to use Beverton–Holt regulation on mortality as described around Equation (4) (setting $$$ a=1 $$$ in Equation (4) and with $$$ K $$$ defined elsewhere), in each time step, we use this interaction type to set survival probabilities:

3 inds = p1.individuals;4 density = i1.localPopulationDensity(inds);5 u = density / ((1 + f) * K);6 inds.fitnessScaling&#x02009;=&#x02009;1 / (1 + f * u);

Note that the local population density is measured after reproduction and before death ($$$ {N}_t^{+} $$$ in Appendix B), so we divide by $$$ 1+f $$$ as well as $$$ K $$$ when converting it to scaled population density $$$ u $$$ from the Equation (1), so it reflects density before reproduction.

On the other hand, offspring are produced by a “reproduction() callback,” a chunk of code that is executed each time step for each individual, and produces any desired new offspring for the focal individual. For instance, if we would like to use Ricker regulation on fecundity (see Equation (19) with $$$ \alpha =1 $$$ and $$$ \beta =0 $$$), we might for efficiency pre-compute each individual's number of offspring:

7 inds.tag = rpois(length(inds), lambda=exp(-density/K));

and then in the reproduction() callback produce the offspring:

8 p1.addCrossed(individual, mate, count=individual.tag);

Spatial Scales and Neighborhood Sizes

Now that we understand how to obtain stable simulations, we can move on to a more spatial topic: that of spatial scale. It is convenient to use spatial scale parameters to describe distance within the model; some are depicted in Figure 2: (i) interaction scale $$$ {\sigma}_X $$$, the typical distance over which individuals affect each other ecologically, (ii) dispersal scale $$$ {\sigma}_D $$$, the typical distance between parent and offspring, (iii) mate choice scale, $$$ {\sigma}_M $$$, the typical distance between mates, and (iv) movement scale $$$ {\sigma}_V $$$, the typical displacement of an individual each time step. ($$$ {\sigma}_D $$$ is often called “dispersal distance,” but here we use “dispersal scale” for consistency.) We met $$$ {\sigma}_X $$$ in Equation (1), where it determined which individuals were close enough to each other to affect local density.

[IMAGE OMITTED. SEE PDF]

Each spatial scale parameter determines how many other individuals typically exist in an individual's “neighborhood.” First, the mating neighborhood size $$$ {N}_M=4{\pi \sigma}_M^2K $$$ and the interaction neighborhood size $$$ {N}_X=4{\pi \sigma}_X^2K $$$, measure, respectively, the typical number of other individuals in a circle of radius $$$ 2{\sigma}_M $$$ or $$$ 2{\sigma}_X $$$ at equilibrium.

In practice, these measures are only intended to be order-of-magnitude diagnostics. For example, $$$ {N}_M $$$ should be the same order of magnitude as the number of potential mates; if it is small then mate limitation may be a problem. That said, depending on the mating kernel (i.e., shape of the distribution that determines how “attractive” nearby mates are), individuals as far away as $$$ 3{\sigma}_M $$$ are probably also available to mate (See Section 5).

Less obvious but equally important is $$$ {N}_X $$$, which measures the typical number of other individuals that “count toward” the local density of a given individual, and hence affect its demographic rates (such as survival probability, as in Box 1). $$$ {N}_X $$$ also provides a measure of individual-to-individual variability in local density. Roughly speaking, a larger $$$ {N}_X $$$ means that local density is obtained by averaging over an area with more individuals, meaning that individuals across the landscape experience more similar local densities. When $$$ {N}_X $$$ is small, densities across the landscape can be more noisy—with some individuals experiencing no competition from neighbors, and some experiencing high density. Further discussion of how $$$ {N}_X $$$ and $$$ {N}_M $$$ can help diagnose odd model behavior is given in Appendices B.2 and B.3.

Another “neighborhood size” is a classical one: Wright's neighborhood size $$$ {N}_W=4{\pi \sigma}_{\mathrm{eff}}^2K $$$, where $$$ {\sigma}_{\mathrm{eff}}^2 $$$ is the (squared) effective dispersal distance, the variance of the displacement between parent and offspring along any axis looking back along a lineage (chain of parent–offspring relationships). This quantity appears frequently in work on continuous spatial models in population genetics (e.g., Wright 1943, 1946; Barton et al. 2002; Rousset 1997; Robledo-Arnuncio and Rousset 2010), and is clearly affected by $$$ {\sigma}_D $$$, $$$ {\sigma}_V $$$, and $$$ {\sigma}_M $$$, as well as the mean generation time if adults move, but no explicit expression for $$$ {\sigma}_{\mathrm{eff}} $$$ from these parameters is known. $$$ {N}_W $$$ gives, roughly, the number of “potential parents” of a given individual, and so is a measure of the rate of local genetic drift. If $$$ {N}_W $$$ is small, then local inbreeding (and spatial structure more generally) will be stronger.

Although we encourage basing modeling decisions on empirical understanding, it is not always feasible. For instance, suppose one is simulating a species of fairly common shrub that lives widely across a landscape. In practice, its local density is determined by microhabitat and complex interactions with other species. Any concrete estimate of the interaction scale for a plant is probably quite small—the competitive effect of one shrub on another more than a few meters away is (in the short term at least) usually quite small. However, implementing a spatial model with an interaction scale of only a few meters (and no other species) will likely lead to a population size that is much too large. One option is to somehow include other species and fine-scale habitat suitability, but doing this in a realistic and efficient way can be a major challenge. A simpler option is to set the interaction scale to be on the order of the mean interindividual spacing and adjust the form of density dependence to roughly match the observed population density. The simulated population will probably be more evenly spread out across space than in reality, but it is hopefully at least a better approximation than a nonspatial model. More work is needed to develop appropriate modeling strategies for such situations and to understand their consequences.

2

BOX

Summarizing the state of the population.

SLiM's GUI lets the user visualize the state of the simulation as it unfolds. We can customize the display to see areas of higher fecundity, differences in age structure, or even local adaptation.

One strategy is to set the color of each individual: by default, individuals are colored by survival probability, but this code snippet will set the color to reflect the sex of each individual:

9 females = p1.subsetIndividuals(sex = “F”);10 females.color = “white”;11 males = p1.subsetIndividuals(sex = “M”);12 males.color = “red”;

Another strategy is to summarize the state of the population as a map. This can be done with the summarizeIndividuals() function, which creates a rasterized map for which the value of each pixel is some summary of the individuals within that pixel.

For instance, the map of density shown as the background was made using the following code:

13 density = summarizeIndividuals(p1.individuals, c(25, 25), p1.spatialBounds,14 operation=”individuals.size();”,15 empty=0.0, perUnitArea=T);16 p1.defineSpatialMap(“density”, “xy”, density, T, c(0,K), colors(20, “viridis”));

In this close-up screenshot of a dioecious simulation, there tend to be more females (white points) where the local density is higher (yellow background); and there tend to be more males (red points) where the density is lower (blue background). We discuss the heterogeneous spatial patterns generated in dioecious simulations in Section 5. See Box 5 for more use of defineSpatialMap().

Movement and Dispersal

How individuals move across the landscape influences both their spatial distribution and how related they will be to nearby individuals. We consider parent–offspring dispersal and the movement of organisms during their lifetimes as different processes, but they can be implemented in similar ways. To model dispersal and movement, we need to consider both the overall scales of movement ($$$ {\sigma}_D $$$ and $$$ {\sigma}_V $$$; introduced in Section 3), as well as the shape of the dispersal distribution, or kernel.

The easiest way to implement dispersal in two dimensions is simply to say that an individual at location $$$ \boldsymbol{x}=\left(x,y\right) $$$ will produce an offspring that lives at $$$ \left(x+{\sigma}_DX,y+{\sigma}_DY\right) $$$, where $$$ X $$$ and $$$ Y $$$ are independent draws from a Gaussian. Most simulations in this paper were done in this way, for familiarity rather than any particular reason.

3

BOX

Random dispersal and displacement.

Perturbing a spatial location by adding a displacement drawn from a given kernel is a common operation in simulations—for instance, to choose offspring locations. To do this, SLiM provides the pointDeviated() method, which also needs to know the shape of the kernel, the type of boundary condition, and other parameters. For instance, the following code:

17 locs = subpop.pointDeviated(nOff, individual.spatialPosition, “reprising”,18 INF, “t”, DF, SD);19 offsprings = subpop.addCrossed(individual, mate, count=nOff);20 offsprings.setSpatialPosition(locs);

creates nOff offspring, and sets each offspring's position to a randomly sampled location near the location of the parent (individual). The random displacement is drawn with density $$$ f(x)\propto {\left(1+\parallel x/\mathrm{SD}{\parallel}^2/\mathrm{DF}\right)}^{-\left(\mathrm{DF}+1\right)/2} $$$, that is, from the kernel whose density is formed by rotating the $$$ t $$$ distribution with scale SD (as in $$$ {\sigma}_D $$$) and DF degrees of freedom about the origin. (See Appendix D for other strategies.) The reprising argument conditions the result on falling within the spatial bounds of the simulation—other options include stopping and reflecting. (Or for an absorbing boundary, none, with offspring falling outside the area removed).

Although the scale of movement most strongly affects spatial patterns, the shape of the kernel is also important, and can have surprising effects: even very rare long-range movement can have strong effects on the speed of a range expansion (Mollison 1972; Paulose et al. 2019) or the relationship between genetic and geographic distances (Smith and Weissman 2023). To see the effects of rare, long-range movement, a convenient “fat-tailed” kernel is the Student's $$$ t $$$: the smaller the degrees of freedom parameter, the more likely are extremely long movements.

How to draw from a different dispersal kernel? The first guess—choose $$$ X $$$ and $$$ Y $$$ from a different distribution—does not work: the result will not be rotationally symmetric, and dispersal will tend to fall along the $$$ x $$$ or $$$ y $$$ directions. Options are to either move a random distance at a uniformly chosen angle (which is conceptually simpler), or to multiply a bivariate Gaussian by a random scaling factor (which has other advantages). What to call a given two-dimensional kernel is not standard—for instance, would a “Student's $$$ t $$$ kernel” have a $$$ t $$$-distributed distance? Or, a $$$ t $$$-shaped cross-section? In simulations below we use the latter convention, as described in Box 3, and discuss these choices more in Appendix D.

The behavior of individuals near boundaries must of course be specified. Common choices include “reflecting” and “absorbing”. These differ substantially: imagine the new individual as taking a straight-line path from $$$ \left(x,y\right) $$$ in the direction of $$$ \left(X,Y\right) $$$ and encountering a boundary either bounce off (reflecting) or die (absorbing). The latter clearly reduces the effective fecundity of individuals near the boundary, and so reduces mean density up to a few multiples of $$$ {\sigma}_D $$$ away. Another choice is “reprising”, for which the random draw $$$ \left(X,Y\right) $$$ is chosen conditionally so that $$$ \left(x+\sigma X,y+\sigma Y\right) $$$ stays within the range. (For more details of the consequences of boundary conditions, see Mazzucco et al. 2018).

Movement in practice often depends on the environment, of course: organisms tend to move within particular habitats, and barriers are ubiquitous on all scales. Small-scale heterogeneity may be averaged out across the time scale simulated, and so incorporated (implicitly) in the movement kernel. However, large-scale heterogeneity can be important. Movement on a heterogeneous landscape still relies on some way of randomly choosing nearby points, and hence a movement kernel. One way to incorporate this is discussed in Box 5.

Clumping

In many species, reproduction is local, and so it tends to produce clumps that are spread out by movement and dispersal. Such patterns can be intriguing or puzzling and have real consequences for demography and genetic variation. Spatial clumpiness is sometimes visually obvious (as in Figure 3B), but more generally, the clumping tendency of individuals can be measured by spatial correlations.

[IMAGE OMITTED. SEE PDF]

One informative measure of clumping is the pair correlation function, which shows for each distance $$$ x $$$ how likely an individual is to have a neighbor at that distance relative to the average density. (See Baddeley et al. (2015) other useful descriptors of point data.) Concretely, the pair correlation function estimates the mean density at distance $$$ x $$$ away from an individual divided by the overall mean density, averaged across individuals; if the points are independently placed, it is constant at 1.0. Figure 3C shows pair correlation functions for simulations in which adults do not move and dispersal follows either a $$$ t $$$ or a Gaussian distribution, or in which both dispersal and adult movement are Gaussian (with $$$ {\sigma}_D={\sigma}_V $$$).

Individuals are more likely to be around a distance of $$$ {\sigma}_D $$$ from each other than otherwise, but this tendency is reduced with more long-distance dispersal ($$$ t $$$ dispersal with lower degrees of freedom; Figure 3C). In other words, a little long-range dispersal reduces clumping. Note, however, that in these examples the scale of clumping is quite narrow: correlations only extend out to 2 or 3 multiples of $$$ {\sigma}_D $$$. For another visualization of this relatively subtle clumping, see Figure S3. Unsurprisingly, adding adult movement ($$$ {\sigma}_V\ne 0 $$$) to a dispersal-only model reduces correlations as well.

In spatial models, neighbors tend to be more related to each other than to distant individuals, a pattern known as “isolation by distance” (Wright 1943). The scale at which this correlation appears is determined by how far individuals move, and is also affected by the shape of the dispersal kernel: Figure 3C,D show that simulations with more long-range movement—but with comparable mean dispersal scale—tend to have a weaker relationship between geographic and genetic distance. (See also Smith and Weissman 2023).

Mating and Other Pairwise Interactions

Mating is a crucial interaction for biological simulations, and there are numerous aspects and choices to consider. We probably don't want to simulate the detailed movement of individual pollen grains or the meanderings of a male moth seeking a female, and instead would like to skip to the realized outcome, that is, “choose a mate nearby.” (We follow the literature on mating systems in calling this “mating,” even when referring to plants or broadcast spawners).

As with dispersal (Section 4), it is easiest to specify what “nearby” means with a kernel: roughly, an individual can be chosen with probability proportional to the kernel. More concretely, if the kernel is $$$ \rho $$$ and the mating scale is $$$ {\sigma}_M $$$, then an individual at $$$ \boldsymbol{x} $$$ would assign weight $$$ {w}_i=\rho \left(\left({\boldsymbol{x}}_i-\boldsymbol{x}\right)/{\sigma}_M\right) $$$ to another individual at location $$$ {\boldsymbol{x}}_i $$$. The probability individual $$$ i $$$ is chosen is equal to $$$ {w}_i $$$ divided by the sum of weights across all nearby individuals. This same method can be used for other individual-to-individual interactions, such as predators choosing prey individuals.

We can adjust the typical distance between mates with the mating scale, $$$ {\sigma}_M $$$, which is directly analogous to parameters used to describe interaction ($$$ {\sigma}_X $$$) and dispersal ($$$ {\sigma}_D $$$) scales (Section 3), although each uses an underlying kernel in slightly different ways. Of course, the criteria determining potential mates for a given individual differ widely among species (Shuster 2003). Relevant questions about the mating system include: How often does selfing occur, and under what circumstances? Are sexes separate (dioecy/gonochory) or not (monoecy/hermaphrodity)? Are there distinct mating types or self-incompatibility systems? An important note is that for dioecious species, calculations to determine stability of population density are easiest if done using only the reproducing sex.

4

BOX

Interactions, and mating.

The mechanism that SLiM uses to mediate most effects that some individuals have on others is called an “interaction type” (see Box 1). We used a symmetric interaction type in Box 1 to compute local density: every individual affected every other. Some interactions are not symmetric: we might, for instance, want each female to be able to find nearby males. To do this, we first set up a sex-specific interaction (again using a Gaussian kernel, with standard deviation SM as in $$$ {sigma}_M $$$):

21 initializeInteractionType(2, "xy", maxDistance=SM * 3, sexSegregation = "FM");22 i2.setInteractionFunction("n", 1.0, SM);

The sexSegregation parameter value of “FM” means that females will receive the interaction and males will exert the interaction; it is asymmetric. Then, we can use the interaction in a reproduction(NULL, “F”) block as follows (the NULL and “F” arguments imply that it applies to all females) to produce a single offspring:

23 mate = i2.drawByStrength(individual, 1);24 subpop.addCrossed(individual, mate, count=1);

This chooses a male mate from among the neighbors of the focal female, individual, with probability proportional to a Gaussian kernel with standard deviation SM. The interaction type itself guarantees that the chosen mate will be male. It is possible to set up other constraints on interaction types as well, such as a minimum or maximum ages, to represent other constraints on the reproductive eligibility of individuals.

Sex-Specific Spatial Structure

The combination of density dependence and mating system can have surprising consequences: for instance, in dioecious simulations, clustering is dependent on sex. To illustrate this, Figure 4 shows the mark connection function (Baddeley et al. 2015, §14.6.4.2) for female/male, female/female, and male/male neighbors, which shows the proportion of pairs of points at distance $$$ r $$$ that are one female and one male, two females, or two males, respectively. Curiously, we see that the probability that an individual's neighbor is of the other sex does not depend on that neighbor's distance. However, the probability of having a same-sex neighbor does change with distance. As shown in Figure 4, within $$$ 3{\sigma}_D $$$ neighbors of a female are more likely than expected to be female, and nearby neighbors of a male are less likely to be male than expected.

[IMAGE OMITTED. SEE PDF]

Why does this happen? A simple reason is that we are modeling a dioecious scenario with no adult movement and where offspring are only generated by females and placed nearby. Local density-dependent mortality means that all individuals tend to kill their neighbors, but only females can replace them. Correspondingly, the spatial range of a dioecious system where offspring disperse from a particular sex is determined by the range of that sex—others on the periphery cannot extend the range because offspring do not disperse from them. This is a disadvantage for colonizing new areas (Obbard et al. 2006) and may explain the spatial distribution of dioecious individuals away from a range front (Mirski et al. 2017).

Indeed, Shuster (2003) showed that female aggregation (as observed in Figure 4) is a universal consequence in mating systems with female choice (Shuster 2003, ch. 2). Other aspects of the simulation may differ by sex: for instance, sex-biased dispersal is pervasive (Trochet et al. 2016) and generates detectable spatial patterns of genetic relatedness between sexes (Aguillon et al. 2017; Broquet and Petit 2009; Laporte and Charlesworth 2002). Social mating structures can induce sex-biased dispersal and thereby create similar patterns (Pusey and Packer 1987; Hammond et al. 2005).

Maps: Spatial Heterogeneity

5

BOX

Defining and manipulating maps.

SLiM provides support for defining and manipulating spatial maps. We can read the values for a spatial map (say, of elevation) from a .csv file containing a rectangular grid of values using the code below:

25 mapValues = readCSV(“elevation.csv”).asMatrix();26 map = p1.defineSpatialMap(“elevation”, “xy”, mapValues);

(Images can also be read in as .png files.) If this is a low-resolution raster then we may wish to smoothly interpolate it to higher resolution, done below with bicubic interpolation:

27 map.interpolate(factor=20, “cubic”);

Once we have this map, we can extract the values of the map at arbitrary locations—for instance, the elevations at which individuals live:

28 elevs = map.mapValue(inds.spatialPosition);

Further operations are available, including blurring and algebraic manipulations of the values. The dispersal method in Box 3 allowed individuals to move equally well in any direction. To guide movement with a map—for instance, to induce a preference for moving uphill, for this example where map values indicate elevation—we can use the map.sampleNearbyPoint() method:

29 inds = p1.individuals;30 pos = map.sampleNearbyPoint(inds.spatialPosition, INF, “n”, SM);31 inds.setSpatialPosition(pos);

This will move each individual to a new location sampled nearby and weighted by map value: if the original location is $$$ \boldsymbol{x} $$$, $$$ \rho $$$ is a Gaussian kernel (specified as type “n”) of width SD, as in $$$ {\sigma}_D $$$, and the value of the map at location $$$ y $$$ is $$$ m(y) $$$, then the new location is chosen with density proportional to $$$ m\left(\boldsymbol{y}\right)\rho \left(\left(\boldsymbol{y}-\boldsymbol{x}\right)/{\sigma}_D\right) $$$.

Most real populations are far from uniformly distributed in space, and in most cases, the underlying cause is thought to be environmental heterogeneity. Until this point, we have considered simulations of homogeneous landscapes only. Such “flat” landscapes are useful for developing models and/or theory, but incorporating aspects of real landscapes can make simulations more realistic. Spatial heterogeneity can be introduced simply by making some parameter of the model, such as fecundity or mortality ($$$ f $$$ or $$$ \mu $$$ as defined in Section 2), vary across space—in which case we can visualize the parameter as a map. Such maps might represent specific environmental conditions, habitat boundaries, or abstract habitat quality.

Raster-based images provide a convenient way to introduce a map of spatial heterogeneity into a simulation framework that cannot directly read geospatial data formats. A monochrome .png file consists of a rectangular grid of pixels with integer values between 0 and 255. These values can be shifted and scaled to lie within a useful range. See Case Studies (Section 9) for several examples of the use of images, and Box 5 for an example using .csv data. High-resolution maps to use as source images are publicly available from various sources, including NASA's Earthdata platform (2024), ESA's Earth System Data Lab (2024), or PRISM Climate Group (2024). Open-source tools for processing remote-sensing data are also available (see Montero et al. 2023).

To use a raster-based image, we must decide how the bounds of the image map onto geographic space. In general, we match the (rectangular) image with the (rectangular) spatial area to be simulated. The image is represented as discrete pixels, but we need to obtain values of the map at arbitrary locations (not just at the center each raster pixel). This can be done by either associating each pixel's value with the rectangle it would visually cover (as seen in Box 2), or associating each pixel's value with the corresponding cell's midpoint and interpolating between.

In the simple example shown in Figure 5, an image represents the altitude map of a mountain, with darker red shades indicating higher elevation. The local carrying capacity, $$$ K\left(\boldsymbol{x}\right) $$$, is modified by the value of the map at $$$ \boldsymbol{x} $$$ so that high-elevation locations can support more organisms. Comparing the provided altitude map (background of Figure 5A) with a summary of realized average density (Figure 5B), demonstrates that population density roughly matches the underlying elevation map.

[IMAGE OMITTED. SEE PDF]

A more precise comparison (Figure 5C) shows that $$$ n\left(\boldsymbol{x}\right) $$$, the local density at $$$ \boldsymbol{x} $$$ defined in Equation (1), is not exactly $$$ K\left(\boldsymbol{x}\right) $$$. In most situations, the equilibrium density is below $$$ K\left(\boldsymbol{x}\right) $$$. This is probably for two reasons: first, stochasticity usually reduces equilibrium density (see Section 7.3 and Appendix B.4); second, since the mountain is conical, most locations $$$ \boldsymbol{x} $$$ are surrounded by more low elevation area (where $$$ K $$$ is lower than $$$ K\left(\boldsymbol{x}\right) $$$) than higher area (where $$$ K $$$ is higher than $$$ K\left(\boldsymbol{x}\right) $$$), making the net flux of migrants at $$$ \mathbf{x} $$$ negative.

The degree of deviation of $$$ K\left(\boldsymbol{x}\right) $$$ from $$$ n\left(\boldsymbol{x}\right) $$$ can depend on other factors such as $$$ {\sigma}_D $$$. When dispersal is short ($$$ {\sigma}_D=0.25 $$$, black points), areas with carrying capacity below about 5 individuals per unit area are not self-sustaining. (See Appendix B.3 for more discussion.) On the other hand, with long-range dispersal ($$$ {\sigma}_D=4 $$$, green points), the overall relationship between density and carrying capacity is flattened as offspring from high-fecundity areas end up across the entire range. However, most offspring are still produced near the top of the mountain, and lower elevations are maintained by source-sink dynamics.

6

BOX

Adaptive empirical tuning for emergent parameters.

We often want a simulation to match a given estimated or observed density. However, this is not as simple as setting the value of $$$ K $$$ (local carrying capacity, introduced in Section 2) in the code, because population density is an emergent quantity—a complex consequence of births, deaths, movement, and local interactions. Fortunately, matching even emergent quantities is possible in SLiM. The code in this box dynamically adjusts parameters to make population density match a desired value. The same general technique may be applied to other quantities such as mean age or degree of clustering.

We first define a global parameter ADJ that will be adjusted:

32 defineGlobal(“ADJ”, 1.0);

Then we modify the density regulation code from Box 1 so that ADJ adjusts the carrying capacity:

33 inds = p1.individuals;34 density&#x02009;=&#x02009;11.localPopulationDensity(inds);35 u = density / ((1 + f) * K);36 inds.fitnessScaling&#x02009;=&#x02009;1 / (1 + ADJ * f * u);

In each time step, ADJ is updated by a factor $$$ \exp \left[\alpha \left(Y-K\right)\right] $$$, where $$$ Y $$$ is global population density (population size divided by total area) and $$$ \alpha $$$ is an update rate (e.g., $$$ \alpha =0.01 $$$):

37 obsDensity = p1.individualCount / (WIDTH * HEIGHT);38 defineGlobal(“ADJ”, ADJ * exp((obsDensity - K) * ALPHA));

When $$$ Y>K $$$, we want the density to be lowered in the next cycle, so ADJ is increased, decreasing carrying capacity. On the other hand, when $$$ Y<K $$$, ADJ is decreased so that we get a higher realized global population density $$$ Y $$$ in the next cycle.

Importantly, this tuning should stop after an initial “burn-in” period, and the appropriate value for ADJ, once found, should be hard-coded into the final model. If changes occur in the simulation that affect population size (e.g., a reduction in habitat), this adaptive code will force the density back to the chosen value, which is generally not desirable—the population size should, in fact, change in response to such changes in the simulated conditions.

Density Dependence, Life History, and Stochasticity

In our initial model, we used Beverton–Holt density-dependent feedback on mortality to control global population size (Section 2.2). Here we give examples that have the same Beverton–Holt form for the net effect of density on population regulation but differ in other ways.

Fecundity Regulation

The Beverton–Holt density-dependent regulation in Equation (4) has the probability of death increase with local density while fecundity stays constant. Alternatively, we can set the probability of death to a constant: $$$ \mu (u)={\mu}_0 $$$, and then solve for $$$ f(u) $$$ to obtain the Beverton–Holt form $$$ F(u)=\alpha \left(\left(1+a\right)/\left(1+ au\right)-1\right) $$$. (It turns out we will need the scaling factor $$$ 0<\alpha <1 $$$ to make a model with positive birth rates and death probabilities between 0 and 1.) Then, fecundity should depend on scaled density $$$ u $$$ as follows:5$$$ f(u)=\frac{1}{\left(1-{\mu}_0\right)}\left(\alpha \frac{\left(1+a\right)}{\left(1+ au\right)}+\left({\mu}_0-\alpha \right)\right). $$$

This functional form is shown in the middle panel of Figure 1.

Compensatory Regulation on Juvenile and Adult Mortality

To stabilize the population, the effect of density on net reproductive rate must be negative, but the effects on individual demographic components (such as birth or death) can be positive, as long as they are compensated for by other components. For instance, the probability of survival $$$ 1-\mu (u) $$$ can increase with density as long as fecundity $$$ f(u) $$$ decreases faster. One way to do this is to set6$$$ f(u)=\frac{f_0}{1+ bu}\kern1em \mathrm{and}\kern1em 1-\mu (u)=\frac{1+F(u)}{1+f(u)}. $$$where $$$ F(u)=\alpha \left(\left(1+a\right)/\left(1+ au\right)-1\right); $$$ that is, the Beverton–Holt form again. This produces valid survival probabilities if $$$ {f}_0\ge a\alpha $$$. Depending on the choice of the parameters, $$$ f(u) $$$ may increase or decrease with $$$ u $$$, and is not necessarily monotone. One choice of valid parameters is shown in the rightmost panel of Figure 1.

Our three examples (4), (5), and (6) (illustrated in Figure 1) might all be reasonably called “Beverton–Holt” models, although they differ substantially in the underlying mechanism of regulation. Although they have similar behavior around equilibrium density, they have quite different life-history implications. Most strikingly, in the first model, mortality increases with density; in the second, mortality is constant, while in the third, mortality decreases with density. There are corresponding differences in age structure among the models (as shown in Figure S5), although the dynamics of total population size are similar. More examples along these lines are given in Appendix C.

Stochastic Effects

Despite all this theory, in practice, equilibrium density is usually not $$$ K $$$—it is often lower. (See Figure 5 and Figure S6, for instance.) This is due to various forms of stochasticity. One is random lack of mates (when $$$ {N}_M $$$ is small); another is local extinction caused by temporal population fluctuations (when $$$ {N}_X $$$ is small). These effects can be quite troublesome when setting up computational experiments across a range of parameters, especially if we want constant total population size.

Apart from those issues, the most common reasons for a significant discrepancy between the realized density and the “desired” density (set by $$$ K $$$) have to do with when and how density is measured. First: when is density measured? In each time step, density increases after birth and decreases after death; $$$ K $$$ can match at most one of these times. Second: where is density measured? We naturally look at total population size divided by total area; however, on average individuals experience a higher density, since they themselves count toward their own local density. So, the correct comparison is of $$$ K $$$ to local population density averaged across the location of all individuals. Since this oversamples areas of higher density, this density will be higher than “total population size divided by area.”

All of these issues are discussed in much greater detail in Appendix B. For practical reasons, it often suffices to simply be aware that population density is fundamentally an emergent property, determined in complex ways by nearly all parts of the life cycle. If a precise total population size is desired in a particular simulation, a simple solution is to adjust some parameter (e.g., the birth rate) until the desired value is achieved (see Box 6), but of course, that will alter the dynamics of the simulation in other respects.

[IMAGE OMITTED. SEE PDF]

Natural Selection

It turns out that there are many different ways to implement natural selection in an individual-based simulation that all map to the same abstract models from population genetics, yet produce distinct outcomes. Modelers may choose to have an individual's genotype affect its choice of mate, its preference for certain habitats, its longevity, or any number of other life-history traits. These effects may differ across a landscape, leading to local adaptation. In this section, we demonstrate selection that acts on birth and death, building on Section 2. Just as we could regulate population size in our individual-based simulation through birth rate, death rate, or both, one may achieve the same increase in “fitness” by increasing fecundity or decreasing mortality, and we will see that the selected allele thus affects equilibrium population density. The resulting models might be classified as having “hard” selection (Wallace 1975); “soft” selection might act through increased or decreased chance of being chosen for mating (and would usually cause little to no change in equilibrium population density). For a view of the literature on these ecological-evolutionary models, see Mallet (2012) or Travis et al. (2023).

Much of population genetics selection theory is expressed in terms of the “selection coefficient” of a variant, usually denoted $$$ s $$$, which is usually defined in the context of a Wright–Fisher model as the change in relative probability it confers on the individual to be “selected” to provide offspring to the next generation. The point in the life cycle when the variant confers its selective effect is known to produce differences even in nonspatial models (Bodmer 1965; Nagylaki and Crow 1974), and so analogies between the dynamics of such alleles and those in a Wright–Fisher model are necessarily approximate. This is important when attempting to match simulation results to theory: just because a variable that affects survival or fecundity is named $$$ s $$$ does not mean that using its value in expressions derived from the Wright–Fisher model correctly predicts the probability of establishment, mean frequency, or other classical quantities, even in a nonspatial simulation.

Below, we demonstrate two models in which selection acts on either mortality or fecundity. We choose our parameter $$$ s $$$ so that a variant's frequency changes on average by a factor of $$$ 1+ hs $$$ per time step when rare, where $$$ h $$$ is the dominance coefficient. (Note that this may not be the best definition for use with classical formulas, which often measure change per generation.) For a historical and philosophical review of definitions and measures of selection see Endler (1986). Implementation in SLiM is described in Box 7.

Defining Allelic Growth Rate

Recall from Section 2 that $$$ f(u) $$$ and $$$ \mu (u) $$$ are, respectively, the per-capita mean number of offspring and probability of death per time step for an individual experiencing scaled density $$$ u $$$. For selection, these vital rates should depend on the individual's genotype. So, let $$$ {f}_k(u) $$$ be the mean fecundity of an individual with $$$ k $$$ copies of a focal allele, and similarly $$$ {\mu}_k(u) $$$ the probability of death. When the focal allele is rare, most of the population use $$$ {f}_0(u) $$$ and $$$ {\mu}_0(u) $$$, but a few individuals use $$$ {f}_1(u) $$$ and $$$ {\mu}_1(u) $$$. So, then the per-capita rate at which the number of focal alleles grows when at scaled density $$$ u $$$ in an outcrossing species is7$$$ G(u)=\left(1-{\mu}_1(u)\right)\left(1+\left({f}_0(u)+{f}_1(u)\right)/2\right)-1. $$$

This is almost the same as the expression for $$$ F(u) $$$, the expected change in the number of individuals from Equation (2), except that the fecundity term is $$$ \left({f}_0(u)+{f}_1(u)\right)/2 $$$; this is because the offspring of a heterozygous parent will only inherit the focal allele from that parent half the time.

Parameterizing Selection

First we consider mortality-based selection. Suppose an additive allele (with $$$ h=1/2 $$$) increases survival by $$$ s/2 $$$ per copy ($$$ k $$$), so that8$$$ 1-{\mu}_k(u)=\left(1-\mu (u)\right)\left(1+ ks/2\right). $$$

If the allele doesn't affect fecundity (meaning $$$ {f}_k(u)=f(u) $$$), then we plug Equation (8) into Equation (7) with $$$ k=1 $$$ to find that per-capita allelic growth rate in an outcrossing species is9$$$ G(u)=\left(1-\mu (u)\right)\left(1+f(u)\right)\left(1+s/2\right)-1. $$$

We have set things up here so that $$$ s $$$ means what we want: at the population's equilibrium density, the allele's growth rate when rare is $$$ G(1)=s/2 $$$.

Now suppose that the allele increases fecundity and does not affect survival (meaning $$$ {\mu}_k(u)=\mu (u) $$$). We define10$$$ {f}_k(u)=\left(1+f(u)\right)\left(1+ ks\right)-1. $$$

With this form for fecundity selection, the allelic growth rate when rare is the same as the rate we calculated for mortality-based selection (Figure 6). We can check this by plugging in our fecundity expression (10) into Equation (7) with $$$ k=1 $$$ and verifying that we get the same expression as Equation (9). Therefore, at equilibrium, the per-capita growth rate of the number of alleles is again $$$ s/2 $$$. Roughly speaking, each allele increases fecundity by $$$ s $$$, rather than the $$$ s/2 $$$ in Equation (8) for mortality selection, because here the effects of fecundity are only affecting half of the parents (the offspring-bearing ones).

Spatial Selection Mechanisms in Practice

Above we laid out deterministic, large-population-size theory that suggests that a variant that affects survivorship or fecundity might have similar frequency dynamics. However, in practice, there are some meaningful differences. Figure 7 shows selective sweeps in spatial and nonspatial simulations. All simulations use density-dependent feedback on mortality (via Equation (4)), but selection acts on either fecundity (via Equation (10)) or mortality (via Equation (8)). In nonspatial simulations, the global density (total number of individuals divided by total area) is substituted for local density when computing the probability of survival.

[IMAGE OMITTED. SEE PDF]

Some features of the sweep experiments are expected. Selective sweeps progress more slowly in spatially structured populations. The allele-frequency trajectories of de novo mutations under either mortality- or fecundity-based selection (Figure 7A), at first increase at the same speed in spatial and nonspatial sweeps, as predicted (i.e., $$$ G(u) $$$ is the same for both when the allele is rare). However, the selection mechanism introduces a few differences. As the beneficial allele fixes, equilibrium population density increases by a greater amount under fecundity-based selection (Figure 7B) than mortality-based selection. To understand why this is, consider the population after the fixation of the beneficial allele: all individual birth and death rates are now $$$ {f}_2(u) $$$ and $$$ {\mu}_2(u) $$$. We are at a new equilibrium density $$$ {u}_{\ast } $$$, which occurs when net per-capita change in population size is zero. Solving Equation (2) with our post-fixation vital rates $$$ \left(F\left({u}_{\ast}\right)=0\right) $$$, the new equilibrium density $$$ {u}_{\ast } $$$ is $$$ {u}_{\ast }=1+s\left(1+\frac{1}{f_0}\right) $$$ for mortality-based selection, and $$$ {u}_{\ast }=1+2s\left(1+\frac{1}{f_0}\right) $$$ for fecundity-based selection. In other words, the increase in equilibrium density for fecundity-based selection should be roughly twice that of the increase in mortality-based selection.

Fecundity-based selection also causes the mean individual age to drop as the beneficial allele increases in frequency (Figure 7C). Population size increases as the fecundity-based sweep progresses, but because the beneficial allele does not confer protection against mortality, individuals are subject to the negative effects of increased population density.

Further differences between the spatial and nonspatial simulations are likely due to the spatial density regulation processes that cause realized population sizes in spatial simulations to be smaller than in their nonspatial counterparts. (See Section 7.3 and Appendix B.4 for discussion).

7

BOX

Selection.

In nonWF SLiM models, the “fitness” property of an individual is the probability that the individual survives to the next time step. To make each mutation affect survival by a factor of S_MORT, we simply declare:

39 initializeMutationType(“m1”, 0.5, “f”, S_MORT);

(The factor of 0.5 is the dominance coefficient: heterozygotes will have fitness multiplied by 1 + S_MORT*0.5, and homozygotes by 1 + S_MORT).

To have the same type of mutations also affect fecundity with selection coefficient S_FEC, we need to account for genotype when setting up offspring numbers (as in Box 1):

40 indiv_s&#x02009;=&#x02009;2*S_FEC * inds.countOfMutationsOfType(m1) / 2;41 inds.tag = rpois(length(inds), mean=FECUN * (1 + indiv_s) + indiv_s);

As written, the effects of any mutation of type m1 on fecundity are additive across loci and copies; other arrangements are possible.

Case Studies

How can we combine all the strategies we have seen thus far into a useful model of a living system? Here we illustrate how the spatial modeling framework we have introduced can be used to model complex scenarios such as (1) temporal change, (2) complex life cycles, (3) continental-scale systems, and (4) competition for resources.

We motivate each scenario with an organism. The scenarios are not meant to be complete or accurate models of the population biology of these organisms; rather, they illustrate how to apply the concepts presented in this paper to a study system. Additional methodological information for each vignette can be found in Appendix E and our code example repository .

Changes Over Time: Pikas and Environmental Change

Modeling spatial heterogeneity (as introduced in Section 6) makes simulations more realistic and informative, but sometimes temporal change in landscape variables is just as important. Here we explore how both seasonal fluctuations and globally rising temperatures affect the population dynamics of an alpine organism: the pika.

Pikas (Ochotona daurica) are adapted to mountainous habitat at relatively high elevation and cannot survive the extreme heat (or cold) at lower (or higher) elevations (Beever et al. 2010). Such ecology makes pikas potentially vulnerable as global temperatures are expected to increase over time.

Modeling Approach

We incorporate three types of temporal change in temperature: (1) within-year seasonal change, (2) random, autocorrelated fluctuations between years, and (3) steady global temperature rise. Because temperature varies with elevation and elevation varies dramatically in mountainous regions, we model temperature as a function of elevation. We used a topographic map for a region of Rocky Mountain National Park in Colorado to calculate temperatures for each point in space and time during our simulation.

We connect temperature to fitness by killing offspring that are exposed to temperatures beyond the minimum ($$$ -{5}^{\circ}\mathrm{C} $$$) or maximum ($$$ {28}^{\circ}\mathrm{C} $$$) sustained by pikas (Beever et al. 2010). Each tick of our simulation represents 1 year. To account for within-year seasonal variation (i.e., winter cold and summer heat), we narrow the viability range by the yearly variation in temperature ($$$ {s}^{\mathrm{seas}} $$$), defining the probability of survival of a pika at location $$$ \mathbf{x} $$$ in year $$$ t $$$ to be$$$ 1-\mu \left(\boldsymbol{x},t\right)=\left\{\begin{array}{ll}0,& \mathrm{if}\ T\left(\boldsymbol{x},t\right)<\left(-5+{s}^{\mathrm{seas}}/2\right)\\ {}0,& \mathrm{if}\ T\left(\boldsymbol{x},t\right)>\left(28-{s}^{\mathrm{seas}}/2\right)\\ {}\frac{1}{1+ fu\left(\boldsymbol{x}\right)},& \mathrm{otherwise}\end{array}\right. $$$where the last term is our usual Beverton–Holt mortality regulation with local scaled density $$$ u\left(\boldsymbol{x}\right)=n\left(\boldsymbol{x}\right)/K $$$ (see Equation (4)). In this equation, $$$ T\left(\boldsymbol{x},t\right) $$$ represents the midpoint of seasonal temperatures in year $$$ t $$$, defined below, and the within-year variability $$$ {s}^{\mathrm{seas}} $$$ broadens the temperatures experienced by the pikas around that midpoint. For example, if $$$ T\left(\boldsymbol{x},t\right) $$$ were $$$ -{4}^{\circ}\mathrm{C} $$$ (inside the viability range), but seasonal temperatures varied by $$$ {4}^{\circ}\mathrm{C} $$$, the winter would be $$$ -{6}^{\circ}\mathrm{C} $$$—cold enough to kill.

Between-year temperature changes are modeled by the temperature function$$$ T\left(\boldsymbol{x},t\right)={T}^{\mathrm{elev}}\left(\boldsymbol{x}\right)+{T}^{\mathrm{fluct}}(t)+0.016t $$$which combines elevation-related temperature change $$$ {T}^{\mathrm{elev}} $$$, correlated year-to-year fluctuations $$$ {T}^{\mathrm{fluct}} $$$, and global directional climate change.

$$$ {T}^{\mathrm{elev}}(x) $$$ is calculated from the elevation map following fitted models from Antonio-Juan Collados et al. (2020). To create autocorrelated year-to-year fluctuation, we let $$$ {T}^{\mathrm{fluct}}(0)=0 $$$ and for $$$ t\ge 1 $$$,$$$ {T}^{\mathrm{fluct}}(t)={p}^{\mathrm{fluct}}T\left(t-1\right)+\upvarepsilon (t) $$$where $$$ \upvarepsilon (t) $$$ is Normally distributed noise with mean zero and standard deviation ($$$ {s}^{\mathrm{fluct}} $$$). The parameters $$$ {p}^{\mathrm{fluct}} $$$ and $$$ {s}^{\mathrm{fluct}} $$$ are known as persistence and shock, respectively, and determine how correlated the noise is between years. Finally, we increase the global temperature over time by adding 0.016 C per year (Foster and Rahmstorf 2011).

Observations and Extensions

The resulting simulation is a habitat suitability model in which the population's geographic distribution moves toward higher elevation as the global temperature increases (Figure 8).

[IMAGE OMITTED. SEE PDF]

We found that if the magnitude of random variation around the expected annual temperature is large ($$$ {s}^{\mathrm{fluct}}={5}^{\circ } $$$C) the probability of extinction increases significantly, particularly for early generations. In other words, one or a few bad years was devastating for pika populations. This result may provide a useful lesson: even if a species appears to be thriving, the long-term success of the species is not guaranteed.

Surprisingly, simulated population size increased in the intermediate-term future, since for the first few hundred years, habitat losses at lower elevation were more than compensated for by habitat gains at higher elevation. This illustrates that results from simulation on a specific map may not be generalizable.

Here, we use an empirically informed linear trend with Gaussian noise as our global temperature change schedule. An alternative strategy for simulating temperature changes would be to pre-process multiple temperature maps reflecting different years and continually feed the simulation new maps over time. Such maps could be generated directly from a climate model rather than using a direct function of elevation.

Life Cycle Stages: Mosquitos in Burkina Faso

So far, our models have not taken life cycle stages into account; individuals have been able to mate immediately after they are born, and their survival has not been age-dependent. However, for many organisms, their patterns of mobility, competition for resources, and mating capability are age-dependent. Here, we demonstrate how to simulate a population with distinct juvenile and adult phases. Specifically, we set up a simulation of mosquitos in Burkina Faso in West Africa, inspired by North et al. (2019).

Modeling Approach

In this model, individuals begin life as juveniles, and mature into adults after a fixed time span. Population regulation and individual behavior depends on the life stage. For adults this is a constant survival probability of $$$ 0.875 $$$ per day (time step), independent of the landscape map. The population is regulated by density-dependent survival of the larvae, which varies across the map.

Following the model in North et al. (2019), we set the carrying capacity of larvae based on water availability. Larval carrying capacity is only non-zero at the outlines of water features extracted from GIS data of inland water in Burkina Faso. There, the carrying capacity fluctuates seasonally to mimic rainy and dry seasons:$$$ K\left(\boldsymbol{x},t\right)={K}_{\mathrm{base}}\left(\boldsymbol{x}\right)+{K}_{\mathrm{rain}}(t) $$$where $$$ {K}_{\mathrm{base}}\left(\boldsymbol{x}\right) $$$ is obtained from the map of inland water, and $$$ {K}_{\mathrm{rain}}(t) $$$ is a sinusoidal function with a period of 365 days and a minimum value of $$$ 0 $$$.

We then use a Beverton–Holt form for the survival probability of larvae. Since there are ten age classes for larvae (and survival of each depends on the total number across all ages), parameterizing the model so that local density is (roughly) $$$ K $$$ involves solving a system of equations (in a matrix model) described in Appendix E.2.

A female adult mosquito mates with an adult male within the maximum mating distance and lays eggs by sampling a location within her dispersal radius weighted by carrying capacity. Each day, adult mosquitos move by a random displacement sampled from a Gaussian distribution, whereas larvae do not move from their original location until they mature.

Observations and Extensions

This model simulates a mosquito population with a structured life history. The population size of larvae and adults fluctuates periodically, following precipitation levels with a slight lag (Figure 9). Recall that adult mortality is simply constant, so the periodic fluctuations of adults are mediated through larval carrying capacity dynamics.

[IMAGE OMITTED. SEE PDF]

There are a few immediate extensions to our simulation. Here, we did not distinguish between perennial and fluctuating water sources. Simulating water sources that appear and disappear (adding a $$$ t $$$ argument to $$$ {K}_{\mathrm{base}}\left(\boldsymbol{x}\right) $$$) could result in extinction dynamics such as those observed in the pika example above. Similarly, the amount and duration of rainfall could be location-dependent (adding an $$$ \boldsymbol{x} $$$ argument to $$$ {K}_{\mathrm{rain}}(t) $$$).

We set the maximum of the $$$ {K}_{\mathrm{base}} $$$ and $$$ {K}_{\mathrm{rain}} $$$ functions to $$$ 0.002 $$$ individuals per square meter, which is probably much lower than in nature, to keep the memory usage of our simulation low. To run a simulation with a realistic (huge) number of mosquitoes we would need to make some efficiency improvements. Fortunately, an example of the necessary techniques follows.

Continental-Scale Systems: Invasion of the Cane Toads

There is a perception that individual-based spatial models are so slow that they are meaningfully limited in the population size that can be modeled. Our goal with this vignette is to demonstrate that this is not necessarily the case by showing how to model large-scale populations and landscapes with relative ease and efficiency.

Cane toads (Bufo marinus or Rhinella marina) are native to Central and South America and were intentionally introduced to the Northeast coast of Australia in 1935 as a pest-control measure. Since their introduction, cane toads have experienced explosive population growth, with hundreds of millions of individuals spreading over several million square kilometers, resulting in considerable negative economic and ecological impacts (Shine 2010; Urban et al. 2008).

Modeling Approach

We use information about the biology of cane toads, when available, to parameterize the model. Toad population densities in established populations have been estimated around 8000 per square kilometer (Freeland 1986). Based on this, and then tailored heuristically over several trial runs to be computationally feasible and to produce a reasonably realistic pattern of range expansion, we settled upon a local carrying capacity of 1000 per square kilometer. Telemetry data shows cane toads can travel up to 0.2 km per day (Shine et al. 2021). In order to interpret simulation time steps roughly as years, we set the spatial scale parameters $$$ {\sigma}_D,{\sigma}_X $$$, and $$$ {\sigma}_M $$$ to 20 km. For dispersal, we use a Student's $$$ t $$$-distribution, which provides more long-range dispersal events than a Gaussian kernel (Figure 3C).

We model survival probability as a function of precipitation due to its similar distribution to empirical occurrence data, which seems reasonable given cane toads' known sensitivity to moisture conditions (Child et al. 2009; Cohen and Alford 1996). Specifically, we multiplied our typical density-dependent survival probability $$$ 1-\mu (u) $$$ by $$$ 1/\left(1+\exp \left(-\left(\alpha +\beta P\left(\boldsymbol{x}\right)\right)\right)\right) $$$, where $$$ \alpha $$$ and $$$ \beta $$$ control the intercept and slope of the survival curve, and $$$ P\left(\boldsymbol{x}\right) $$$ is the amount of yearly precipitation at location $$$ \boldsymbol{x} $$$.

The invasion began with 10,000 individuals (though the number of individuals actually released in 1935 was likely greater (Shine et al. 2020)) at locations randomly sampled from the first 4 years of the observed occurrence data.

Indeed, mating and interaction neighborhoods were large, which initially led to prohibitively long runtimes. To make the simulation more efficient, we followed the approach described in Box 8, modified to measure the local density of individuals per unit of habitable area. Mate choice was also modified as described in Box 9.

Observations and Extensions

We were able to approach the true spatial and population scale of the cane-toad invasion, with a final census size of about 120 million individuals, nearing the estimated modern census size of Australian cane toads (WWF Australia 2019), after running for about 5 days and using 200 GB of RAM at maximum. We visually compared empirical occurrence data for cane toads to our simulations, with and without annual precipitation's effect on survival (Figure 10). While there are obvious differences in densities and locations between the simulated and observed data, it is clear that modeling annual precipitation's effect on survival greatly improves the likeness.

[IMAGE OMITTED. SEE PDF]

The approach described here is similar to classical niche modeling (Peterson 2001), which has been used extensively to model cane-toad distributions (Shine 2010), but is somewhat simplified by only using one environmental variable. This would be a straightforward extension: multiple environmental variables could be combined into a single habitability map for use in the model, with no impact on runtime performance. The benefit of this vignette's approach over previous approaches is the way that it combines information from the explicit individual-based simulation and environmental data, producing a simulation at a realistic scale with respect to both landscape size and population size.

8

BOX

Using maps for faster spatial interactions.

In Box 1, we estimated the local population density for each individual with Equation (1), using the localPopulationDensity() method. For each individual, this method sums the “interaction strength” (i.e., kernel density) for every other individual within the provided maximum distance. So, if the total number of individuals is $$$ T $$$ and the typical number of neighboring individuals is $$$ {N}_X $$$ (defined in Section 3), then the complexity of this operation is $$$ {TN}_X $$$. If each individual has a large number of neighbors, this can be quite costly. However, if the number of neighbors is large, it should work just as well to (a) create a (discretized) map of local density, and (b) look up the value of the local density experienced by each individual on that map. Map lookup is quick, so if the cost of creating the map is smaller than $$$ {TN}_X $$$, we will have a more efficient model.

We can create a map whose value at $$$ \boldsymbol{x} $$$ is (approximately) given by Equation (1) in two steps: (1) use summarizeIndividuals() to measure the number of individuals per unit area in each cell in a grid (see Box 2), and (2) smooth() this map using the appropriate kernel, so that the density value at a given point in the resulting map depends upon an appropriate weighted average of the individuals per unit area across nearby cells of the grid. Using the same kernel as in Box 1:

42 grid_dims = ceil(2 * p1.spatialBounds / SX);43 raw = summarizeIndividuals(p1.individuals, grid_dims, p1.spatialBounds,44 operation=”individuals.size();”, perUnitArea=T);45 density_map = p1.defineSpatialMap(“density”, “xy”, raw);46 density_map.smooth(SX * 3, “n”, SX);47 defineGlobal(“DENSITY”, density_map);

Then, we can modify the code of Box 1 to use the map instead:

48 inds = p1.individuals;49 density = p1.spatialMapValue(DENSITY, inds.spatialPosition);50 u = density / ((1 + f) * K);51 inds.fitnessScaling&#x02009;=&#x02009;1 / (1 + f * u);

This obtains (nearly) the same value as would localPopulationDensity() if the resolution of the map should be finer than the scale over which the density kernel changes. In this example, that scale is SX, as in $$$ {\sigma}_X $$$, so we have ensured that the map has cells of size smaller than SX/2. The approximation is examined in Appendix F.

9

BOX

Using maps for faster mate choice.

A similar problem as in Box 8 arises when choosing mates: even though only one mate needs to be chosen, the underlying operation is of order $$$ {TN}_M $$$, where $$$ T $$$ is the global population size. The same map of density can be used to solve this problem as well: instead of choosing an individual with probability proportional to a kernel, it is (nearly) equivalent to: (1) choose a point in space nearby, with probability proportional to the map multiplied by the kernel, and then (2) take the individual closest to that point. Recall the number of possible mates scales as $$$ {N}_M=4{\sigma}_M^2K $$$ from Section 3. This implies that the number of neighbors grows linearly with $$$ K $$$. Thus, we can keep the number of potential mate roughly constant regardless of how large local density is by rescaling the maximum distance in the InteractionType used for mate choice by $$$ 1/\sqrt{K} $$$ from the code in Box 4:

52 initializeInteractionType(2, “xy”, maxDistance=5/sqrt(K), sexSegregation = “FM”);53 i2.setInteractionFunction(“f”, 1.0);

We have set the maximum distance in the interaction kernel to be a value that should give us around 25 neighbors for each individual; however, if density varies significantly across the landscape, this may make some individuals in low density areas fail to mate. Then, we choose the mate as follows:

54 mate_location = DENSITY.sampleNearbyPoint(individual.spatialPosition,55 3*SM, “n”, SM);56 mate = i2.nearestNeighborsOfPoint(mate_location, p1, 1);

Here the specification of a Gaussian mate choice kernel with standard deviation SM has moved from the definition of the InteractionType to the sampleNearbyPoint call: given a location $$$ x $$$, a map with value $$$ m(y) $$$ at $$$ y $$$, and a kernel $$$ \rho $$$, this returns a random point $$$ z $$$ sampled with probability proportional to $$$ m(z)\rho (z) $$$. We then choose the mate as the individual nearest to that point. The approximation is examined in Appendix F.

Resource-Explicit Competition: Monarchs and Milkweed

In the preceding examples, we have regulated populations through competitive interactions between individuals: either explicitly, as in Box 1, or in a space-averaged manner, as in Box 8. Population regulation in this model is managed quite differently. This is a “resource-explicit” model, wherein the population is extrinsically regulated by the availability of an external resource (Champer et al. 2024), as outlined in Box 10.

We simulate monarch butterflies (Danaus plexippus) and the milkweed (Asclepias spp.) plants on which they lay their eggs. Though adult monarchs feed on nectar from numerous types of plants, monarch caterpillars are specialists that eat only milkweed (Oberhauser et al. 2004).

10

BOX

Resource-explicit foraging and mortality.

In this box, we outline the “resource-explicit” modeling approach described more fully in Champer et al. (2024). This approach implements density-dependent population regulation that is mediated indirectly through the availability of a resource, rather than directly through competitive interactions between individuals. A simple resource-explicit simulation contains two species in a multispecies SLiM model: the focal species, in subpopulation p1, and a species representing the resource, in subpopulation p2, in which individuals represent patches of the modeled resource.

During the foraging phase, individuals collect resources from nearby patches. Each patch can support a certain number of individuals per time step. If a patch can support 10 individuals, but there are 100 individuals nearby at a given time step, each individual will receive 10% of the amount of resource necessary to guarantee survival (but might also forage from other nearby patches). The amount of resource that each individual collects during the foraging phase is stored in the tagF property.

1 i1.evaluate(c(p1, p2));2 for (patch in p2.individuals) {3 customers = i1.nearestNeighbors(patch, p1.individualCount, p1);4 customers.tagF = customers.tagF + inds_fed_per_patch / size(customers);5 }

After foraging, individuals die if they have not consumed enough resources. In the code example below, this occurs when individuals reach age 2.

57 at_risk_indivs = p1.subsetIndividuals(minAge=2);58 mortality_indices = runif(size(at_risk_indivs)) > at_risk_indivs.tagF;59 dead = at_risk_indivs[mortality_indices];60 species.killIndividuals(dead);

If an individual collected at least one unit of the resource, its survival to the next time step is guaranteed. Otherwise, it survives with a probability equal to the total amount it collected. See Appendix F for discussion of the relationship to the method of Box 8.

Modeling Approach

Monarchs in the model progress through three life-cycle stages: caterpillar, pupa, and butterfly. Mortality is modeled differently in each stage. The food resource, milkweed, is directly included in the model as a second species.

During the first 2 weeks of their lives (time steps 0 and 1), caterpillars interact with nearby milkweed plants and accumulate resources. The amount of resources collected from a milkweed plant each week is inversely proportional to the number of competing caterpillars within a given interaction scale ($$$ {\sigma}_X $$$) of that plant. This means that a plant fed upon by more caterpillars will be depleted more quickly. When individuals reach their third week, they enter the pupa phase. At this time, survival is proportional to the amount of milkweed eaten as a caterpillar. (The survival probability is calculated only once when it becomes a pupa.) Surviving pupae become butterflies during their fifth week. Butterflies disperse, reproduce, and experience mortality at an age-dependent rate. Mortality is not density-dependent during the butterfly stage. Figure 11 shows a snapshot of the simulation, highlighting the dispersal and reproduction of individuals across the landscape.

[IMAGE OMITTED. SEE PDF]

Reproduction also involves a resource-explicit interaction, since monarchs only lay their eggs on milkweed. To accomplish this, a spatial interaction is evaluated between adult males and milkweed (with a much longer range than the interaction between caterpillars and milkweed); each milkweed caches a list of nearby males from that interaction. Next, a similar spatial interaction is evaluated between adult females and milkweed; females are iterated through, with each female selecting a mate from the males cached at the milkweed plants within the female's interaction range. Finally, the females randomly distribute their eggs at the plants where matings occurred.

Observations and Extensions

In addition to reflecting the life history of monarch butterflies, the resource-explicit modeling approach is highly performant. Each female monarch can lay several hundred eggs, very few of which go on to reach adulthood. As a result, the number of caterpillars in the model far exceeds the number of milkweed plants. Thus, regulating the size of the caterpillar population by evaluating a spatial interaction between caterpillars and plants is far more efficient than regulating the population by evaluating an interaction directly between caterpillars.

This approach also provides an intuitive way to investigate questions related to resource availability. For example, this model could probe the effects that milkweed habitat loss, caused by urbanization or climate change, could have on the viability of the local monarch population.

Discussion

The first spatial models intended to represent landscapes (as opposed to, say, two-deme island-mainland models) were based on partial differential equations and so did not explicitly represent individuals (Skellam 1951; Beverton and Holt 1957; Cantrell and Cosner 2004). Many others have used arrays of discrete populations (e.g., metapopulation models, Hanski 1997), or had individuals living on a regular grid (Epperson 2003). Although there are use cases for these, the main reason that individual-based continuous-space models are not more common may simply be convenience, both for mathematical analysis and software programming. The continuous-space formulation we use was introduced by Pacala (1986) and Bolker and Pacala (1997), and has been used in many theoretical studies (Dieckmann and Law 1999; Snyder and Chesson 2004; Etheridge et al. 2024).

A more modern extension to matrix population models, integral projection models, more commonly incorporate density dependence (Ellner et al. 2016), and may even estimate functional forms—see for example Adler et al. (2010). Much of the remaining literature on density dependence only considers the form of $$$ F(u) $$$, rather than separating out the effects of density on different life-history components (reviewed in Caswell 2000; Eskola and Geritz 2007). There are a wide variety of methods to infer $$$ F(u) $$$ from time series data (Lande et al. 2003), but statistical issues make the problem difficult (Freckleton et al. 2006). However, the “$$$ F(u) $$$” thus inferred is not necessarily the same as ours—many of these methods assume no demographic stochasticity, and hence an infinite population size. In other words, the $$$ F(u) $$$ that these infer is a landscape-scale relationship, averaging the net effects of birth and death over thousands or millions of individuals. Our $$$ F(u) $$$, on the other hand, is local, and describes how a single individual is affected by having a few more or less neighbors. These are related, but need not have the same functional forms. In particular, the behavior of $$$ F(u) $$$ for very large or small values of $$$ u $$$ is usually much more important when applied to an individual, because the local density around an individual might be proportionally much larger or smaller than average, due to random fluctuations, than the total population size across a landscape. These can have a strong effect—for instance, increasing population density can have a positive effect on population growth (for at least low enough densities), a dynamic which can have interesting and important effects (Courchamp et al. 2008).

Furthermore, the effects of density dependence are often mediated (at least in part) through interactions with other species, and the nature of these interactions may depend on environmental context. These interactions are often modeled as the net result of pairwise interactions, and various methods are used to estimate these potentially numerous and environmentally dependent effects (for recent examples, see Weiss-Lehman et al. 2022; Malyon et al. 2023). It is beyond the scope of this article to review the full range of approaches and possibilities—but note that there is no obstacle to simulation of multiple species whose interactions vary across space and/or time in SLiM (Haller and Messer 2023).

Discretized Continuous Space

There is not a hard distinction between “continuous” and “discrete” spatial models. However, in practice, discretized models usually only allow density-dependent regulation to act within demes: in other words, vital rates of an individual can only depend on the other individuals in the same deme. This is important because the degree of demographic stochasticity depends on deme size; therefore, different discretizations may produce different amounts of stochasticity. What demographic stochasticity there is shared among all individuals in the same deme, reducing the overall variance in the simulation. The effects of these discretization artifacts are not well described (but see Battey et al. 2020). If demographic rates within each deme instead depended on nearby demes, then this would not be an issue. However, we are not aware of any simulators that make this choice.

Empirical Estimation of Density Dependence

Huge amounts of observational effort have gone into detailed estimates of the demography of particular species, especially those of conservation or management concern. However, such observations by necessity describe a snapshot of demographic rates in particular conditions (or, averaged over a particular range of conditions). Estimating the functional responses to density necessary to describe long-term dynamics is much more difficult, although methodological progress has been made, for instance, by incorporating many sources of information with “integrated population models” (Zipkin et al. 2023). A great deal of ecological work has also sought to quantify the effects of density-dependent demographic feedback. In a spatial model, this population density is usually measured locally—around the individual in question—yielding concepts such as the “competition kernel” (Bolker and Pacala 1997) or “crowding index” (Pacala and Silander 1985). Density-dependent feedback has a key place in the theory of coexistence between species: the Janzen-Connell hypothesis suggests that stronger intra- than interspecific negative density-dependent effects could provide a mechanism for species coexistence (Janzen 1970; Connell 1971; Terborgh 2012; Hülsmann et al. 2021). A substantial number of studies across a variety of organisms (mostly plants) have quantified these effects, both within species (e.g., Weiner 1982; Silander and Pacala 1985; Specht and Arnold 2018; Spotswood et al. 2017) and between (e.g., Mack and Harper 1977; Song et al. 2021; Zaiats et al. 2021). Estimation of density-dependent interaction strengths between many species in a community is challenging, but important for understanding community assembly (e.g., Weiss-Lehman et al. 2022; Malyon et al. 2023). Adding to this complexity, the effects likely often depend strongly on age (Richardson et al. 2024).

Researchers in plant ecology have made the most progress toward empirical understanding of the mechanistic underpinnings of the sort of local density dependence we require. In practice there are a great many possible ways to quantify the cumulative effect of the neighbors of a single individual (Weigelt and Jolliffe 2003). Our formulation here treats all individuals equivalently, but in practice, one could include the effects of age or size. For instance, a common method in forestry modeling defines a “neighborhood competition index” for a given tree as the sum over all neighbors of their diameter (to some power) divided by distance (to another power) (Bella 1971; Daniels 1976; Canham et al. 2004). Empirical studies have estimated these kernels in a variety of situations: for instance, Teller et al. (2016) and Adler et al. (2018) use spline methods to flexibly estimate the effects of total area of nearby plants on a target plant's growth and survival.

Empirical Estimation of Dispersal

The concept of a dispersal distribution (as discussed in Section 4) is perhaps most well defined for plants, which (mostly) have only one opportunity to move during their lifetime, as a seed or other propagule. Seeds and pollen can be moved by gravity, wind (Nurminiemi et al. 1998), water (Murray 2012), or animals (Morales and Morán López 2022; Pons and Pausas 2007), and complex models for these have been developed and estimated from empirical data (Neubert et al. 1995; Tufto et al. 1997; Austerlitz et al. 2004; Katul et al. 2005). Of course, many animals have relatively stable locations as well, perhaps depending on the season: for instance, Paradis et al. (1998) review estimates of post-natal and breeding dispersal in many bird species. The combination of these various processes (which often happen across many different spatial scales) has been referred to as the “total dispersal kernel” (Rogers et al. 2019), and can easily lead to “long-tailed” distributions (Cain et al. 2000; Edwards et al. 2007) in which rare, long-distance events are important. Dispersal is notoriously hard to estimate, as it often requires observing rare events, but has been done in a variety of organisms including kangaroo rats (French et al. 1968), mosquitoes (Estep et al. 2010), Drosophila (Dobzhansky and Wright 1943), Prunus shrubs (Robledo-Arnuncio and García 2007), pine pollen (Robledo-Arnuncio and Gil 2005), and butterflies (Suchan et al. 2024). Dispersal can of course also depend on density (Harman et al. 2020) (e.g., if individuals preferentially disperse out of crowded areas). Clobert et al. (2012) review many aspects of dispersal, from parameterization and estimation to implications for ecology and evolution. See also Saastamoinen et al. (2018); Edelaar and Bolnick (2012).

Conclusion

Individual-based simulations are a powerful method for studying how demographic and population-genetic processes operate over continuous geographic space. Modelers must design rules for how individuals in the simulation interact with others nearby and how forces such as selection operate. Individual-based simulations are well suited to this problem since they are very flexible and can be tailored to a specific research system. However, flexibility can be both a blessing and a curse: it is easy to design a simulation with unstable population dynamics or unrealistic life-history traits. Similarly, stochasticity can cause a population that is intended to equilibrate to instead die out. Such problems likely reveal a flaw in our understanding of the system being modeled.

Here, we provide guidance and connections to the ecological literature for researchers interested in designing stable, efficient, and interpretable spatial simulations. Interpretability is a key advantage of spatial, individual-based simulations, since it can take substantial effort to translate the results of more abstracted models back to the domain of interest.

Realism is not a goal of our case studies, but each illustrates the degree of realism that can be obtained from SLiM without serious effort. Increasing computational efficiency and flexibility of simulation engines are bringing individual-based simulations closer to realistic models of the ecology of specific systems. Careful implementation of ecologically realistic evolutionary models will be important to many applied fields, such as understanding and predicting how climate change affects organisms' ranges, predicting the consequences of a gene-drive release in the wild, and rescuing species close to extinction. As the scale and specificity of in silico models improve, individual-based simulations will become valuable tools beyond academic research for management professionals in conservation, management, and public health.

Modeling Density Dependence

Setting out to write this paper, we hoped to provide a comprehensive yet simple guide to best practices in implementing density-dependent population regulation. Although we have provided one or two paths forward and elucidated many of the issues (see in particular Appendix B), careful empirical practitioners will soon encounter additional questions. What are some flexible and robust families of functional forms, and what aspects of these matter in practice? How can these be parameterized so that parameters naturally correspond to observable/interpretable quantities? How can these be fit to data? How should local habitat quality and density interact? At first, we imagined that the answers would be found in familiar names, and so the relationship between local density and fitness would be described as logistic, Beverton–Holt, Ricker, et cetera. However, we quickly found that these models were developed to describe population-level net changes, and so they not only do not account for individual-level stochasticity, but furthermore do not separate birth from death. This is an area of active work—see, for instance, Aoyama et al. (2022) and Adler et al. (2018) for recent good examples. Full exploration of these questions was too much for this paper.

Spatial Data and Niche Modeling

Several of our case studies use environmental variables to specify where on the simulated landscape organisms are most likely to live. The explosion in remote-sensing data provides many potential data sources for modeling spatial heterogeneity. However, what is often needed in a model is a composite proxy for “suitability” that can be incorporated into local demographics. The process of predicting where a species might or does live is known as Ecological (or, Environmental) Niche Modeling (Booth et al. 1988; Peterson 2001). This can be done in a variety of ways; for instance, one might model either the potential or the realized niche and predict probabilities of occurrence or population densities (reviewed in Sillero 2011). Environmental niche models are often used to predict suitable habitat either in other locations or other time periods (Werkowska et al. 2017; Yates et al. 2018), but resulting estimates can vary widely in quality and there are a number of statistical pitfalls (Sillero and Barbosa 2021) that simulation testing could diagnose and simulation-based inference could potentially help avoid.

While ecological niche modeling uses observational data to predict where organisms might live, a collection of landscape genetics tools tries to use genetic relatedness to predict where organisms move. The methods generate a map of “landscape resistance” that aims to describe how easily individuals move over different parts of the landscape (McRae et al. 2008). However, resistance models rely on correlating genetic distance to an abstract notion motivated by electrical circuits (reviewed in Peterson et al. 2019; Cruzan and Hendrickson 2020). They lack an underlying mechanistic model, so estimates can be problematic in practice (Cushman et al. 2013; Graves et al. 2013), and are expected to mislead in plausible situations such as biased dispersal (Lundgren and Ralph 2019). Again, simulation-based inference provides a promising route forward (Smith et al. 2023, 2024), since it does not rely on explicit likelihoods or other mathematical descriptions.

Sampling

Any simulation study that wishes to make comparisons to real data needs to also consider the sampling effort that led to those data, and realistic simulation of many sampling schemes can be daunting. In practice, sampling can strongly affect results (for an example in population genetics, see Battey et al. 2020). However, it is relatively easy to assess the robustness of results to variations in the sampling scheme. Furthermore, it is often possible to “oversample” simulations. For instance, in Smith et al. (2023), many simulated datasets can be obtained from each costly spatial simulation, simply by repeating the sampling effort (however, if sufficient simulations are not done, model performance will be poor). It would be useful to develop a standard set of tools that implements various sampling schemes for simulated spatial populations.

The Future

Although the spatial simulations we present here incorporate many more aspects of real organisms' lives than does the Wright–Fisher model, there are many things that we have not tried to explicitly model, such as seasonal migration, herding or flocking, territoriality, foraging strategies, microhabitat variation, broadcast spawning, resource storing, pollination, predation, facilitation, and other inter-species interactions. Any of these can be modeled in SLiM with more or less effort, and indeed many are described in the SLiM manual (Haller and Messer 2024). The decision of which aspects of biology to model is in practice made by prioritizing those aspects expected to substantially affect the question being studied. We are excited to see the wide variety of simulations that researchers develop in the future, as we explore these questions and build on each other's work.

Author Contributions

Elizabeth T. Chevy: conceptualization (lead), investigation (lead), methodology (lead), visualization (lead), writing – original draft (lead), writing – review and editing (lead). Jiseon Min: conceptualization (lead), formal analysis (lead), investigation (lead), methodology (lead), visualization (lead), writing – original draft (lead), writing – review and editing (lead). Victoria Caudill: conceptualization (supporting), investigation (supporting), methodology (supporting), visualization (supporting), writing – original draft (supporting), writing – review and editing (supporting). Samuel E. Champer: conceptualization (supporting), investigation (supporting), methodology (supporting), visualization (supporting), writing – original draft (supporting), writing – review and editing (supporting). Benjamin C. Haller: conceptualization (supporting), methodology (lead), software (lead), supervision (supporting), writing – review and editing (supporting). Clara T. Rehmann: conceptualization (supporting), investigation (supporting), methodology (supporting), visualization (supporting), writing – original draft (supporting), writing – review and editing (supporting). Chris C. R. Smith: conceptualization (supporting), investigation (supporting), methodology (supporting), visualization (supporting), writing – original draft (supporting), writing – review and editing (supporting). Silas Tittes: conceptualization (supporting), investigation (supporting), methodology (supporting), visualization (supporting), writing – original draft (supporting), writing – review and editing (supporting). Philipp W. Messer: funding acquisition (equal), supervision (supporting), writing – review and editing (supporting). Andrew D. Kern: conceptualization (supporting), funding acquisition (equal), supervision (supporting), writing – review and editing (supporting). Sohini Ramachandran: conceptualization (equal), funding acquisition (equal), supervision (equal), writing – review and editing (equal). Peter L. Ralph: conceptualization (equal), funding acquisition (equal), investigation (equal), project administration (equal), software (equal), supervision (equal), visualization (equal), writing – original draft (equal), writing – review and editing (equal).

Acknowledgements

Many thanks to Jeremy Collings and two anonymous reviewers for useful and insightful comments on the manuscript, and to Erin Landguth for input on landscape simulation methods. This work was funded in part by NIH/NHGRI grant R01 HG012473 to P.L.R., A.D.K., B.C.H., and P.W.M., NIH/NIGMS grants R35 GM148253 to A.D.K., R35 GM139628 to S.R. and E.T.C., R35 GM152242 to P.W.M., and F32 GM146484 to C.S. E.T.C. also acknowledges support as a trainee under the Brown University Predoctoral Training Program in Biological Data Science (NIH T32 GM128596).

Conflicts of Interest

The authors declare no conflicts of interest.

Data Availability Statement

SLiM scripts suitable for reuse of all simulations used in this paper are at (for SLiM v4.2 Haller and Messer 2023). Scripts to produce the figures in this manuscript are available at .

References

Adler, P. B., S. P. Ellner, and J. M. Levine. 2010. “Coexistence of Perennial Plants: An Embarrassment of Niches.” Ecology Letters 13, no. 8: 1019–1029. https://doi.org/10.1111/j.1461‐0248.2010.01496.x.

Adler, P. B., A. Kleinhesselink, G. Hooker, J. B. Taylor, B. Teller, and S. P. Ellner. 2018. “Weak Interspecific Interactions in a Sagebrush Steppe? Conflicting Evidence From Observations and Experiments.” Ecology 99, no. 7: 1621–1632. https://doi.org/10.1002/ecy.2363.

Aguillon, S. M., J. W. Fitzpatrick, R. Bowman, et al. 2017. “Deconstructing Isolation‐By‐Distance: The Genomic Consequences of Limited Dispersal.” PLoS Genetics 13, no. 8: e1006911. https://doi.org/10.1371/journal.pgen.1006911.

*Alcock, J. 1987. “Leks and Hilltopping in Insects.” Journal of Natural History 21, no. 2: 319–328. https://doi.org/10.1080/00222938700771041.

Antonio‐Juan Collados, L., S. R. Fassnacht, E. Pardo‐Igúzquiza, and D. Pulido‐Velazquez. 2020. “Assessment of High Resolution Air Temperature Fields at Rocky Mountain National Park by Combining Scarce Point Measurements With Elevation and Remote Sensing Data.” Remote Sensing 13, no. 1: 113.

Aoyama, L., L. G. Shoemaker, B. Gilbert, et al. 2022. “Application of Modern Coexistence Theory to Rare Plant Restoration Provides Early Indication of Restoration Trajectories.” Ecological Applications 32, no. 7: e2649. https://doi.org/10.1002/eap.2649.

Austerlitz, F., C. W. Dick, C. Dutech, et al. 2004. “Using Genetic Markers to Estimate the Pollen Dispersal Curve.” Molecular Ecology 13, no. 4: 937–954. https://doi.org/10.1111/j.1365‐294X.2004.02100.x.

Baddeley, A., E. Rubak, and R. Turner. 2015. Spatial Point Patterns: Methodology and Applications With R. Interdisciplinary Statistics Series. Chapman and Hall/CRC.

Barton, N. H. 1979. “The Dynamics of Hybrid Zones.” Heredity 43, no. 3: 341–359. https://doi.org/10.1038/hdy.1979.87.

Barton, N. H., F. Depaulis, and A. M. Etheridge. 2002. “Neutral Evolution in Spatially Continuous Populations.” Theoretical Population Biology 61, no. 1: 31–48. https://doi.org/10.1006/tpbi.2001.1557.

Barton, N. H., A. M. Etheridge, and A. Véber. 2013. “Modelling evolution in a spatial continuum.” Journal of Statistical Mechanics: Theory and Experiment 2013, no. 1: P01002.

Battey, C. J., P. L. Ralph, and A. D. Kern. 2020. “Space Is the Place: Effects of Continuous Spatial Structure on Analysis of Population Genetic Data.” Genetics 215, no. 1: 1943–2631. https://doi.org/10.1534/genetics.120.303143.

Beever, E. A., C. Ray, P. W. Mote, and J. L. Wilkening. 2010. “Testing Alternative Models of Climate‐Mediated Extirpations.” Ecological Applications 20, no. 1: 164–178.

Bella, I. E. 1971. “A New Competition Model for Individual Trees.” Forest Science 17, no. 3: 364–372. https://doi.org/10.1093/forestscience/17.3.364.

Beverton, R. J. H., and S. J. Holt. 1957. On the Dynamics of Exploited Fish Populations. Springer Netherlands.

Bocedi, G., S. C. F. Palmer, A.‐K. Malchow, D. Zurell, K. Watts, and J. M. J. Travis. 2021. “RangeShifter 2.0: An Extended and Enhanced Platform for Modelling Spatial Eco‐Evolutionary Dynamics and Species' Responses to Environmental Changes.” Ecography 44, no. 10: 1453–1462. https://doi.org/10.1111/ecog.05687.

Bodmer, W. F. 1965. “Differential Fertility in Population Genetics Models.” Genetics 51, no. 3: 411–424. https://doi.org/10.1093/genetics/51.3.411.

Bolker, B. M., and S. W. Pacala. 1997. “Using Moment Equations to Understand Stochastically Driven Spatial Pattern Formation in Ecological Systems.” Theoretical Population Biology 52, no. 3: 179–197.

Booth, T. H., H. A. Nix, M. F. Hutchinson, and T. Jovanic. 1988. “Niche Analysis and Tree Species Introduction.” Forest Ecology and Management 23, no. 1: 47–59. https://doi.org/10.1016/0378‐1127(88)90013‐8.

Bradburd, G. S., and P. L. Ralph. 2019. “Spatial Population Genetics: It's About Time.” Annual Review of Ecology, Evolution, and Systematics 50, no. 1: 427–449. https://doi.org/10.1146/annurev‐ecolsys‐110316‐022659.

Broquet, T., and E. J. Petit. 2009. “Molecular Estimation of Dispersal for Ecology and Population Genetics.” Annual Review of Ecology, Evolution, and Systematics 40, no. 1: 193–216. https://doi.org/10.1146/annurev.ecolsys.110308.120324.

Cain, M. L., B. G. Milligan, and A. E. Strand. 2000. “Long‐Distance Seed Dispersal in Plant Populations.” American Journal of Botany 87, no. 9: 1217–1227.

Canham, C. D., P. T. LePage, and K. D. Coates. 2004. “A Neighborhood Analysis of Canopy Tree Competition: Effects of Shading Versus Crowding.” Canadian Journal of Forest Research 34, no. 4: 778–787. https://doi.org/10.1139/x03‐232.

Cantrell, R. S., and C. Cosner. 2004. Spatial Ecology via Reaction‐Diffusion Equations. John Wiley & Sons.

Caswell, H. 2000. Matrix Population Models. Sinauer.

Champer, S. E., B. Chae, B. C. Haller, J. Champer, and P. W. Messer. 2024. “Resource‐Explicit Interactions in Spatial Population Models.” Methods in Ecology and Evolution 15: 2316–2330. https://doi.org/10.1111/2041‐210X.14432.

Champer, S. E., N. Oakes, R. Sharma, P. García‐Díaz, J. Champer, and P. W. Messer. 2021. “Modeling CRISPR Gene Drives for Suppression of Invasive Rodents Using a Supervised Machine Learning Framework.” PLoS Computational Biology 17, no. 12: e1009660.

Charlesworth, B., D. Charlesworth, and N. H. Barton. 2003. “The Effects of Genetic and Geographic Structure on Neutral Variation.” Annual Review of Ecology, Evolution, and Systematics 34, no. 1: 99–125. https://doi.org/10.1146/annurev.ecolsys.34.011802.132359.

Child, T., B. L. Phillips, and R. Shine. 2009. “Does Desiccation Risk Drive the Distribution of Juvenile Cane Toads (Bufo marinus) in Tropical Australia?” Journal of Tropical Ecology 25, no. 2: 193–200.

Clobert, J., M. Baguette, T. G. Benton, and J. M. Bullock. 2012. Dispersal Ecology and Evolution. Oxford University Press.

Cohen, M. P., and R. A. Alford. 1996. “Factors Affecting Diurnal Shelter Use by the Cane Toad, Bufo marinus.” Herpetologica 20: 172–181.

Connell, J. H. 1971. “On the Role of Natural Enemies in Preventing Competitive Exclusion in Some Marine Animals and in Rain Forest Trees.” Dynamics of Populations 298: 312.

Coulson, T., T. H. G. Ezard, F. Pelletier, et al. 2008. “Estimating the Functional Form for the Density Dependence From Life History Data.” Ecology 89, no. 6: 1661–1674. https://doi.org/10.1890/07‐1099.1.

Courchamp, F., L. Berec, and J. Gascoigne. 2008. Allee Effects in Ecology and Conservation. OUP Oxford.

Cruzan, M. B., and E. C. Hendrickson. 2020. “Landscape Genetics of Plants: Challenges and Opportunities.” Plant Communications 1, no. 6: 100100. https://doi.org/10.1016/j.xplc.2020.100100.

Cushman, S. A., A. J. Shirk, and E. L. Landguth. 2013. “Landscape Genetics and Limiting Factors.” Conservation Genetics 14, no. 2: 263–274. https://doi.org/10.1007/s10592‐012‐0396‐0.

Daniels, R. F. 1976. “Simple Competition Indices and Their Correlation With Annual Loblolly Pine Tree Growth.” Forest Science 22, no. 4: 454–456. https://doi.org/10.1093/forestscience/22.4.454.

De Wit, C. T. 1960. On Competition. Technical Report 66.8. Instituot Voor Biologisch en Scheikundig Onderzoek Van Landbouwgewassen ‐ Wageningen.

DeAngelis, D. L., and S. Yurek. 2017. “Spatially Explicit Modeling in Ecology: A Review.” Ecosystems 20, no. 2: 284–300. https://doi.org/10.1007/s10021‐016‐0066‐z.

*DeVore, J. L., M. R. Crossland, R. Shine, and S. Ducatez. 2021. “The Evolution of Targeted Cannibalism and Cannibal‐Induced Defenses in Invasive Populations of Cane Toads.” Proceedings of the National Academy of Sciences 118, no. 35: e2100765118.

Dieckmann, U., and R. Law. 1999. “Relaxation Projections and the Method of Moments.” In The Geometry of Ecological Interactions: Simplifying Spatial Complexity, edited by U. Dieckmann, R. Law, and J. A. Metz, 412–455. Cambridge University Press.

Dobzhansky, T., and S. Wright. 1943. “Genetics of Natural Populations. X. Dispersion Rates in Drosophila Pseudoobscura.” Genetics 28, no. 4: 304–340.

Dunning, J. B., Jr., D. J. Stewart, B. J. Danielson, et al. 1995. “Spatially Explicit Population Models: Current Forms and Future Uses.” Ecological Applications 5, no. 1: 3–11. https://doi.org/10.2307/1942045.

Edelaar, P., and D. I. Bolnick. 2012. “Non‐Random Gene Flow: An Underappreciated Force in Evolution and Ecology.” Trends in Ecology & Evolution 27, no. 12: 659–665. https://doi.org/10.1016/j.tree.2012.07.009.

Edwards, A. M., R. A. Phillips, N. W. Watkins, et al. 2007. “Revisiting Lévy Flight Search Patterns of Wandering Albatrosses, Bumblebees and Deer.” Nature 449, no. 7165: 1044–1048. https://doi.org/10.1038/nature06199.

Ellner, S. P., D. Z. Childs, and M. Rees. 2016. Data‐Driven Modelling of Structured Populations: A Practical Guide to the Integral Projection Model. Springer International Publishing.

Endler, J. A. 1986. Natural Selection in the Wild. Princeton University Press.

Epperson, B. K. 2003. Geographical Genetics. Monographs in Population Biology. Princeton University Press.

*Erb, L. P., C. Ray, and R. Guralnick. 2014. “Determinants of Pika Population Density vs. Occupancy in the Southern Rocky Mountains.” Ecological Applications 24, no. 3: 429–435.

Eskola, H. T. M., and S. A. H. Geritz. 2007. “On the Mechanistic Derivation of Various Discrete‐Time Population Models.” Bulletin of Mathematical Biology 69, no. 1: 329–346. https://doi.org/10.1007/s11538‐006‐9126‐4.

Estep, L. K., N. D. Burkett‐Cadena, G. E. Hill, R. S. Unnasch, and T. R. Unnasch. 2010. “Estimation of Dispersal Distances of Culex erraticus in a Focus of Eastern Equine Encephalitis Virus in the Southeastern United States.” Journal of Medical Entomology 47, no. 6: 977–986. https://doi.org/10.1603/ME10056.

*Estoup, A., S. J. E. Baird, N. Ray, et al. 2010. “Combining Genetic, Historical and Geographical Data to Reconstruct the Dynamics of Bioinvasions: Application to the Cane Toad Bufo marinus.” Molecular Ecology Resources 10, no. 5: 886–901.

Etheridge, A. M., T. G. Kurtz, I. Letter, P. L. Ralph, and T. T. H. Lung. 2024. “Looking Forwards and Backwards: Dynamics and Genealogies of Locally Regulated Populations.” Electronic Journal of Probability 29: 1–85. https://doi.org/10.1214/24‐EJP1075.

Etheridge, A., and S. Penington. 2022. “Genealogies in Bistable Waves.” Electronic Journal of Probability 27: 1–99. https://doi.org/10.1214/22‐EJP845.

European Space Agency. 2024. “Earth System Data Lab.” https://www.earthsystemdatalab.net/.

Felsenstein, J. 1975. “A Pain in the Torus: Some Difficulties With Models of Isolation by Distance.” American Naturalist 109, no. 967: 359–368.

Felsenstein, J. 1976. “The Theoretical Population Genetics of Variable Selection and Migration.” Annual Review of Genetics 10: 253–280. https://doi.org/10.1146/annurev.ge.10.120176.001345.

*Fick, S. E., and R. J. Hijmans. 2017. “WorldClim 2: New 1‐Km Spatial Resolution Climate Surfaces for Global Land Areas.” International Journal of Climatology 37, no. 12: 4302–4315.

Foster, G., and S. Rahmstorf. 2011. “Global Temperature Evolution 1979–2010.” Environmental Research Letters 6, no. 4: 44022.

Freckleton, R. P., A. R. Watkinson, R. E. Green, and W. J. Sutherland. 2006. “Census Error and the Detection of Density Dependence.” Journal of Animal Ecology 75, no. 4: 837–851. https://doi.org/10.1111/j.1365‐2656.2006.01121.x.

Freeland, W. J. 1986. “Populations of Cane Toad, Bufo marinus, in Relation to Time Since Colonization.” Wildlife Research 13, no. 2: 321–329.

French, N. R., T. Y. Tagami, and P. Hayden. 1968. “Dispersal in a Population of Desert Rodents.” Journal of Mammalogy 49, no. 2: 272–280.

Geritz, S. A. H., and É. Kisdi. 2012. “Mathematical Ecology: Why Mechanistic Models?” Journal of Mathematical Biology 65, no. 6: 1411–1415.

Graves, T. A., P. Beier, and J. A. Royle. 2013. “Current Approaches Using Genetic Distances Produce Poor Estimates of Landscape Resistance to Interindividual Dispersal.” Molecular Ecology 22, no. 15: 3888–3903.

Grimm, V. 1999. “Ten Years of Individual‐Based Modelling in Ecology: What Have We Learned and What Could We Learn in the Future?” Ecological Modelling 115, no. 2–3: 129–148.

Haller, B. C., and P. W. Messer. 2023. “SLiM 4: Multispecies Eco‐Evolutionary Modeling.” American Naturalist 201, no. 5: E127–E139. https://doi.org/10.1086/723601.

Haller, B. C., and P. W. Messer. 2024. “SLiM: An Evolutionary Simulation Framework (The SLiM Manual).” https://messerlab.org/SLiM.

Hammond, R. L., L. J. Lawson Handley, B. J. Winney, M. W. Bruford, and N. Perrin. 2005. “Genetic Evidence for Female‐Biased Dispersal and Gene Flow in a Polygynous Primate.” Proceedings of the Royal Society B: Biological Sciences 273, no. 1585: 479–484. https://doi.org/10.1098/rspb.2005.3257.

Hanski, I. 1997. “Metapopulation Dynamics: From Concepts and Observations to Predictive Models.” In Metapopulation Biology, edited by I. Hanski and M. E. Gilpin, 69–91. Academic Press. https://doi.org/10.1016/B978‐012323445‐2/50007‐9.

Harman, R. R., J. Goddard, R. Shivaji, and J. T. Cronin. 2020. “Frequency of Occurrence and Population‐Dynamic Consequences of Different Forms of Density‐Dependent Emigration.” American Naturalist 20: 8156.

Hewitt, G. M. 2011. “Quaternary Phylogeography: The Roots of Hybrid Zones.” Genetica 139, no. 5: 617–638. https://doi.org/10.1007/s10709‐011‐9547‐3.

Hülsmann, L., R. A. Chisholm, and F. Hartig. 2021 01695347. “Is Variation in Conspecific Negative Density Dependence Driving Tree Diversity Patterns at Large Scales?” Trends in Ecology & Evolution 36, no. 2: 151–163. https://doi.org/10.1016/j.tree.2020.10.003.

Janzen, D. H. 1970. “Herbivores and the Number of Tree Species in Tropical Forests.” American Naturalist 104, no. 940: 501–528. https://doi.org/10.1086/282687.

*Jordahl, K., J. Van den Bossche, M. Fleischmann, et al. 2020. “geopandas/geopandas: v0.8.1.” https://doi.org/10.5281/zenodo.3946761.

Katul, G. G., A. Porporato, R. Nathan, et al. 2005. “Mechanistic Analytical Models for Long‐Distance Seed Dispersal by Wind.” American Naturalist 166, no. 3: 368–381. https://doi.org/10.1086/432589.

*Kearney, M., B. L. Phillips, C. R. Tracy, K. A. Christian, G. Betts, and W. P. Porter. 2008. “Modelling Species Distributions Without Using Species Distributions: The Cane Toad in Australia Under Current and Future Climates.” Ecography 31, no. 4: 423–434.

*Kotz, S., and S. Nadarajah. 2004. Multivariate T‐Distributions and Their Applications. Cambridge University Press.

Lande, R., S. Engen, and B.‐E. Saether. 2003. Stochastic Population Dynamics in Ecology and Conservation. Oxford University Press.

Landguth, E. L., A. Bearlin, C. C. Day, and J. Dunham. 2017. “CDMetaPOP: An Individual‐Based, Eco‐Evolutionary Model for Spatially Explicit Simulation of Landscape Demogenetics.” Methods in Ecology and Evolution 8, no. 1: 4–11. https://doi.org/10.1111/2041‐210X.12608.

Landguth, E. L., and S. A. Cushman. 2010. “CDPOP: A Spatially Explicit Cost Distance Population Genetics Program.” Molecular Ecology Resources 10, no. 1: 156–161. https://doi.org/10.1111/j.1755‐0998.2009.02719.x.

Laporte, V., and B. Charlesworth. 2002. “Effective Population Size and Population Subdivision in Demographically Structured Populations.” Genetics 162, no. 1: 501–519. https://doi.org/10.1093/genetics/162.1.501.

*Lloyd, M. 1967. “Mean Crowding.” Journal of Animal Ecology 36, no. 1: 1–30.

Lundgren, E., and P. L. Ralph. 2019. “Are Populations Like a Circuit? Comparing Isolation by Resistance to a New Coalescent‐Based Method.” Molecular Ecology Resources 19, no. 6: 1388–1406. https://doi.org/10.1111/1755‐0998.13035.

Mack, R. N., and J. L. Harper. 1977. “Interference in Dune Annuals: Spatial Pattern and Neighbourhood Effects.” Journal of Ecology 65, no. 2: 345–363.

Malécot, G. 1969. The Mathematics of Heredity. Freeman.

Mallet, J. 2012. “The Struggle for Existence. How the Notion of Carrying Capacity, k, Obscures the Links Between Demography, Darwinian Evolution and Speciation.” Evolutionary Ecology Research 20: 1860.

Malyon, B., M. M. Mayfield, T. E. Martyn, and D. B. Stouffer. 2023. “Estimating Interaction Strengths for Diverse Horizontal Systems Using Performance Data.” Methods in Ecology and Evolution 14, no. 3: 968–980. https://doi.org/10.1111/2041‐210X.14068.

Mazzucco, R., M. Doebeli, and U. Dieckmann. 2018. “The Influence of Habitat Boundaries on Evolutionary Branching Along Environmental Gradients.” Evolutionary Ecology 32, no. 6: 563–585. https://doi.org/10.1007/s10682‐018‐9956‐1.

McRae, B. H., B. G. Dickson, T. H. Keitt, and V. B. Shah. 2008. “Using Circuit Theory to Model Connectivity in Ecology, Evolution, and Conservation.” Ecology 89, no. 10: 2712–2724. https://doi.org/10.1890/07‐1861.1.

*Miller, T. E., A. L. Angert, C. D. Brown, et al. 2020. “Eco‐Evolutionary Dynamics of Range Expansion.” Ecology 101, no. 10: e03139. https://doi.org/10.1002/ecy.3139.

Min, J., M. Gupta, M. M. Desai, and D. B. Weissman. 2022. “Spatial Structure Alters the Site Frequency Spectrum Produced by Hitchhiking.” Genetics 222, no. 3: 139. https://doi.org/10.1093/genetics/iyac139.

Mirski, P., E. Brzosko, I. Jedrzejczyk, J. Kotowicz, B. Ostrowiecka, and A. Wróblewska. 2017. “Genetic Structure of Dioecious and Trioecious Salix myrsinifolia Populations at the Border of Geographic Range.” Tree Genetics & Genomes 13, no. 1: 15. https://doi.org/10.1007/s11295‐016‐1096‐6.

Mollison, D. 1972. “The Rate of Spatial Propagation of Simple Epidemics.” In Proceedings of the Sixth Berkeley Symposium on Mathematical Statistics and Probability (Univ. California, Berkeley, Calif., 1970/1971), Vol. III: Probability Theory, 579–614. University of California Press.

Montero, D., C. Aybar, M. D. Mahecha, F. Martinuzzi, M. Söchting, and S. Wieneke. 2023. “A Standardized Catalogue of Spectral Indices to Advance the Use of Remote Sensing in Earth System Research.” Scientific Data 10, no. 1: 197. https://doi.org/10.1038/s41597‐023‐02096‐0.

Montgomery, S. 1973. “Gene Flow and Selection in a Cline.” Genetics 75, no. 4: 733–756.

Morales, J. M., and T. Morán López. 2022. “Mechanistic Models of Seed Dispersal by Animals.” Oikos 2022, no. 2: 8328. https://doi.org/10.1111/oik.08328.

Murray, D. R. 2012. Seed Dispersal by Water, Chapter 2, 49–86. Academic Press.

Nagylaki, T., and J. F. Crow. 1974. “Continuous Selective Models.” Theoretical Population Biology 5, no. 2: 257–283. https://doi.org/10.1016/0040‐5809(74)90045‐8.

NASA. 2024. “Earthdata.” https://www.earthdata.nasa.gov/.

Neubert, M. G., M. Kot, and M. A. Lewis. 1995. “Dispersal and Pattern Formation in a Discrete‐Time Predator‐Prey Model.” Theoretical Population Biology 48, no. 1: 7–43. https://doi.org/10.1006/tpbi.1995.1020.

Neuenschwander, S., F. Michaud, and J. Goudet. 2018. “QuantiNemo 2: A Swiss Knife to Simulate Complex Demographic and Genetic Scenarios, Forward and Backward in Time.” Bioinformatics 35, no. 5: 886–888. https://doi.org/10.1093/bioinformatics/bty737.

North, A. R., A. Burt, and H. C. J. Godfray. 2019. “Modelling the Potential of Genetic Control of Malaria Mosquitoes at National Scale.” BMC Biology 17, no. 1: 1–12. https://doi.org/10.1186/s12915‐019‐0645‐5.

Nullmeier, J., and O. Hallatschek. 2013. “The Coalescent in Boundary‐Limited Range Expansions.” Evolution 67, no. 5: 1307–1320. https://doi.org/10.1111/evo.12037.

Nurminiemi, M., J. Tufto, N.‐O. Nilsson, and O. A. Rognli. 1998. “Spatial Models of Pollen Dispersal in the Forage Grass Meadow Fescue.” Evolutionary Ecology 12, no. 4: 487. https://doi.org/10.1023/A:1006529023036.

Obbard, D. J., S. A. Harris, and J. R. Pannell. 2006. “Sexual Systems and Population Genetic Structure in an Annual Plant: Testing the Metapopulation Model.” American Naturalist 167, no. 3: 354–366. https://doi.org/10.1086/499546.

Oberhauser, K. S., M. J. Solensky, Journey North, and Monarch Larva Monitoring Project. 2004. Monarch Butterfly Biology & Conservation. Cornell University Press.

Pacala, S. W. 1986. “Neighborhood Models of Plant Population Dynamics. 2. Multi‐Species Models of Annuals.” Theoretical Population Biology 29, no. 2: 262–292. https://doi.org/10.1016/0040‐5809(86)90011‐0.

Pacala, S. W., and J. A. Silander. 1985. “Neighborhood Models of Plant Population Dynamics. I. Single‐Species Models of Annuals.” American Naturalist 125, no. 3: 385–411.

Paradis, E., S. R. Baillie, W. J. Sutherland, and R. D. Gregory. 1998. “Patterns of Natal and Breeding Dispersal in Birds.” Journal of Animal Ecology 67, no. 4: 518–536. https://doi.org/10.1046/j.1365‐2656.1998.00215.x.

Paulose, J., J. Hermisson, and O. Hallatschek. 2019. “Spatial Soft Sweeps: Patterns of Adaptation in Populations With Long‐Range Dispersal.” PLoS Genetics 15, no. 2: 1–29 02 2019. https://doi.org/10.1371/journal.pgen.1007936.

Pearson, O. P. 1960. “A Mechanical Model for the Study of Population Dynamics.” Ecology 41, no. 3: 494–508.

Peterson, A. T. 2001. “Predicting Species' Geographic Distributions Based on Ecological Niche Modeling.” Condor 103, no. 3: 599–605.

Peterson, E. E., E. M. Hanks, M. B. Hooten, J. M. Ver Hoef, and M.‐J. Fortin. 2019. “Spatially Structured Statistical Network Models for Landscape Genetics.” Ecological Monographs 89, no. 2: e01355. https://doi.org/10.1002/ecm.1355.

Pons, J., and J. G. Pausas. 2007. “Acorn Dispersal Estimated by Radio‐Tracking.” Oecologia 153, no. 4: 903–911. https://doi.org/10.1007/s00442‐007‐0788‐x.

PRISM Climate Group. 2024. “PRISM Climate Data.” https://www.prism.oregonstate.edu/.

Pusey, A. E., and C. Packer. 1987. “The Evolution of Sex‐Biased Dispersal in Lions.” Behaviour 101, no. 4: 275–310. https://doi.org/10.1163/156853987X00026.

Ramachandran, S., O. Deshpande, C. C. Roseman, N. A. Rosenberg, M. W. Feldman, and L. L. Cavalli‐Sforza. 2005. “Support From the Relationship of Genetic and Geographic Distance in Human Populations for a Serial Founder Effect Originating in Africa.” Proceedings of the National Academy of Sciences 102, no. 44: 15942–15947. https://doi.org/10.1073/pnas.0507611102.

Rebaudo, F., A. Le Rouzic, S. Dupas, J.‐F. Silvain, M. Harry, and O. Dangles. 2013. “SimAdapt: An Individual‐Based Genetic Model for Simulating Landscape Management Impacts on Populations.” Methods in Ecology and Evolution 4, no. 6: 595–600. https://doi.org/10.1111/2041‐210X.12041.

Richardson, L. K., S. W. Nordstrom, A. Waananen, et al. 2024. “Juvenile Survival Increases With Dispersal Distance and Varies Across Years: 15 Years of Evidence in a Prairie Perennial.” Ecology 20: 4331. https://doi.org/10.1002/ecy.4331.

Robledo‐Arnuncio, J. J., and L. Gil. 2005. “Patterns of Pollen Dispersal in a Small Population of Pinus sylvestris L. Revealed by Total‐Exclusion Paternity Analysis.” Heredity 94, no. 1: 13–22. https://doi.org/10.1038/sj.hdy.6800542.

Robledo‐Arnuncio, J. J., and F. Rousset. 2010. “Isolation by Distance in a Continuous Population Under Stochastic Demographic Fluctuations.” Journal of Evolutionary Biology 23, no. 1: 53–71. https://doi.org/10.1111/j.1420‐9101.2009.01860.x.

Robledo‐Arnuncio, J. J., and C. García. 2007. “Estimation of the Seed Dispersal Kernel From Exact Identification of Source Plants.” Molecular Ecology 16, no. 23: 5098–5109. https://doi.org/10.1111/j.1365‐294X.2007.03427.x.

Rogers, H. S., N. G. Beckman, F. Hartig, et al. 2019. “The Total Dispersal Kernel: A Review and Future Directions.” AoB Plants 11, no. 5: 42. https://doi.org/10.1093/aobpla/plz042.

Rosenberg, N. A., S. Mahajan, S. Ramachandran, C. Zhao, J. K. Pritchard, and M. W. Feldman. 2005. “Clines, Clusters, and the Effect of Study Design on the Inference of Human Population Structure.” PLoS Genetics 1, no. 6: e70. https://doi.org/10.1371/journal.pgen.0010070.

Rousset, A. 2000. “Genetic Differentiation Between Individuals.” Journal of Evolutionary Biology 13, no. 1: 58–62. https://doi.org/10.1046/j.1420‐9101.2000.00137.x.

Rousset, F. 1997. “Genetic Differentiation and Estimation of Gene Flow From F‐Statistics Under Isolation by Distance.” Genetics 145, no. 4: 1219–1228.

Saastamoinen, M., G. Bocedi, J. Cote, et al. 2018. “Genetics of Dispersal.” Biological Reviews of the Cambridge Philosophical Society 93, no. 1: 574–599.

*Sasaki, A. 1997. “Clumped Distribution by Neighbourhood Competition.” Journal of Theoretical Biology 186, no. 4: 415–430.

Savolainen, O., T. Pyhäjärvi, and T. Knürr. 2007. “Gene Flow and Local Adaptation in Trees.” Annual Review of Ecology, Evolution, and Systematics 38, no. 1: 595–619. https://doi.org/10.1146/annurev.ecolsys.38.091206.095646.

Schumaker, N. H., and A. Brookes. 2018. “HexSim: A Modeling Environment for Ecology and Conservation.” Landscape Ecology 33, no. 2: 197–211. https://doi.org/10.1007/s10980‐017‐0605‐9.

Sedghifar, A., Y. Brandvain, and P. Ralph. 2016. “Beyond Clines: Lineages and Haplotype Blocks in Hybrid Zones.” Molecular Ecology 25, no. 11: 2559–2576. https://doi.org/10.1111/mec.13677.

Shine, R. 2010. “The Ecological Impact of Invasive Cane Toads (Bufo marinus) in Australia.” Quarterly Review of Biology 85, no. 3: 253–291.

Shine, R., R. A. Alford, R. Blennerhasset, et al. 2021. “Increased Rates of Dispersal of Free‐Ranging Cane Toads (Rhinella marina) During Their Global Invasion.” Scientific Reports 11, no. 1: 23574.

Shine, R., G. Ward‐Fear, and G. P. Brown. 2020. “A Famous Failure: Why Were Cane Toads an Ineffective Biocontrol in Australia?” Conservation Science and Practice 2, no. 12: e296.

Silander, J. A., and S. W. Pacala. 1985. “Neighborhood Predictors of Plant Performance.” Oecologia 66, no. 2: 256–263.

Sillero, N. 2011. “What Does Ecological Modelling Model? A Proposed Classification of Ecological Niche Models Based on Their Underlying Methods.” Ecological Modelling 222, no. 8: 1343–1346. https://doi.org/10.1016/j.ecolmodel.2011.01.018.

Sillero, N., and A. M. Barbosa. 2021. “Common Mistakes in Ecological Niche Models.” International Journal of Geographical Information Science 35, no. 2: 213–226. https://doi.org/10.1080/13658816.2020.1798968.

Skellam, J. G. 1951. “Random Dispersal in Theoretical Populations.” Biometrika 38: 196–218.

*Smith, A. T. 1974. “The Distribution and Dispersal of Pikas: Consequences of Insular Population Structure.” Ecology 55, no. 5: 1112–1119.

Smith, C. C. R., G. Patterson, P. L. Ralph, and A. D. Kern. 2024. “Estimation of Spatial Demographic Maps From Polymorphism Data Using a Neural Network.” bioRxiv 20: 5300. https://doi.org/10.1101/2024.03.15.585300.

Smith, C. C. R., S. Tittes, P. L. Ralph, and A. D. Kern. 2023. “Dispersal Inference From Population Genetic Variation Using a Convolutional Neural Network.” Genetics 224, no. 2: iyad068. https://doi.org/10.1093/genetics/iyad068.

Smith, T. B., and D. B. Weissman. 2023. “Isolation by Distance in Populations With Power‐Law Dispersal.” G3 (Bethesda, Md.) 13, no. 4: 23. https://doi.org/10.1093/g3journal/jkad023.

Snyder, R. E., and P. Chesson. 2004. “How the Spatial Scales of Dispersal, Competition, and Environmental Heterogeneity Interact to Affect Coexistence.” American Naturalist 164, no. 5: 633–650. https://doi.org/10.1086/424969.

Song, X., J. Y. Lim, J. Yang, and M. S. Luskin. 2021. “When Do Janzen‐Connell Effects Matter? A Phylogenetic Meta‐Analysis of Conspecific Negative Distance and Density Dependence Experiments.” Ecology Letters 24, no. 3: 608–620. https://doi.org/10.1111/ele.13665.

Specht, H. M., and T. W. Arnold. 2018. “Banding Age Ratios Reveal Prairie Waterfowl Fecundity Is Affected by Climate, Density Dependence and Predator–Prey Dynamics.” Journal of Applied Ecology 55, no. 6: 2854–2864. https://doi.org/10.1111/1365‐2664.13186.

Spotswood, E. N., P. Mariotte, E. C. Farrer, L. Nichols, and K. N. Suding. 2017. “Separating Sources of Density‐Dependent and Density‐Independent Establishment Limitation in Invading Species.” Journal of Ecology 105, no. 2: 436–444. https://doi.org/10.1111/1365‐2745.12686.

Shuster, S. M. 2003. “Mating Systems and Strategies.” In Monographs in Behavior and Ecology. University Press.

Stillman, R. A., S. F. Railsback, J. Giske, U. Berger, and V. Grimm. 2015. “Making Predictions in a Changing World: The Benefits of Individual‐Based Ecology.” Bioscience 65, no. 2: 140–150. https://doi.org/10.1093/biosci/biu192.

Suchan, T., C. P. Bataille, M. S. Reich, et al. 2024. “A Trans‐Oceanic Flight of Over 4,200 Km by Painted Lady Butterflies.” Nature Communications 15, no. 1: 5205.

*Taylor, P. 2010. “Birth–Death Symmetry in the Evolution of a Social Trait.” Journal of Evolutionary Biology 23, no. 12: 2569–2578. https://doi.org/10.1111/j.1420‐9101.2010.02122.x.

Teller, B. J., P. B. Adler, C. B. Edwards, G. Hooker, and S. P. Ellner. 2016. “Linking Demography With Drivers: Climate and Competition.” Methods in Ecology and Evolution 7, no. 2: 171–183. https://doi.org/10.1111/2041‐210X.12486.

Terborgh, J. 2012. “Enemies Maintain Hyperdiverse Tropical Forests.” American Naturalist 179, no. 3: 303–314. https://doi.org/10.1086/664183.

Travis, J., R. D. Bassar, T. Coulson, D. Reznick, and M. Walsh. 2023. “Density‐Dependent Selection.” Annual Review of Ecology, Evolution, and Systematics 54, no. 1: 85–105. https://doi.org/10.1146/annurev‐ecolsys‐110321‐055345.

Trochet, A., E. A. Courtois, V. M. Stevens, et al. 2016. “Evolution of Sex‐Biased Dispersal.” Quarterly Review of Biology 91, no. 3: 297–320. https://doi.org/10.1086/688097.

Tufto, J., S. Engen, and K. Hindar. 1997. “Stochastic Dispersal Processes in Plant Populations.” Theoretical Population Biology 52, no. 1: 16–26.

Tuljapurkar, S. 1989. “An Uncertain Life: Demography in Random Environments.” Theoretical Population Biology 35, no. 3: 227–294. https://doi.org/10.1016/0040‐5809(89)90001‐4.

Tuljapurkar, S. 2013. Population Dynamics in Variable Environments. Lecture Notes in Biomathematics. Springer.

Uecker, H., S. P. Otto, and J. Hermisson. 2014. “Evolutionary Rescue in Structured Populations.” American Naturalist 183, no. 1: 17–35. https://doi.org/10.1086/673914.

Urban, M. C., B. L. Phillips, D. K. Skelly, and R. Shine. 2008. “A Toad More Traveled: The Heterogeneous Invasion Dynamics of Cane Toads in Australia.” American Naturalist 171, no. 3: E134–E148.

Wallace, B. 1975. “Hard and Soft Selection Revisited.” Evolution 29, no. 3: 465–473.

Weigelt, A., and P. Jolliffe. 2003. “Indices of Plant Competition.” Journal of Ecology 91, no. 5: 707–720.

Weiner, J. 1982. “A Neighborhood Model of Annual‐Plant Interference.” Ecology 63, no. 5: 1237–1241.

Weiss‐Lehman, C. P., C. M. Werner, C. H. Bowler, et al. 2022. “Disentangling Key Species Interactions in Diverse and Heterogeneous Communities: A Bayesian Sparse Modelling Approach.” Ecology Letters 25, no. 5: 1263–1276. https://doi.org/10.1111/ele.13977.

Werkowska, W., A. L. Márquez, R. Real, and P. Acevedo. 2017. “A Practical Overview of Transferability in Species Distribution Modeling.” Environmental Reviews 25, no. 1: 127–133. https://doi.org/10.1139/er‐2016‐0045.

Wright, S. 1943. “Isolation by Distance.” Genetics 28: 114–138.

Wright, S. 1946. “Isolation by Distance Under Diverse Systems of Mating.” Genetics 31, no. 1: 39–59.

WWF Australia. 2019. “10 Facts About Cane Toads.” https://wwf.org.au/blogs/10‐facts‐about‐cane‐toads/.

Yates, K. L., P. J. Bouchet, M. J. Caley, et al. 2018. “Outstanding Challenges in the Transferability of Ecological Models.” Trends in Ecology & Evolution 33, no. 10: 790–802. https://doi.org/10.1016/j.tree.2018.08.001.

Zaiats, A., M. J. Germino, M. D. Serpe, B. A. Richardson, and T. Trevor Caughlin. 2021. “Intraspecific Variation Mediates Density Dependence in a Genetically Diverse Plant Species.” Ecology 102, no. 11: e03502.

Zipkin, E. F., J. W. Doser, C. L. Davis, W. Leuenberger, S. Ayebare, and K. L. Davis. 2023. “Integrated Community Models: A Framework Combining Multispecies Data Sources to Estimate the Status, Trends and Dynamics of Biodiversity.” Journal of Animal Ecology 92, no. 12: 2248–2262. https://doi.org/10.1111/1365‐2656.14012.