Stochastic diffusion using mean-field limits to approximate master equations

Authors: Laurent Hébert-Dufresne, Matthew M. Kling, Samuel F. Rosenblatt, Stephanie N. Miller, P. Alexander Burnham, Nicholas W. Landry, Nicholas J. Gotelli, Brian J. McGill

Abstract

Stochastic diffusion is the noisy process through which dynamics like epidemics, or agents like animal species, disperse over a larger area. These processes are increasingly important to better prepare for pandemics and as species ranges shift in response to climate change. Unfortunately, modelling is mostly done with expensive computational simulations or inaccurate deterministic tools that ignore the randomness of dispersal. We introduce ‘mean-FLAME’ models, tracking stochastic dispersion using approximate master equations to follow the probability distribution over all possible states of an area of interest, up to states active enough to be approximated using a mean-field model. In the limit where we track all states, this approach is locally exact, and in the other limit collapses to traditional deterministic models. In predator–prey systems, we show that tracking a handful of states around key absorbing states is sufficient to accurately model extinction. In disease models, we show that classic mean-field approaches underestimate the heterogeneity of epidemics. And in nonlinear dispersal models, we show that deterministic tools fail to capture the speed of spatial diffusion. These effects are all important for marginal areas that are close to unsuitable for diffusion, like the edge of a species range or epidemics in small populations.

Introduction

Self-organizing structures (e.g. populations, metapopulations, communities) and their complex dynamics (e.g. epidemics, species diffusion, culture) involve a combination of deterministic trends and stochastic effects. Although most ecological modelling efforts have focused on the deterministic component, the stochastic component can be critical when considering the movement of entities through time and space [1]. Models of spread or diffusion, in particular, often rely on approximations that are built from deterministic assumptions [2,3], like supposing well-mixed populations [4] or by simply tracking the average state of a system [5,6]. These approximation methods work well if the system is large enough or if the probabilities of rare events are low enough. Still, they often do a poor job representing non-deterministic processes near the limits of systems (e.g. small populations, populations near the edge of their geographic range). Here, we present a novel method that can better approximate the state of ecological systems under near-limit conditions.

First, consider the spread of an emergent disease under differently structured populations. In a large, well-connected population, a virulent disease could spread quickly, moving among dense population centres and eventually among sparser rural areas once transmission has exceeded a high enough rate [7]. Alternatively, in a network of small, isolated populations (e.g. bark beetles moving among stands of their host trees), transmission to new areas occurs more slowly, if at all, and relies on rare events. Models that ignore stochasticity will capture the dynamics of the large connected population but fail to capture the sparse stochastic transmission chains critical to the small isolated populations. With close to half of the global population living in rural areas, this is a large omission for many state-of-the-art epidemic models and the ecology of emerging diseases [8].

Second, under climate change, organisms are expected to track suitable conditions by shifting their geographic ranges [9]. For the range of a species to shift, individuals need suitable conditions to effectively survive and reproduce, and need to reach locations with those conditions. In cases where distances are short and there are no physical barriers between current occurrences and future suitable locations, modelling these range shifts using deterministic mean-field averages may be appropriate. However, in cases where long-distance dispersal is required to reach suitable conditions or get beyond barriers (e.g. mountain ranges, rivers, highways), stochastic events of dispersal are more important. Rare and random dispersal events have shaped species ranges as we know them today [10,11]. Capturing rare and random dispersal is therefore critical to understanding how ecosystems will shift in response to climate change [12,13].

To better model the stochastic components like those critical to species range shifts and the spread of emergent diseases, we propose using a new approach based on the master equation framework. The objective of this study is to explore the potential of using a coupled system of mean-field and master equations to efficiently capture the discrete and stochastic nature of fluctuations in spreading and dispersion dynamics. This hybrid approach aims to maintain the precision of master equations around important absorbing states while leveraging the fact that the complexity of mean-field models does not scale with population size, unlike classic master equation or simulation frameworks. We do so by introducing mean-field limits to approximate master equations (or mean-FLAME). The methods we propose can help improve the computational and mathematical modelling of many applications, some of which we tackle

• —population dynamics,
• —consumer–resource or predator–prey model,
• —metapopulation models,
• —epidemic dynamics,
• —age- and stage-structured models,
• —species range shifts,
• —competition or cooperation among species,
• —evolutionary dynamics,
• —trophic interactions,
• —and invasive species dynamics.

and others left for future

• —competition or cooperation among species,
• —evolutionary dynamics,
• —trophic interactions,
• —and invasive species dynamics.

The rest of the paper is structured as follows. We first introduce the concept of master equations and mean-field limits without space, using a simple birth–death process for a single population and a coupled two-population system, in §2. We then extend this method to a two-dimensional process using Lotka–Volterra predator–prey dynamics, again without space, in §3. We then apply the method to spatially embedded multidimensional examples in §§4 and 5, using the above scenarios of metapopulation epidemic dynamics and species dispersal, respectively. We conclude by sharing some details of our software implementations of mean-FLAME and suggesting improvements for future work.

Modelling approach with a birth–death process example

In this section, we use a stochastic birth–death process for a single population to describe three components of our proposed methodology. First, we present the traditional master equation framework which tracks the probability distribution of a system over all of its possible states by accounting for transition probability between discrete states. Next, we use a mean-field limit to approximate the master equation far from any absorbing state where dynamics die (for population size far above zero in this birth–death example). Then we show how approximating correlations across sub-systems or locations can simplify the processing of the master equation for coupled or spatially embedded systems.

For this example, we are interested in the number n of individuals found in a population (system) at time t. In this example, births and deaths are stochastic events that follow a Poisson process with births occurring at a rate μn and deaths at rate νn2. With this formulation, n∗=0 is an absorbing state that systems cannot leave (i.e., extinction) and there exists a meta-stable state at n∗∗=μ/ν.

In a classic mean-field description, we would track the expected number of individuals n¯ at time t. This would traditionally be done with the following ordinary differential

With this equation, the absorbing state n∗=0 is always unstable if μ>0 such that any system with initial conditions n¯(0)>0 will end up in the stable state n∗∗=μ/ν. However, consider an actual stochastic simulation of a system with n=1 active individual; there is a probability ν/(μ+ν) that the system falls into the absorbing state instead of moving to n=2. Mean-field approximations do not capture these stochastic extinction events. In the ecological literature, stochastic outcomes like this are referred to as demographic or environmental stochasticity and can be important drivers of population change [14]. In fact, all populations could eventually fall into the absorbing state, and the rate at which these extinction events occur will depend on how close n∗∗=μ/ν is to zero. We wish to capture these different paths to extinction.

Master equation

Master equations are systems of equations that track the occupation probability (or occupation number) for every possible state of a system by calculating the exact transition rates between states [15]. Consider Py→(t), the occupation probability of state y→ at time t. The state y→ in this example is a simple count of individuals but it could include any description of the state of the system (e.g. location, epidemiological status or age of individuals, local policies), as long as y→ is sufficient to exactly specify the transition rates ω(y→,z→) from state y→ to other states z→ and vice versa. We can then define a general master equation for Py→(t) as,

where positive terms correspond to probability flowing into state y→, from z→ and negative terms to probability flowing from y→ to other states. Often, the sums in equation (2.2) involve only a few terms as the transition matrix across all states can be quite sparse. The complexity of the approach is instead related to the number of states available to the system, since the approach is only exact if all states are tracked explicitly.

In general, the master equation describing a given stochastic process can be written by answering a simple set of questions.

• (1) What is the set of discrete states available to the system? The answer to this questions fixes the number of variables and differential equations needed.
• (2) What are the possible transitions between states? The answer to this question dictates how many terms with other states z→ should appear in a given differential equation.
• (3) At what rates do these transitions occur? These rates fix the transition factors ω in the master equation.

To write the master equation of this birth–death process, we follow our simple recipe. The system consists of a discrete number of active individuals n which can take any value from 0 to infinity. The system can go from state n to n+1 through reproduction or to state n−1 when an individual dies. The reproduction transition from n to n+1 occurs at rate ω(n,n+1)=μn and the death transition to n−1 occurs at rate ω(n,n−1)=νn2 for all possible values of n; all other transitions are impossible. The corresponding master equation is therefore

Master equation with mean-field limit

One possible issue with the master equation framework is that we need to track all possible states for the system. What if the population of active individuals goes to infinity? We can track thousands of equations without issue, but systems of several millions or billions of equations are computationally expensive.

Therefore, we introduce a mean-field limit in our master equation. We collapse all states with intense activity as a mean-field quantity I(t). Technically, we define intense activity as a system far enough from any absorbing state. In this birth–death example, this would be any population with size greater than some condition, n≥nc, far from extinction. For the approximation to be accurate, nc just needs to not be too close to zero. We then track both the occupation PI(t)(t) of that mean-field quantity and as well as its value I(t). In our birth–death example, I(t) represents the expected population size for populations with at least nc individuals, while PI(t)(t) represents the probability that a given population has at least nc individuals at time t. This approximation is illustrated in figure 1.

To include this mean-field limit in our master equation framework, we need to track the transitions in and out of that mean-field limit by assuming some distribution around its average state I(t) for the fraction PI(t)(t) of systems therein. Indeed, based on whatever distribution we assume around the average state, we can calculate the probability that a system in the mean-field limit is actually right on the edge of the limit at n=nc. Given that probability, we calculate the rate at which a system might leave the mean-field limit back to an explicitly tracked state at n=nc−1 in the master equation.

Let us write the system of equations for the birth–death process using the scheme shown in figure 1 where we use a mean-field I(t) to describe systems with nc=5 or more active individuals (which is an arbitrary choice for the mean-field condition nc). We need four types of general master equation states, the last master equation state at nc−1 which is coupled to the mean-field limit, the occupation probability or occupation number of the mean-field limit, as well as the position of the mean-field limit itself. In the same order, these four types of equations can be written as follows for the birth–death

The equations for all occupation numbers follow the same logic as before where every transition from state y to z (or arrows in figure 1) corresponds to a negative term for the state y and an equal positive term for the state z.

Importantly, the transition out of the mean-field limit and back to the master equation involves a coupling function ρ[nc,I(t)], where we need to approximate the distribution of states around the expected I(t) mean-field limit. In general, we recommend using a truncated Poisson distribution. It is not the right distribution in this particular birth–death example, because of the quadratic death term, but still works relatively well. Thus, ρ[nc,I(t)] is the probability that a system drawn from a Poisson distribution with mean I(t) and support [nc,∞] is found at exactly nc, on the edge of that distribution and susceptible to leaving the mean-field limit.

Note that equation (2.7) governing the position of the mean-field limit I(t) is exactly the same as the standard mean-field description used in equation (2.1). Therefore, if we define the mean-field limit as a system with nc=0 or more active individuals, we fall back on a standard mean-field model. Conversely, as we push that mean-field limit to infinity, i.e. nc→∞, we recover an exact master equation description. Between these limits, we approximate the true distribution of states, capturing stochastic extinctions around the absorbing state while maintaining a much smaller system of equations than in a true master equation framework.

We compare the predictions of the different models introduced so far in figure 2.

Coupled subsystems with approximate master equations with mean-field limits

Consider now that there are two weakly coupled subsystems (e.g. neighbouring populations) undergoing the same birth–death process. Active individuals in one subsystem give birth to new entities that migrate to the other subsystem at a rate λ.

In a master equation framework, we would then want to track the overall state of the system using y→=[n1,n2] capturing the number of active individuals in both subsystems. It would therefore be natural to extend this idea to our master equation with mean-field limit. However, the complexity of this approach obviously grows exponentially with the number of subsystems.

We rely on a simpler approximation to describe coupled subsystems, called an approximate master equation [16–18]. The idea is to describe every system separately with its own master equation, and couple them using a mean-field quantity. Using this idea on master equations with mean-field limit, we can write the following system of equations for the two subsystems birth–death process, with subsystems labelled i and j ∈ [1,2],

where we introduced an approximate mean-field coupling representing the expected number of births in i coming from j,

which is the expected migration rate based on the average state of the neighbour subsystem. The average state is calculated as a sum over the number k of individuals in the subsystem j, plus its own mean-field limit I(j)(t). In what follows, we often drop the subsystem label on mean-field limits to avoid complicating the notation (i.e. I(j)(t) would simply be I(t)).

The mean-field coupling could be improved by including known correlations using a clever moment closure in equation (2.12), such as network effects or other correlations [16,17]. Yet, taking a simple average is a good first-order approximation for most processes.

Importantly, this approach does not exactly capture joint extinctions where both subsystems have 0 active individuals. For example, the state where subsystem i goes extinct, P0(i), is no longer an absorbing state due to the average coupling θj although it does capture the vanishing internal dynamics of system i. The joint absorbing state can still be reached, but is not exactly described since we ignore correlations between the subsystems such that θi or θj do not go to zero. This is the cost of using a mean-field coupling across subsystems for the resulting model to scale linearly with the number of subsystems rather than exponentially.

We compare how well mean-FLAME captures coupled stochastic birth-death processes in figure 3.

Two-dimensional generalization with a Lotka–Volterra example

Species interactions, modelled with Lotka–Volterra systems of equations, are one of the most ubiquitous examples of mean-field models in ecology. They capture a simplified version of predator–prey or competition dynamics. But in many ways, the approximations made are extreme. Consider the following predator–prey dynamics in a population of F(t) prey fish and S(t) sharks at time t. Every prey reproduces logistically at rate μ[KF−F(t)] in an environment with carrying capacity KF (which only affects the dynamics of the prey fish, and not that of the sharks). Predators can come into contact and eat any prey at a given rate β. Predation allows the predator to have enough energy to reproduce (at rate βF) and die at a background rate, ν. Assuming that contacts between predators and prey are well mixed, we can write the following mean-field

These classic Lotka–Volterra equations yield some interesting insights about cycles in the population of predators and prey. However, populations in these cycles can go arbitrarily close to zero, which are absorbing states; the model does not allow extinctions because it ignores the fact that the dynamics are stochastic and that animals are discrete quantities [19].

In this section, we focus on a single system (no spatial aspects or coupling between subsystems) to demonstrate how to generalize the internal dynamics of mean-FLAME models to multi-state or multidimensional systems with more than one degree of freedom (or variable of interest). We do so for pedagogical reasons. The complexity of mean-FLAME models grows linearly when adding subsystems but grows exponentially when adding internal degrees of freedom to each subsystem.

Two-dimensional master equation with mean-field limit

We now generalize our mean-FLAME approach to multidimensional systems, which gets conceptually quite complicated. Readers can refer to the state and transition schema shown in figure 4. Assuming the number of fish is a vertical dimension and the number of sharks is a horizontal dimension, we can explain the different transitions as follows. Vertical state transitions correspond to prey reproduction while transitions to the left correspond to predator death. The act of predation, in the general master equation states, corresponds to the diagonal transitions that simultaneously decreases the number of prey by one and increases the number of predator by one. However, in the mean-field limit, this transition is diagonal only if the number of prey is such that the systems lies at the edge of the mean-field approximation and predation takes the prey system out of its limit. Otherwise, it is a horizontal move to the right that increases the number of predators while also leaving the prey system in the mean-field limit of prey numbers.

We now have to introduce many mean-field limits for both prey and predator numbers given all possible combinations of population sizes for prey and predator populations. There is a mean-field quantity for the number of prey given zero predators in the system, or one predator, or two predators, and so on; likewise for mean-field limits for the number of predators. Finally, there is a double mean-field limit where counts of both predator and prey are approximated by mean-field quantities.

In two dimensions, we therefore have nine types of occupation (i) the general master equation, (ii, iii) states that are next to one mean-field limit but not the other, (iv) a single state that is next to both mean-field limits, (v, vi) mean-field limit in one dimension but far from the limit in the other, (vii, viii) mean-field limit in one dimension and next to the other limit, and (ix) the double mean-field regime. This is without counting the mean-field quantities that also need to be followed.

First, the general master equation of a Lotka–Volterra predator–prey system can be written as

For large numbers of prey and a small number of predators, the system is eventually coupled to a mean-field limit. If the limits start at nF=ncF and nS=ncS, we write the following close to the limit in nF but with nS<ncS−1:

A similar equation governs the occupation number of the corresponding mean-field limit in the number of

Notice that some subtle quantities appear in these equations, such as the mean-field quantity for prey, FnS−1(t), which corresponds to the mean-field expectation for the number of prey fish in a system with nS−1 sharks at time t. That quantity in the predation term as the mean-field quantity itself mediates the rate at which predation occurs in the mean-field limits. The important distinction in structure with previous equations is that predation only adds to this equation number if the system in which predation occurs does not leave the mean-field limit (diagonal versus lateral move); hence the 1−ρ term. These last equations are then a closed system once we write out the corresponding mean-field quantities,

which are typical mean-field equations relating the change in expected number of prey FnS(t) for a system with carrying capacity KF, birth rate μ, predation rate β and a precise known number of predators nS.

The same approach can be used to describe the mean-field limits of the number of predators. The full system of equations for this case study is presented in our electronic supplementary material, appendix S1, along with all our other case studies.

Extinction probability

To illustrate the usefulness of the approach, we rely on the ability of the master equations to exactly capture stochastic dynamics, and therefore extinctions, in a Lotka–Volterra system (see figure 5 for a comparison of master equation and mean-field equations). Here, we define extinction as a system with no prey or predator left by time t, and the probability of that happening is captured by P0,0(t) in our approach. We then investigate how mean-field limits allow us to capture the dynamics of extinction with fewer equations. The results of this case study are shown in figure 6. We use systems with n explicit master equation states per species. With n=0, we have a fully mean-field description which is deterministic and thus predicts no extinctions. Interestingly, a system with n=1 actually overpredicts extinctions, most likely because the joint-Poisson approximation for the double mean-field regime underestimates correlations and overestimates the likelihood of having low numbers of both predator and prey. That being said, the system gives near perfect predictions with n as low as three. In this particular case, compared to a full master equation framework, that is 20 times fewer equations in the mathematical model without much loss of precision.

Metapopulation epidemic models

Having introduced all of our modelling techniques, we now turn our focus to the two spatially explicit complex models described in the introduction (figure 7) . For the first example, we will show how our proposed technique can be used to track the spread of an emergent pathogen through a heterogeneously distributed metapopulation where most local communities are very sparsely populated.

Susceptible-Infectious-Recovered-Susceptible (SIRS)

In this example, we want to model an epidemic process where individuals can take one of three susceptible (s), infected (i) or recovered (r). Individuals can cycle through all states and back based on a general set of mechanisms referred to as the SIRS process [20]. Susceptible individuals become infected by contact with infected individuals at a infection rate β. Infected individuals recover and gain immunity at a recovery rate α. And recovered individuals lose their immunity and become susceptible again according to a waning rate γ. In general, the epidemic can be summarized through its reproduction number R0=βN/α, which represents the expected number of infections caused by a first infected case. An epidemic would die out with R0<1, but grow exponentially when R0>1. We can then also define a herd immunity threshold generally estimated as h=1−1/R0=1−α/βN such that h is the fraction of the population that needs to be immune for the reproduction number R0=βN(1−h)/α of the infectious disease to drop below one in a population with N(1−h) susceptible individuals. If there are more susceptible individuals, the section of the population of immune individuals is below the herd immunity threshold and can cause a supercritical epidemic (defined as R0>1) until the the herd immunity threshold is reached.

Current epidemic models of this process that are amenable to simpler mathematical techniques work best in large, densely populated areas (e.g. large population centres) with well-mixed contact patterns between individuals. However, in sparsely populated areas (e.g. rural areas) these methods are less effective since the influence of small population numbers and stochastic fluctuations can have greater impact and need to be considered more explicitly.

In fact, stochastic effects can have surprising impacts in small populations which are never considered in models used for large, high density populations. For example, consider a stochastic epidemic process that starts in a small local population of 50 individuals. Even with R0>1, transmission might sputter out quickly [21] and die out as a small cluster of five cases, for example. The settlement is then left with 45 susceptible people, a 10% reduction in global susceptibility. Small populations can therefore easily reach herd immunity without ever seeing supercritical outbreaks since even small outbreaks can build immunity in a large fraction of the population. This can, in turn, affect how easily or quickly an epidemic might be able to spread across a sparsely populated area.

In a population of size N, we need to know the number of individuals in at least two out of three epidemiological states to know the complete state of the system. For example, if we know the numbers i and r of infectious and recovered individuals, we know that we have s=N−i−r susceptible individuals given that we have a closed population. We can therefore use a two-dimensional mean-FLAME approach to describe the system.

The approximate master equation for a single SIRS process without mean-field limits is given by

where θ is again the mean-field coupling estimate for one local population to other populations of the metapopulation. We then introduce two-dimensional mean-field limits for both the i and r dimensions to simplify the description of large epidemics in large populations. We can denote these Ir(t) (or Ri(t)) for mean-field numbers of infectious (or recovered) individuals in populations with exactly r recovered (or i infectious) individuals; and PIr(t),r(t) (or Pi,Ri(t)(t)) represent the occupation probabilities of these states. For populations in the double mean-field limits, we track IR(t)(t) and RI(t)(t). The schema of the model is presented in figure 4, and the resulting system of equations is again written in electronic supplementary material, appendix S1.

To account for spatially distributed metapopulations of potentially infected individuals, we embed n by n such systems in space where every local population is given a spatial position (x,y), a set of states as in figure 4, and a connectivity kernel to other local populations j∈N(x,y) which represents neighbouring populations. The set of all neighbourhoods N(x,y) could be any arbitrary metapopulation network. Over these coupled metapopulations, we can calculate the expected coupling θ(x,y) specific to the (x,y) local population as

This mean-field quantity is again a basic coupling term with a λ factor (akin to a migration rate). This factor multiplies a sum over all neighbours j of the (x,y) local population where we calculate their expected number of infected individuals by summing over the three types of general master equation (first sum in braces), infectious mean-field with exact recovered state (second sum in braces) and double mean-field state (final term in braces). To clarify the notation of the second sum and last term, recall that PIr(t),r(j)(t) and PIR(t),RI(t)(j) represent the probability of finding population j in a state in the mean-field limit of infectious individuals and, respectively, exactly r recovered individuals or also in the mean-field limit of recovered individuals RI(t). Therefore, Ir(t) and IR(t) represent the expected number of infectious individuals for a system in a mean-field limit with either exactly r recovered individuals or in the mean-field limit of recovered individuals (represented by a capital R). Likewise, RI(t) represent the expected number of recovered individuals in this double mean-field regime.

Note that this mean-field coupling does not depend on the state of the (x,y) local populations, which makes calculations easier but means that we ignore all correlations induced by the dynamics. That being said, a future version of the model could include known correlations in this coupling, since if we know that the population in (x,y) is in a state with low or high infectious counts, then neighbouring populations are likely to be in a similar state.

Spread of epidemics in heterogeneous spatially distributed populations

We evaluated how well a mean-FLAME model captured the stochastic diffusion of an SIRS epidemic process through a metapopulation system by comparing it to a classic mean-field description. We quantified the temporal patterns of stochatistic diffusion by looking at the time to introduction of the epidemic in specific local populations. The results of this case study are shown in figure 8.

As expected, a classic mean-field description overestimates the homogeneity of potential disease spread. We obtain this mean-field description by tracking the expected number of infected and recovered individuals for each local population in the metapopulation system, and then by assuming a Poisson distribution around that average. Doing so creates two big issues. First, the mean-field approach overestimates how likely it is for a given population to have had a first infected individual. In other words, even if the expected number of infected in a given local population is less than one, internal dynamics drive the expected number of infections to increase when it should not. Second, the mean-field approach assumes that the number of internal contacts between susceptible and infectious individuals is given by the product of the expected numbers of susceptible and infectious individuals, ignoring potential correlations between these two numbers. This slows down the spreading dynamics for systems in the tail of the distribution of active cases. If the number of active cases is higher than average, the number of recovered cases might be lower than average since some cases are taking longer to recover, leading to a higher total number of contacts between susceptible and infected. However, mean-FLAME does not suffer from this same issue. This allows the mean-FLAME model to produce more heterogeneous distributions of potential cases per population than classic mean-field models (see figure 8).

Of note, a full approximate master equation would have 12 321 equations (population size squared) per local population to capture the full range of potential number of infectious and recovered cases. A real master equation framework would essentially be impossible, with a total number of equations scaling exponentially with the number of populations. Yet, the mean-FLAME model used in figure 8 runs more efficiently by using only 142 equations per population by approximating high numbers with a mean-field limit.

Species dispersal model

For our final example, we approximated the range shift of a tree in response to a slowly shifting landscape of suitable environments driven by climate change. Specifically, we model a tree species that moves by seed dispersal over a discrete space (x,y) where every grid cell in that space has local conditions that specifies the cell’s local carrying capacity K(x,y,t) for the tree species at time t.

Handling of complex life cycles

We increased the realism of our model by considering a four-stage life cycle for the tree (i) seeds (s), (ii) seedlings or small greens (g), (iii) saplings or young trees (y), and (iv) adult trees (a). Each stage dies at a rate, νx, and transitions into the next stage of the cycle at rate μx, where all rates are stage dependent. For instance, adult trees produce seeds at a fast rate of μa while young trees transition into adults at a slower rate equal to μy times the square root of the number of young trees times the difference between the carrying capacity and the current population. This transition from young to adult tree stages captures the non-linear growth expected when increasing biomass leads to an increase in competition for resources among individuals [23]. The inclusion of nonlinear terms in a model increases the potential error of using deterministic mean-field models. Since the population of trees at any stage of the life cycle is unconstrained, this model has four degrees of freedom for each spatial cell.

In general, the complexity of mean-FLAME models grows exponentially with the degrees of freedom for each spatial cell because of the underlying master equation structure. However, instead of systematically moving to master equation models of higher dimensionality, we can also use purely mean-field descriptions with no explicit master equation states for degrees of freedom that are systematically far from absorbing states. For example, when modelling a tree species, the number of seeds S(t) might not need to be explicitly tracked as it tends to be high and is therefore not as sensitive to fluctuations near zero. This allows us to focus on stages of the life cycles which are population bottlenecks. In this example, the numbers of young trees Y(t) and adult trees A(t) are more important to track explicitly as they are fewer in numbers, have more influence and drive the numbers of seeds and seedlings at lower stages of the life cycles. We therefore track Y(t) and A(t) with master equations and mean-field limits, and these quantities will influence the purely mean-field descriptions for the number of seeds SY(t),A(t)(t) and seedlings or small greens GY(t),A(t)(t). We can think of this generalization as adding third and fourth dimensions to our existing model, but one where there are no explicitly tracked master equation states, only mean-field equations. The resulting state transition scheme is illustrated in figure 9 for a tree lifecycle and the full system of equations can be found in our electronic supplementary material, appendix S1.

Speed of diffusion across changing climate

In figure 10, we show an example of how a deterministic mean-field approach can overestimate the speed of diffusion of a species over a landscape. The key factor is that the mean-field approach does not capture nonlinear terms correctly. For example, sites with very low carrying capacity might spend most of their time without adult trees yet will consistently produce seeds in the mean-field version. Likewise, we used a model where young trees are competing for resources (e.g. light) such that the growth potential of the population scales with the square-root of the number of young trees rather than linearly. The mean-FLAME approach captures the fact that young trees are discrete quantities and systems will have 0 young trees, or 1, or 2, etc. In a mean-field model, the growth rate of the first young trees is artificially boosted when their expected quantity can be any continuous number above zero.

Finally, note that the results of figure 10 used four mean-field equations per grid cell in the mean-field model, and a total of 141 equations per grid cell in the mean-FLAME implementation. The latter is composed of 5×5=25 exact master equation states describing the number of young and adult trees, followed by 10 mean-field limits (two equations each) and 1 double-mean-field regime (two equations), multiplied by 3 to add mean-field limits of seeds and greens. A similar but smaller mean-FLAME implementation is represented in figure 9.

Conclusion

Standard mathematical models that rely on deterministic or mean-field approximations struggle to capture scenarios where we need to account for the stochastic and discrete nature of spreading or diffusion processes. We developed a master equation description close to absorbing states coupled with mean-field limits far from the absorbing states. This hybrid mean-FLAME modelling approach allows us to succinctly capture local extinctions mechanisms as well as stochastic bottlenecks in life cycles or spreading processes. In a nutshell, when tracking only ϵ states around an absorbing state in a d-dimensional system of m populations of size N, mean-FLAME appears to perform as well as a full approximate master equation while requiring mϵd equations instead of mNd equations. That represents a gain of (N/ϵ)d in computational complexity. Conversely, a full master equation approach is often impossible as it would require Nmd equations. This speed-up can be incredibly useful in high-dimensional systems.

Beyond the mathematical details, our main contribution is the idea that mathematical models could track individual system states when important (i.e. around absorbing states) and approximate them when possible (i.e. far from absorbing states). We implemented this idea in a hybrid mean-field and master equation framework, which means it can be directly used in applications where master equations are common but can be computationally dynamics on higher-order networks [24], heterogeneous group-structured populations [25], co-evolution of individuals and institutions [26] or noise-induced transitions [27]. The core idea could potentially be used to simplify other mathematical frameworks that track individual states, such as probability generating functions or belief propagation frameworks that track marginal distributions of potential outcomes [28,2930].

Implementation of our methods is not without challenges. First, as made clear in our case studies, one must carefully walk through multiple types of states. There are three types of states when the system has a single dimension—the general master equation, the state coupled to the mean-field limit, and the mean-field limit—but 3d types of states in d-dimensional systems. Second, numerical integrators might struggle with a single system of equations where some values like the occupation number have small values and derivatives, while others like the position of the mean-field limits can have gigantic values. Interestingly, the mean-field limits behave independently of their occupation numbers such that one can integrate the mean-field equations first and then use their solutions as an input to integrate the system of occupation numbers. As a starting point, we offer different numerical approaches in both C++ and Python as part of our software release. Some pedagogical examples are available at https://bit.ly/meanFLAME, and we have a complete repository of our different implementations at https://github.com/nwlandry/mean-FLAME.

Future work could improve on the presented hybrid approach by creating stages of mean-field limits, essentially binning non-explicit states in different ways. Going back to a birth–death process, the number of active individuals could diverge if the death rate was not quadratic and we could then imagine creating a mean-field buffer for systems with more than 10 individuals but fewer than 100 and a mean-field limit beyond that. Following a master equation description with a series of mean-field compartments with carefully calibrated coupling might help capture the heterogeneity of growth processes that do not converge to a homogeneous distribution around a steady-state (e.g. rich-get-richer processes).

Recall that the mean-FLAME approach introduced here can fall back on our standard mean-field models in one limit (mean-field limit going to the absorbing state) and on standard master equations in the other (mean-field limit going to infinity). And despite ongoing implementation challenges, the resulting modelling tool can streamline the process of quickly exploring scenarios or parameter space using a simple mean-field model, then adding precision around important stochastic events as desired by adding explicit states as needed and tuning the position of the mean-field limit. We therefore hope that this modelling approach will be used to explore possible interventions to better control complex systems around their absorbing states, especially in rural or marginal areas where system sizes can be much smaller.

Finally, we have presented use cases of the mean-FLAME approach that focused on extinction dynamics, the spread of infectious diseases and dispersal of individuals into previously unoccupied areas. However, it is worth considering other contexts that would benefit from our near absorbing state approximation of stochastic processes. Specifically, our approach might allow more explicit population viability analysis particularly for small, isolated local populations found in increasingly fragmented landscapes [31]. Moreover, in the coming decades we expect a reshuffling of most ecological systems. Modelling techniques like we have presented here that simplify the number of states to be tracked while maintaining the multidimensional complexity of the system may improve our ability to consider the consequences of new combinations of organisms, potential priority effects that may influence the resiliency of the new ecosystem, or even the potential for new pandemics [32]. More generally, we hope that the approach we have proposed helps to highlight critical gaps in how we are approaching the development of ecological models when we assume deterministic dynamics or that the expected state of a system represents its full distribution of possible states. Life at the margins of systems may in fact be more critical when modelling complex systems and processes.