Interpreting how nonlinear diffusion affects the fate of bistable populations using a discrete modelling framework

Understanding whether a population will survive or become extinct is a central question in population biology. One way of exploring this question is to study population dynamics using reaction–diffusion equations, where migration is usually represented as a linear diffusion term, and birth–death is represented with a nonlinear source term. While linear diffusion is most commonly employed to study migration, there are several limitations of this approach, such as the inability of linear diffusion-based models to predict a well-defined population front. One way to overcome this is to generalize the constant diffusivity, D, to a nonlinear diffusivity function D(C), where C>0 is the population density. While the choice of D(C) affects long-term survival or extinction of a bistable population, working solely in a continuum framework makes it difficult to understand how the choice of D(C) affects survival or extinction. We address this question by working with a discrete simulation model that is easy to interpret. This approach provides clear insight into how the choice of D(C) either encourages or suppresses population extinction relative to the classical linear diffusion model.

Interpreting how nonlinear diffusion affects the fate of bistable populations using a discrete modelling framework Yifei  Understanding whether a population will survive or become extinct is a central question in population biology. One way of exploring this question is to study population dynamics using reaction-diffusion equations, where migration is usually represented as a linear diffusion term, and birth-death is represented with a nonlinear source term. While linear diffusion is most commonly employed to study migration, there are several limitations of this approach, such as the inability of linear diffusion-based models to predict a well-defined population front. One way to overcome this is to generalize the constant diffusivity, D, to a nonlinear diffusivity function D(C), where C > 0 is the population density. While the choice of D(C) affects long-term survival or extinction of a bistable population, working solely in a continuum framework makes it difficult to understand how the choice of D(C) affects survival or extinction. We address this question by working with a discrete simulation model that is easy to interpret. This approach provides clear insight into how the choice of D(C) either encourages or suppresses population extinction relative to the classical linear diffusion model.

Introduction
Predicting whether a population will survive or go extinct is a key question in population biology [1][2][3][4][5]. For example, predicting whether a species released into a wild area will survive is crucial in protecting endangered animals [6]. Similarly, whether cancer spreads to a different body part from a primary tumour site depends on the survival of small numbers of tumour cells growing successfully in new locations [7,8]. A classical continuum model for studying the survival of biological populations is the strong Allee effect model, based on an ordinary differential equation (ODE), where C(t) ≥ 0 is the population density at time t ≥ 0, λ > 0 is the intrinsic growth rate, K > 0 is the carrying capacity density and 0 < A < K is the Allee threshold density [9][10][11][12][13][14][15]. The fate of a population described by (1.1) depends solely upon the initial density, C(0). Extinction occurs if C(0) < A, leading to C(t) → 0 as t → ∞. By contrast, the population survives if C(0) > A, leading to C(t) → K as t → ∞. Such dynamics, leading to either eventual survival or extinction, are sometimes called bistable population dynamics. The strong Allee effect model belongs to a broader class of population models, called bistable population dynamics models, and the key feature of these models is that they involve three equilibrium points: C = 0 and C = K > 0 are stable equilibrium points, and C = A, where 0 < A < K is unstable. There are many ODE models of this kind, not just the classical cubic form in (1.1) [8,[16][17][18].
To investigate spatial effects, such as moving invasion fronts, some studies consider incorporating equation (1.1) into a reaction-diffusion equation, where the population density depends upon both position and time [19][20][21][22][23][24][25][26][27][28]. In reaction-diffusion models, the dynamics of bistable populations involve a more complicated interaction between the bistable source term and the diffusion term. Unlike ODE models where the fate of a bistable population is solely determined by the initial density, many factors influence whether the population will survive or go extinct in reaction-diffusion models [20,[29][30][31][32]. For example, on an infinite domain, the initial area occupied by a bistable population needs to be greater than a threshold, called the critical initial area, so that the population avoids extinction [20].
Most reaction-diffusion models in population biology consider a constant diffusivity associated with Fick's first law of diffusion, which states that the diffusive flux is proportional to the spatial gradient of density [10,19,[21][22][23]25,33,34]. In one spatial dimension, the diffusive flux is J = −D∂C(x, t)/∂x, where D > 0 is the constant diffusivity. Linear diffusion is popular in modelling biological populations, since this model is very simple, and has a straightforward connection with a range of underlying stochastic models, such as conceptualizing the motion of individuals in the population as a simple unbiased random walk in the dilute limit, where interactions between individuals are weak [10,[34][35][36]. However, despite the immense popularity of linear diffusion, there are well-documented circumstances where population dynamics cannot be described by this simple model. For example, sharp moving fronts in cell migration assays cannot be represented by linear diffusion, and so there has been great interest in modelling the motion of well-defined fronts using degenerate nonlinear diffusion terms [37][38][39][40][41]. Similar to cell biology applications, mathematical models of insect dispersal with linear diffusion are unable to replicate observations where well-defined sharp fronts play an important role [34,42]. Therefore, reaction-diffusion equations with nonlinear diffusion are considered in a variety of applications where, in one spatial dimension, the flux is J = −D(C(x, t))∂C(x, t)/∂x, with the key difference that the constant linear diffusivity D is now generalized to a nonlinear function D(C) > 0 [34,37,[39][40][41][42][43][44][45][46][47][48][49][50]. Unlike the constant diffusivity that can be interpreted as undirected random motion of individuals without interaction [35,36], identifying the behaviour of individuals corresponding to a given nonlinear diffusion term is less clear. Therefore, it is not always obvious which nonlinear diffusion term is appropriate to model a particular situation [34,40,44,51]. Since it is known that nonlinear diffusion can impact the conditions required for survival of a bistable population [52,53], exploring the behaviour of individuals is helpful to provide a simple interpretation of how D(C) affects the fate of bistable populations subject to a nonlinear diffusion migration mechanism. Therefore, it is valuable to study the connection between the behaviour of individuals and the nonlinear diffusion term in reaction-diffusion equations. To connect continuum models with the behaviour of individuals, we work with a physically intuitive discrete framework. The discrete model incorporates straightforward crowding effects into birth, death and movement of individuals on a two-dimensional hexagonal lattice [32]. In particular, we quantify the influence of crowdedness on the motility of individuals by using a crowding function G(C) > 0, which explicitly describes how the local crowding affects the ability of individuals to move. The continuum limit of the discrete model is a reaction-diffusion equation with a strong Allee effect source term, and a general nonlinear diffusivity function D(C). This framework allows us to investigate population dynamics through repeated simulation of the discrete model, as well as solving the associated reaction-diffusion continuum limit model numerically. Through the mathematical relationship between the nonlinear diffusivity function D(C) and the underlying crowding function G(C), we develop an intuitive understanding of how different choices of D(C) affect the extinction or survival of the population. To improve our understanding, we derive expressions for the density-dependent flux of populations associated with the discrete model. The expression for the flux can be rewritten as the flux associated with a linear diffusion mechanism plus a term, which we interpret as a correction that is associated with the effects of nonlinear diffusion. Writing the flux in this way allows us to directly relate how different choices of D(C) either encourage or suppress extinction. All interpretations of our modelling framework are supported by a suite of stochastic simulations and numerical solutions of the associated continuum limit reaction-diffusion equation. All numerical algorithms required to replicate our work are available on Github.

The discrete model and the continuum limit
In this section, we introduce a lattice-based discrete model and the corresponding continuum limit model description that is closely related to our previous work [32]. Unlike the work in [32], which only considers examples where the motility of individuals is given by a linear diffusion mechanism, here we focus on a more broad range of motility mechanisms that include a range of choices of nonlinear diffusivity functions.
In the discrete model individuals are represented as agents on a two-dimensional hexagonal lattice with spacing > 0. A lattice site s, indexed by (i, j) with a unique Cartesian coordinate (x, y), is either occupied C s = 1, or vacant C s = 0. Stochastic simulations include birth, death and movement events, and we will now explain the details of these mechanisms.
If there are Q(t) agents on the lattice, we use a random sequential updating method to evolve the discrete model from time t to time t + τ . To achieve this we select Q(t) agents at random, with replacement, and give those agents an opportunity to undergo a movement event. We then select another Q(t) agents, at random, with replacement, and give those agents an opportunity to undergo a birth/death event. Once these two sets of events have been assessed, we advance time from t to t + τ and repeat until the desired output time is reached [54].
For a potential motility event, if the agent in question is at site s, that agent will move with probability M = MG(K s ), where M is the probability that an isolated agent will attempt to move during a time interval of duration τ , and G(K s ) ∈ [0, 1] is a movement crowding function which quantifies how crowding in a small neighbourhood of s influences motility. We interpret G(K s ) to be a measure of the influence of the local density upon movement since K s is a simple measure of density around site s, given by  That is, the agent will be removed out of the lattice with probability P. If a movement event occurs, the agent at site s will move into a randomly chosen vacant site in N r {s}. Therefore, the probability for the agent at site s moving to one of the vacant sites is M/(|N r |(1 − K s )). We show the movement mechanism with r = 1 leading to |N r |= 6 in figure 1b. In this particular configuration the agent at site s has two neighbouring agents, giving K s = 1/3. The probability of undergoing a movement event is M = MG(1/3). As there are four vacant sites in N 1 {s}, the probability of moving to one of the vacant sites is M/4. Note that we require G(1) = 0 as individuals have no space to move if their neighbouring sites are all occupied.
For a potential growth event, if the agent in question is at site s, that agent will undergo a growth event with probability P = PF(K s ), where P is the probability that an isolated agent will attempt to undergo a growth event during a time interval of duration τ , and F(K s ) ∈ [−1, 1] is the growth crowding function which quantifies how crowding in a small neighbourhood of s influences the propensity of agents to proliferate or die. Since the growth mechanism includes both proliferation and death as potential outcomes, we define F(K s ) ∈ [−1, 1] such that a proliferation event takes place when F(K s ) > 0 and a death event takes place when F(K s ) < 0. In the case of a proliferation event, a new daughter agent will be placed on a randomly chosen vacant site in N r (s), whereas if a death event takes place the agent at site s will be removed from the simulation. After Q(t) potential growth events have been assessed, the value of Q(t) is updated accordingly.
We show the growth mechanism with r = 1 in figure 1c. As the agent at site s has two neighbouring agents, we have K s = 1/3. Therefore, the probability of undergoing a growth event is as there are four vacant sites in N 1 {s}, the agent will place a new agent on one of the vacant sites with probability P/4. If F(1/3) < 0, the agent will be removed from the lattice with probability P.
A key feature of the discrete model is that we use G(K s ) to explicitly describe the influence of crowding effects on the movement of individuals. We provide several examples of G(K s ) and show how they influence the movement of agents on the spatial template with r = 1 in figure 2. We first consider G(K s ) = 1 − K s , which has a constant slope, as shown in figure 2a. The probability of a movement event with K s = 0, K s = 1/3 and K s = 2/3 is given in figure 2bd, respectively. With this simple movement crowding function, the probability of the agent at site s moving to one of its neighbouring vacant sites is always M/6, which is independent of the local density, K s . We then consider a concave down function, figure 2e. The probability of a movement event with K s = 0, K s = 1/3 and K s = 2/3 is given in figure 2f -h, respectively. Compared with the simplest crowding function G(K s ) = 1 − K s , the agent has a larger net movement probability for K s ∈ (0, 1). Finally, we consider a concave up function, figure 2i. The probability of a movement event with K s = 0, K s = 1/3 and K s = 2/3 is shown in figure 2j-l, respectively. In this case, the agent has a reduced probability of movement for K s ∈ (0, 1), relative to the simplest case G(K s ) = 1 − K s . The growth mechanism of agents applies a similar way of incorporating the influence of crowding effects into the discrete model [32]. We present the pseudo-code for implementing a single realization of the discrete model in the electronic supplementary material.
If we consider the spatial template with r = 1 for the movement mechanism and r ≥ 1 for the growth mechanism of agents, the continuum limit of the discrete model is with nonlinear diffusivity function and a source term   Here, D 0 > 0 is a constant in the limit that → 0 and τ → 0 with the ratio 2 /τ held constant, and λ > 0 is constant when P = O(τ ), which implies that the continuum limit is valid when P M.
Note that in the discrete model we have K as the argument of the crowding function, and that in the continuum limit the argument of the crowding function is C. This difference can be reconciled through carrying out the full details of the discrete-to-continuum averaging arguments. Full algebraic details of the intermediate steps required to derive the continuum limit are given in the electronic supplementary material. Throughout this study we work with dimensionless simulations by setting = τ = 1 in the discrete model, which leads to D 0 = M/4 and λ = P in the continuum limit. In cell biology, experimental observations imply that cell motility is reasonably well approximated by a nearest neighbour random walk whereas cell proliferation involves the disposition of daughter agents at some distance from the mother agent [54]. Therefore, throughout this work we set r = 1 for the motility mechanism and r = 4 for the proliferation mechanism, which is consistent with previous modelling [32,54]. Moreover, as we are interested in the survival and extinction of populations, we choose a growth crowding function , which leads to a cubic source term, R(C), associated with the strong Allee effect. With this choice of growth crowding function, we have F(0) = −1 indicating that isolated agents have the largest probability of dying.

Relationship between D(C) and G(C)
Based on equation (2.3), we are now in a unique position where we can specify an intuitive crowding function for the discrete model, G(C), and use the discrete-to-continuum framework to understand how this translates into a population-level nonlinear diffusivity function, D(C). This approach is very different to the more usual approach of simply specifying some   There are two ways of taking advantage of (2.3) to study population dynamics. First, for a given movement mechanism of individuals described by G(C), we can simply substitute into this expression to give the corresponding D(C). To demonstrate this first approach, we present three examples of G(C), which were examined in figure 2, and study the corresponding D(C), see figure 3a,b. In each of these three crowding functions, we always have G(0) = 1, which is reasonable since this condition implies that isolated agents are unaffected by crowding [33]. We first consider G(C) = 1 − C, which has a constant slope and leads to a constant diffusivity D(C) = D 0 . As we mentioned in §2, the probability of an agent moving to one of its neighbouring vacant sites is always M/6, which is independent of density. This is consistent with the continuum model where the standard linear diffusion mechanism means that the diffusivity is independent of density. Next, we consider the concave down crowding function G(C) = (1 − C)(1 + C/2), which has the property that G(C) > 1 − C for all C ∈ (0, 1). For this crowding function we obtain an increasing nonlinear diffusivity function D(C) > D 0 , which is reasonable since the motility of of individuals is reduced more by crowding than in the case where G(C) = 1 − C, corresponding to linear diffusion. The first approach to use (2.3) involves specifying a physically reasonable crowding function, G(C), and using the discrete-to-continuum conservation argument to understand the corresponding population-level nonlinear diffusivity function, D(C). Of course, we can also view (2.3) as allowing us to specify a particular nonlinear diffusivity function D(C), and to understand which particular crowding function is associated with that choice of D(C). One of the challenges associated with this second approach is that a given D(C) may not lead to a physically realistic crowding function, as we will now explore. One of the most standard choices of nonlinear diffusivity is a power-law diffusivity, D(C) = D 0 C m , where m is some constant exponent. This model has played a very important role in population biology models [41,55], since combining this nonlinear diffusivity term with a logistic source term gives rise to the well-studied Porous-Fisher model [34,40,51,56,57]. If we consider m ≥ 0 and G(C) ∈ [0, 1] we can solve (2.3) to give is the hypergeometric function [58]. We present three examples with m = 1, 2 and 3, given by Here, we see that each D(C) is associated with some particular crowding function, but some of the properties of these crowding functions are not as physically reasonable as those in figure 3a,b. One attractive property of the crowding function in figure 3d is that G(1) = 0 for each case, and this is reasonable since we expect motility to cease when the lattice is packed to maximum density. One less appealing feature of the crowding function in figure 3d is that each case has G(0) = 0, which means that isolated agents do not move, and this is at odds with our intuition, and experimental evidence, that crowding reduces motility [59,60]. However, while we refer to, and interpret G(C) as a crowding function, it is possible to interpret G(C) more broadly. For example, if G(C) represents the effect of intraspecific competition in population dispersal, the modelling framework in this work allows us to explain how isolated individuals die due to a lack of competition [61]. Above all, it is insightful and interesting that we are able to take canonical choices of nonlinear diffusivity function D(C), and to explore what choice of crowding function G(C) leads to those nonlinear diffusivities.

Nonlinear diffusion influences population dynamics
In this section, we quantify various population dynamics using both the discrete and associated continuum models. In all simulations, we consider an L × L domain where L = 100, and we impose periodic boundary conditions along all boundaries. Agents are initially located in a central vertical strip of width w, which may represent a species along a river [62] or a population of cells in a scratch assay [40]. For the continuum model, since the initial distribution is independent of the vertical location and evolves with periodic boundary conditions, the population density remains independent of the vertical position for all t > 0 [63]. Therefore, equation (2.2) simplifies to where C(x, t) represents the average column density of population [54,63]. We numerically solve equation (4.1) and compute which is the total population density across the whole domain. We apply the method of lines to solve (4.1) numerically and full details of the numerical method are given in the electronic supplementary material.
To quantify results from the discrete model we always consider performing V identically prepared realizations, and use this data to calculate the average occupancy of site s, where C (v) 1} is the occupancy of site s at time t in the vth identically prepared realization. As the initial occupancy is independent of the vertical position, we can denote the average column density at time t = nτ as which corresponds to C(x, t) in the continuum model, where indexes i and j, indicating the position of site s, relate to position (x, y). We also compute the total population density across the whole domain at time t = nτ as which corresponds to C(t) in the continuum model. Figure 4a,b shows results from both the discrete and continuum models with G(C) = 1 − C and D(C) = D 0 , corresponding to linear diffusion. Discrete simulations are performed with M = 1 and P = 6/1000, leading to D 0 = 1/4 and λ = 6/1000 in the continuum model. For the discrete model, we initially locate agents in the central vertical strip with width w = 10, which means that the initial condition for (4.1) is C(x, 0) = 1 for x ∈ [45,55] and C(x, 0) = 0 elsewhere. Comparing C(x, t) and C(x, t) shows that the match between the discrete and continuum results is excellent. Moreover, the density eventually reaches zero everywhere, which suggests that the population goes extinct. We then consider a larger width w = 30, and compare C(x, t) with C(x, t) at t = 800, 1600, 2400 in figure 4b. The solutions of (4.2) also match well with the averaged data from discrete simulations. In this case, the column density eventually reaches the carrying capacity everywhere, which suggests that the population survives.
Next, we consider the same initial conditions except that we use G(C) = (1 − C)(1 + C/2) in the discrete model and the increasing D(C), given by (3.2), in the continuum model. Results in figure 4c,d correspond to w = 10 and w = 30, respectively. Again, the continuum and discrete results match well. The population goes extinct with w = 10, but survives with w = 30. Comparing results in figure 4a,b with those in figure 4c,d shows that the evolution of C(x, t) is different, and this difference is due to the role of nonlinear diffusion. Traditionally, if we were working with the continuum model alone, it would be difficult to provide a physical interpretation of these differences, but in our framework we can explain these differences through our simple crowding function, G(C). We then use G(C) = (1 − C)(1 − C/2) in the discrete model and the decreasing D(C) given by (3.3) in the continuum model in figure 4e,f. In this case the continuum and discrete results reasonably match. Again, the population goes extinct with w = 10, but survives with w = 30.
Results in figure 4a-f indicate the continuum limit of our discrete model provides an accurate approximation of the stochastic population dynamics for these three choices of crowding functions. In the electronic supplementary material we show that the discrete and continuum results also match well when we consider the power-law diffusivity. A natural question that    arises when confronted with these results is the following: introducing a nonlinear diffusivity function changes the rate at which the population spreads across the domain, and we wish to understand how these differences affect the long-term survival or extinction of the bistable population. To begin to explore this question we now vary the initial width of the vertical strip show that solutions of the continuum model with these three D(C) functions lead to different critical values of w that separate long-term survival from long-term extinction. This numerical exploration demonstrates that the fate of bistable populations depends upon the choice of D(C). We will now take advantage of our discrete-to-continuum framework to interpret how different choices of D(C) either enhance or suppress population extinction.

Interpretation of how D(C) affects extinction
To understand how different choices of D(C) affect long-term survival or extinction, we now derive mathematical expressions for the average flux of agents in the discrete model in a general setting. We consider an agent at site s, at location (x, y), where the occupancy is C (x, y, t). Assuming that the density is sufficiently smooth, the densities at the neighbouring sites can be obtained by expanding C(x, y, t) in a truncated Taylor series, and where, for convenience, we denote C(x, y, t) as C and index the densities at neighbouring sites with subscripts as shown in figure 5. Note that here we are treating the density as being independent of vertical position, which is consistent with the situation in figure 4, and we will generalize this assumption later. The transition probability of an agent moving out of site s to one of its neighbouring sites s i , for i = 1, 2, 3, . . . , 6, is Similarly, the transition probability of an agent moving from site s i into site s is Therefore, combining these expressions for the transition probabilities with the geometry of the lattice in figure 4, allows us to write down an expression for the horizontal component of the net flux of agents at site s, Substituting the expressions for the transition probabilities into (5.1) and then expanding the resulting terms in a truncated Taylor series about site s gives We note that (5.2) can be written as   figure 6 show that the change of flux influences the speed at which the population spreads in space. Due to the change of spreading speed, the initial width of the vertical strip needed for a population to survive changes as well. This suggests that the nonlinear diffusivity function affects the fate of bistable populations through influencing the flux of populations. However, in figure 6, we fix the ratio of growth and movement, P/M, while P/M also influences the fate of bistable populations. Next, we are going to vary P/M and study the influence of nonlinear diffusion on the fate of a broader range of populations.

Results in
Our results so far show that the long-term survival of the population involves a complicated relationship between the width of the initial condition w, the time scale of migration M, the time scale of growth P, as well as the particular nonlinear diffusivity function D(C). We now systematically explore this relationship by taking the (w, P/M) phase space and discretizing it uniformly for w ∈ [0, 40] and P/M ∈ [1/1000, 40, 1000], as shown in figure 7. We vary P and fix M = 1 leading to λ = P and D 0 = 1/4 in the continuum model. As we are interested in the longterm outcomes of bistable populations, we solve (4.1) and calculate C(T) for a sufficiently long period of time T, so that the outcome is either C(T) = 1 or C(T) = 0. With these simulation outcomes we identify the boundaries that separate survival and extinction on the phase diagram. Overall, results in figure 7 show that large w encourages survival, and that for each value of w there is a threshold value of P/M that determines the eventual survival or extinction of populations. In the electronic supplementary material we show that the boundaries generated from the discrete model are consistent with the boundaries identified using the continuum model. The main result in figure 7 is that the curves that delineate the survival/extinction boundary depend upon the choice of D(C), and we plot three curves for the linear diffusion model, the increasing D(C) given by (3.2) and the decreasing D(C) given by (3.3). The horizontal dashed line at P/M = 0.006 highlights results shown previously in figures 4 and 6, but this phase diagram summarizes the long-term survival/extinction patterns for a much wider choice of parameters than we explored in these previous cases.

Conclusion and future work
In this work, we consider the question of population survival or extinction, with a focus on understanding how various migration mechanisms either encourage or suppress extinction. In the population biology modelling literature, the most common way to study bistable population dynamics with spatial effects is to use a reaction-diffusion model with a linear diffusion term to represent migration and a bistable source term to represent birth-death processes. While most studies employ a linear diffusion mechanism for simplicity, there are many cases where linear diffusion is inadequate. For example, mathematical models based on a linear diffusion mechanism do not predict a well-defined front that is often observed experimentally or in the field [37]. This limitation of linear diffusion is typically overcome by generalizing the constant diffusivity, D, to a nonlinear diffusivity function D(C). One of the main challenges of working with a nonlinear diffusion framework is the important question of how to choose the functional form of D(C), and there are conflicting results in the literature. For example, in the cell migration modelling literature some studies have found that using a power-law diffusivity function D(C) = D 0 C m , with m ≥ 1 can lead to a good match to experimental data [39,40,51]. One of the features of these models it that this nonlinear diffusivity is an increasing function of density. Curiously, other researchers working in precisely the same field have suggested that a decreasing nonlinear diffusivity function is appropriate, D(C) = 1/(α + C), with α > 0 [44]. This highlights the fact that choosing an appropriate nonlinear diffusivity function is not always straightforward.
In addition to understanding how to choose an appropriate nonlinear diffusivity function, a related challenge is to understand how different forms of D(C) affect the long-term survival or extinction of bistable populations. While it has been established that different choices of D(C) impact the long-term survival of populations [53], an intuitive understanding of why different choices of D(C) encourage or suppress extinction has been lacking. In this work, we address this question by working with a very simple discrete modelling framework on a two-dimensional  hexagonal lattice, where migration and birth/death events are controlled through relatively simple, easy-to-interpret crowding functions. In particular, we work with a migration crowding function G(C), which provides a very simple measure of how the ability of an individual agent to move is reduced as a function of density, C. Our discrete-to-continuum averaging arguments provide a mathematical relationship that allows us to either: (i) specify G(C) and determine the associated nonlinear diffusivity function D(C); or (ii) specify D(C) and determine the associated crowding function G(C). This new relationship allows us to explore how the averaged populationlevel flux of agents varies relative to the classical linear diffusion model for a particular crowding function, G(C). We find that choices of G(C) that increase the flux encourage population extinction relative to the classical linear diffusion model, whereas choices of G(C) that decrease the flux suppress population extinction. These results are summarized in the conceptual diagram in figure 8 showing that, for the initial conditions considered, increasing the flux of agents tends to reduce the density across the domain as the population spreads, meaning that a greater proportion of the domain has C < A, where the bistable source term acts to reduce the population and encourage extinction. There are many ways that our work can be extended. For example, all the simulation results presented here consider a very simple one-dimensional vertical strip initial condition. These results can be generalized to other initial shapes, such as circular, square or more complicated initial populations and the mathematical and computational tools presented in this work can be applied directly to this generalization. We show that nonlinear diffusion plays the same role as linear diffusion on the fate of populations when we consider a simple well-mixed initial distribution in the electronic supplementary material. There are also many ways that the current modelling framework can be extended. For example, the discrete model can be extended to consider multiple interacting subpopulations, and the same discrete-to-continuum averaging approach could be used to construct a continuum limit model, which would take the form of a system of coupled partial differential equation models [43,[64][65][66][67][68]. We leave this extension for future consideration.
Data accessibility. The code used to generate the results presented here can be found on Github. Additional material is provided in the electronic supplementary material [69].