Dimensionality affects extinction of bistable populations
Yifei Li, Stuart T. Johnston, Pascal R. Buenzli, Peter van Heijster, Matthew J. Simpson
DDimensionality affects extinction of bistable populations
Yifei Li , Stuart T. Johnston , Pascal R. Buenzli , Peter van Heijster , Matthew J. Simpson School of Mathematical Sciences, Queensland University of Technology, Brisbane, Australia. Systems Biology Laboratory, School of Mathematics and Statistics, and Department of Biomedical Engineering,Melbourne School of Engineering, University of Melbourne, Parkville, Victoria, Australia. Biometris, Wageningen University and Research, Wageningen, The Netherlands.
Abstract
The question of whether biological populations survive or are eventually driven to extinctionhas long been examined using mathematical models. In this work we study populationsurvival or extinction using a stochastic, discrete lattice-based random walk model whereindividuals undergo movement, birth and death events. The discrete model is defined on atwo-dimensional hexagonal lattice with periodic boundary conditions. A key feature of thediscrete model is that crowding effects are introduced by specifying two different crowdingfunctions that govern how local agent density influences movement events and birth/deathevents. The continuum limit description of the discrete model is a nonlinear reaction-diffusionequation, and in this work we focus on crowding functions that lead to linear diffusion anda bistable reaction term that is often associated with the strong Allee effect. Using boththe discrete and continuum modelling tools we explore the complicated relationship betweenthe long-term survival or extinction of the population and the initial spatial arrangementof the population. In particular, we study three different types of initial conditions: (i) azero-dimensional initial condition where the initial density is independent of position in thedomain; (ii) a one-dimensional initial condition where the initial density is independent ofvertical position in the domain; and, (iii) a two-dimensional initial condition. Our resultshighlight the often overlooked role of dimensionality in determining the long-term survival orextinction of bistable populations. These results are significant since we show that populationsof individuals subjected to identical stochastic rules can either lead to survival or extinctionas a result of differences in the initial spatial arrangement of the population. All softwarerequired to solve the discrete and continuum models used in this work are available on GitHub.
Keywords: population dynamics, birth, death, movement, reaction–diffusion, survival
Preprint submitted to arXiv January 20, 2021 a r X i v : . [ q - b i o . P E ] J a n . Introduction The classical logistic growth model is widely adopted in mathematical biology and math-ematical ecology (Kot 2001; Murray 2002; Edelstein-Keshet 2005). In the logistic model,small initial population densities increase over time to approach a maximum carrying-capacitydensity (Maini et al. 2004a,b). An implicit assumption in using the logistic growth modelis that any population, no matter how small, will always grow and survive. To address thislimitation, more complicated models have been developed, including models based on thestrong Allee effect (Allee and Bowen 1932; Lewis and Kareiva 1993; Stephens et al. 1999;Courchamp et al. 1999; Taylor and Hastings 2005; Courchamp et al. 2008; Arroyo-Esquiveland Hastings 2020). In the strong Allee model, initial densities greater than a threshold, calledthe
Allee threshold , grow to eventually reach the carrying capacity, whereas initial densitiesless than the Allee threshold eventually go extinct (Allee and Bowen 1932; Courchamp et al.1999; Taylor and Hastings 2005; Courchamp et al. 2008; Fadai and Simpson 2020). This kindof population dynamics, also referred to as bistable population dynamics (Kot 2001), is oftenadopted to model situations where the potential for population extinction is thought to beimportant (Saltz and Rubenstein 1995; Courchamp et al. 1999; Drake 2004; B¨ottger et al. 2015;Vortkamp et al. 2020). Bistable population dynamics are often studied using mathematicalmodels that take the form of an ordinary differential equation (ODE). In this case, the eventualextinction or survival of the population is dictated solely by whether the initial density isgreater than, or less than, the Allee threshold density. Such ODE models assume that thepopulation is well-mixed, and hence neglect spatial effects. Spatial effects, such as movinginvasion fronts, can be incorporated by considering partial differential equation (PDE) modelswhere the density of individuals depends explicitly upon position and time (Lewis and Kareiva1993; Hastings et al. 2005). A common PDE framework is to consider a reaction-diffusionPDE with a cubic bistable source term (Neufeld et al. 2017; Johnston et al. 2017).In this work we develop a mathematical modelling framework for studying bistable popu-lation dynamics on two-dimensional domains with periodic boundary conditions. Using thisframework we explore how different initial densities and different initial spatial distributionsof individuals influence the eventual survival or extinction of the population. Our modellingframework is based on a two-dimensional stochastic discrete random walk model on a hexagonallattice (Jin et al. 2016; Fadai et al. 2020). The discrete model is an exclusion process, sothat each lattice site can be occupied by no more than one agent. Individuals in the model2ndergo a birth-death process that is modulated by localised crowding effects (Jin et al. 2016;Johnston et al. 2017). The continuum limit of the discrete model leads to a two-dimensionalreaction-diffusion equation (RDE) with a bistable source term. This framework allows us toexplore discrete simulations together with solutions of the RDE. This approach is convenientbecause the discrete model is more realistic in the sense that it incorporates fluctuations, butthis benefit incurs additional computational overhead (Macfarlane et al. 2018; Chaplain et al.2020; West et al. 2016). In contrast, the continuum RDE model can be solved very efficiently,but does not provide any information about the role of stochasticity (Macfarlane et al. 2018;Chaplain et al. 2020; West et al. 2016).In all cases we study population dynamics on a square domain of side length L , withperiodic boundary conditions along all boundaries. We explore the role of dimensionality byconsidering three different initial conditions. For the two-dimensional initial condition wedistribute agents uniformly within a rectangle region of area w × w , as shown in Figure 1(a).For the one-dimensional initial condition we distribute agents uniformly within a column ofwidth w , as shown in Figure 1(b). For the zero-dimensional initial condition we distributeagents uniformly across the entire domain, as shown in Figure 1(c). Using both the discreteand continuum models we study the population dynamics associated with these three initialconditions, showing that the question of whether the population survives or goes extinct ismore subtle than simply considering whether the initial density is above or below the Alleethreshold. In fact, we observe different outcomes in terms of population extinction or survivaldepending on the dimensionality of the initial condition alone.This work is organised as follows. In Section 2 we describe the discrete individual-basedmodel, paying particular attention to incorporating realistic movement and growth mechanisms.For simplicity, we use the generic term growth to refer to the birth/death process in the discretemodel. The reason why we make a distinction between birth and death will become clear whenwe describe the modelling framework. In Section 3 we explain how to analyse the discretemodel using a mean-field assumption to arrive at an approximate continuum limit descriptionsin terms of a classical RDE. Results in Section 4 show how both the discrete and continuummodels can be used to study the role of dimensionality by considering the three different initialconditions in Figure 1. In Section 5, we systematically explore how population survival orextinction depends upon the dimensionality of the problem. All software required to solve thediscrete and continuum models used in this work are available on GitHub.3 w L L (a) two-dimensional w L L (b) one-dimensional 0 L L (c) zero-dimensional Figure 1:
Initial population distributions with different dimensions on an L × L rectangular domain. In (a), individuals are distributed uniformly in the central shaded region, of length w and width w . In (b),individuals are distributed uniformly in a vertical strip of width w , and height L . In (c), individuals aredistributed uniformly across the entire L × L domain.
2. Discrete model
We consider a lattice-based discrete model describing movement, birth and death events ina population of individuals on a hexagonal lattice, with lattice spacing ∆ >
0. Each latticesite is indexed by ( i, j ), and has a unique Cartesian coordinate,( x, y ) = i ∆ , j ∆ √ ! , if j is even, (cid:18) i + 12 (cid:19) ∆ , j ∆ √ ! , if j is odd. (1)In any single realisation of the stochastic model, a lattice site s is either occupied, C s = 1, orvacant, C s = 0. If there are Q ( t ) agents on the lattice at time t , we advance the stochasticsimulation from time t to time t + τ by randomly selecting Q ( t ) agents, one at a time, withreplacement, so that any particular agent may be selected more than once, and allowing thoseagents to attempt to move. Once the Q ( t ) potential movement events have been assessed, wethen select Q ( t ) agents at random, one at a time, with replacement, to attempt to undergo agrowth event, which could be either a birth or death event depending upon the local crowdingconditions.We now explain some features of the discrete model in terms of the schematic in Figure 2.In this initial description of the discrete model we consider nearest-neighbour movement andgrowth events only, and we will relax this assumption later. Figure 2(a) shows a potentialmovement event for an agent at site s , where all nearest-neighbour sites are vacant. In this case,4 (cid:35)(cid:86) s t s t + τ (cid:45) (cid:31) P (cid:81)(cid:96) s t + τ (cid:45) (cid:31) P (cid:85)(cid:47)(cid:86) s t s t + τ (cid:45) (cid:31) P (cid:81)(cid:96) s t + τ (cid:45) (cid:31) P (cid:74)(cid:81)(cid:112)(cid:50)(cid:75)(cid:50)(cid:77)(cid:105) (cid:84)(cid:96)(cid:81)(cid:35)(cid:28)(cid:35)(cid:66)(cid:72)(cid:66)(cid:105)(cid:118)(cid:44) (cid:31) M ( K ( (cid:75) ) s ) = MG ( K ( (cid:75) ) s ) (cid:34)(cid:66)(cid:96)(cid:105)(cid:63)(cid:102)(cid:46)(cid:50)(cid:28)(cid:105)(cid:63) (cid:84)(cid:96)(cid:81)(cid:35)(cid:28)(cid:35)(cid:66)(cid:72)(cid:66)(cid:105)(cid:118)(cid:44) (cid:31) P ( K ( (cid:59) ) s ) = P | F ( K ( (cid:59) ) s ) | (cid:85)(cid:28)(cid:86) s t s t + τ (cid:45) (cid:31) M (cid:85)(cid:43)(cid:86) s t s t + τ (cid:45) (cid:31) M Figure 2:
Movement and birth/death mechanisms.
In each lattice fragment site s is occupied and shadedin grey, and occupied neighbouring sites are shaded in blue, while vacant neighbouring sites are unshaded(white). In (a) the agent at site s moves with probability b M and moves to the target site, highlighted with agreen circle, with probability b M/
6. In (b) the agent at site s undergoes a birth event with probability b P andplaces a new agent on the target site with probability b P /
F >
0. In contrast,it dies with probability b P if F <
0. In (c) the agent moves with probability b M and moves to the target sitewith probability b M/
4. In (d) the agent undergoes a birth event with probability b P and places a new agent onthe target site with probability b P /
F >
0. In contrast, it dies with probability b P if F < the probability of attempting to move during the next time step of duration τ , is M ∈ [0 , M ≤ M . Here we notethat the two probabilities, M and ˆ M are, in general, different. This difference is a result of thelocal crowding effects. The special case in Figure 2(a) where the agent at site s is uncrowdedwe have ˆ M = M . If the attempted motility event is successful, the agent at site s moves toa randomly-chosen vacant site chosen among the set of vacant nearest-neighbour sites. Inthis case, as all six neighbour sites are vacant, the probability of moving to the target site,highlighted with a green circle, is ˆ M / s , where again allnearest-neighbour sites are vacant. Here, the probability of attempting to grow in the next timestep of duration τ is P ∈ [0 , P ≤ P . Again, the difference between P and ˆ P is caused by local crowding effects, and sincethis agent is uncrowded we have ˆ P = P . If the attempted growth event is successful, thereare two possible outcomes. First, the growth event is a birth event. In this case a daughteragent is placed at a randomly-chosen vacant site within the set of nearest-neighbour sites withprobability ˆ P . As there are six vacant neighbour sites, the probability of placing a daughter5gent at the target site, highlighted in green, is ˆ P /
6. Second, the growth event is a deathevent, and the agent is removed from the lattice, with probability ˆ P . The distinction betweenthe birth and death events is governed by the sign of the growth crowding function, F , thatwill be explained later.To illustrate how crowding effects are incorporated into the movement component of themodel, we now consider the schematic in Figure 2(c), where the agent at site s is surroundedby two agents, highlighted in purple. The probability of attempting to move is M ∈ [0 , M = M G ( K (m) s ). Here, K (m) s is a measure of the local density of site s , and G ( K (m) s ) ∈ [0 ,
1] is the movement crowdingfunction that specifies how the local density influences the probability of this agent to undergoa movement event. If this attempt is successful, as there are four vacant neighbour sites, theprobability of moving to the target site, highlighted in green, is ˆ
M / s is surrounded by two agents. The probability ofattempting to grow is P ∈ [0 , P = P | F ( K (g) s ) | . Here, K (g) s is again a measure of the local density of site s and the function F ( K (g) s ) ∈ [ − ,
1] is called the growth crowding function that specifies how the local densityinfluences the probability of this agent to undergo a growth event. If this attempt is successful,there are two possible outcomes reflected by the sign of F . If F >
0, the growth event is abirth event, and a daughter agent is placed at a randomly-chosen vacant site with probabilityˆ P . As there are four vacant neighbour sites, the probability of placing a daughter agent at thevacant target site, highlighted in green, is ˆ P /
4. Second, if
F <
0, the growth event is a deathevent, and the agent is removed from the lattice with probability ˆ P . The special case where F = 0 leads to no change.A key feature of our model is in the way that the local density about each site affectsmovement and growth events through the movement and growth crowding functions. Todescribe this we take N r { s } to denote the set of neighbouring sites around site s , where r ≥ s , so that |N r | = 3 r ( r + 1)(Jin et al. 2016; Fadai et al. 2020). The probability that any potential movement or growthevent is successful depends upon the crowdedness of the local region surrounding site s . We6ount neighbouring agents in N r , and consider K s ( r ) = 1 |N r | X s ∈N r { s } C s ∈ [0 , , (2)as a simple measure of the crowdedness of the local region surrounding site s . While in Figure 2we explain the model with r = 1 and |N | = 6, it is possible to use different-sized templates,depending on the choice of r . Sometimes it is useful for us to use different-sized templatesfor the movement and growth mechanisms. For example, Simpson et al. (2010) argue thatcell movement can be modelled using a nearest-neighbour random walk with r = 1, whereascell proliferation often involves non nearest-neighbour interactions since daughter cells areoften deposited several cell diameters away from the location of the mother cell. To simulatethis, Simpson et al. (2010) introduce proliferation mechanisms where daughter agents areplaced up to four lattice sites away from the mother agent to faithfully capture this biologicaldetail into their model. This would be similar to setting r = 1 for movement and r = 4 forgrowth in our model. It is thus convenient for us to make a notational distinction between thesize of the templates for motility and growth. Therefore, we denote the motility template as K (m) s = K s ( r ) and the growth template as K (g) s = K s ( r ) where r ≥ r ≥ s is measuredusing a nearest-neighbour template with r = 1, and the growth crowding function is givenin Figure 3(d) such that the probability of undergoing a birth event is ˆ P = P | F ( K (g) s ) | . InFigure 3(a) where K (g) s = 0, we have F (0) = 1 and ˆ P = P . As there are six vacant sites in N ,the probability of placing a daughter agent at the target site, highlighted in green, is ˆ P /
6. InFigure 3(b), where the agent at site s is surrounded by two neighbour agents, the probabilityof undergoing a birth event is ˆ P = 2 P/
3, since K (g) s = 1 / F (1 /
3) = 2 /
3. As there arefour vacant sites in N , the probability of placing a daughter agent at the target site is ˆ P / P = P/ K (g) s = 2 / F (2 /
3) = 1 /
3. As there aretwo vacant sites in N , the probability of placing a daughter agent at the target site is ˆ P / s using a larger spatial template with r = 2. Therefore, if7 = r = r = r = (cid:58)(cid:96)(cid:81)(cid:114)(cid:105)(cid:63) (cid:43)(cid:96)(cid:81)(cid:114)(cid:47)(cid:66)(cid:77)(cid:59) (cid:55)(cid:109)(cid:77)(cid:43)(cid:105)(cid:66)(cid:81)(cid:77) s F (1 /
3) = 2 / , (cid:31) P = 2 P/ (cid:85)(cid:35)(cid:86) s F (2 /
3) = 1 / , (cid:31) P = P/ (cid:85)(cid:43)(cid:86) s F (1 /
3) = 2 / , (cid:31) P = 2 P/ (cid:85)(cid:55)(cid:86) s F (2 /
3) = 1 / , (cid:31) P = P/ (cid:85)(cid:59)(cid:86) s F (1 /
3) = − / , (cid:31) P = 2 P/ (cid:85)(cid:68)(cid:86) s F (2 /
3) = 1 / , (cid:31) P = P/ (cid:85)(cid:70)(cid:86) s F (1 /
3) = − / , (cid:31) P = 2 P/ (cid:85)(cid:77)(cid:86) s F (2 /
3) = 1 / , (cid:31) P = P/ (cid:85)(cid:81)(cid:86) (cid:82) / / (cid:82)(cid:121) F ( K ( (cid:59) ) s ) = 1 − K ( (cid:59) ) s K ( (cid:59) ) s F ( K ( (cid:59) ) s ) (cid:85)(cid:47)(cid:86) (cid:82) / / (cid:82)(cid:121) F ( K ( (cid:59) ) s ) = 1 − K ( (cid:59) ) s K ( (cid:59) ) s F ( K ( (cid:59) ) s ) (cid:85)(cid:63)(cid:86) (cid:82) / / (cid:121)(cid:88)(cid:107)(cid:64)(cid:82)(cid:121) K ( (cid:59) ) s F ( K ( (cid:59) ) s ) = 2(1 − K ( (cid:59) ) s )( K ( (cid:59) ) s − / F ( K ( (cid:59) ) s ) (cid:85)(cid:84)(cid:86) (cid:82) / / (cid:121)(cid:88)(cid:107)(cid:64)(cid:82)(cid:121) K ( (cid:59) ) s F ( K ( (cid:59) ) s ) = 2(1 − K ( (cid:59) ) s )( K ( (cid:59) ) s − / F ( K ( (cid:59) ) s ) (cid:85)(cid:72)(cid:86) K ( (cid:59) ) s = 0 K ( (cid:59) ) s = 1 / K ( (cid:59) ) s = 2 / s F (0) = − , (cid:31) P = P (cid:85)(cid:66)(cid:86) s F (0) = − , (cid:31) P = P (cid:85)(cid:75)(cid:86) s F (0) = 1 , (cid:31) P = P (cid:85)(cid:50)(cid:86) s F (0) = 1 , (cid:31) P = P (cid:85)(cid:28)(cid:86) Figure 3:
Growth mechanisms with different-sized spatial templates and growth crowding functions . In each lattice fragment site s is shaded grey,occupied sites within the template are shaded blue, and vacant sites within the template are unshaded (white). Each subfigure shows a potential outcome for anagent at site s . The crowdedness of N is shown in (a)–(c) and (i)–(k). The crowdedness of N is shown in (e)–(g) and (m)–(o). The agent at site s can undergo abirth event when F > s can undergo a death event when F < he agent at s undergoes a successful birth event, the daughter agent is able to be placed atany vacant site within N . The probability of undergoing a birth event is ˆ P = P | F ( K (g) s ) | . Forthe agent in Figure 3(e) where K (g) s = 0 and F (0) = 1, we have ˆ P = P . In this configurationthere are 18 vacant sites in N and the probability of placing a daughter agent at the targetsite, highlighted in green, is ˆ P /
18. In Figure 3(f), where the agent at site s is surrounded bysix neighbour agents, the probability of undergoing a birth event is ˆ P = 2 P/
3, as K (g) s = 1 / F (1 /
3) = 2 /
3. Since there are 12 vacant sites in N , the probability of placing a daughteragent at the target site is ˆ P /
12. Similarly, in Figure 3(g), we have ˆ P = P/
3, as K (g) s = 2 / F (2 /
3) = 1 /
3. The probability of placing a daughter agent at the target site is ˆ
P /
6. Allresults in Figures 3(a)–(c) and Figures 3(e)–(g) consider the simplest linear crowding function F ( K (g) s ) = 1 − K (g) s , which does not allow for agents to die.We now choose a nonlinear growth crowding function F ( K (g) s ) that can take on bothpositive and negative values, as shown in Figure 3(l). In this case we make a distinctionbetween a birth event when F ( K (g) s ) >
0, a death event when F ( K (g) s ) <
0, and no event when F ( K (g) s ) = 0. We first consider a nearest-neighbour template with r = 1 in Figures 3(i)–(k).In Figure 3(i), the agent at site s dies with probability ˆ P = P | F ( K (g) s ) | . Here, K (g) s = 0 and F (0) = −
1, thus ˆ P = P . In Figure 3(j) the agent at site s dies with probability ˆ P = 2 P/ K (g) s = 1 / F (1 /
3) = − /
9. In Figure 3(k) the agent at site s undergoes a birth eventwith probability ˆ P = P/ K (g) s = 2 / F (2 /
3) = 1 /
9. As there are two vacant sites in N , the probability of placing a daughter agent at the target site is ˆ P / N in Figures 3(m)–(o). In Figure 3(m), theagent at site s dies with probability ˆ P = P | F ( K (g) s ) | . Here K (g) s = 0 and F (0) = −
1, thusˆ P = P . In Figure 3(n), the agent at site s dies with probability ˆ P = 2 P/ K (g) s = 1 / F (1 /
3) = − /
9. In Figure 3(o) the agent at site s undergoes a birth event with probabilityˆ P = P/ K (g) s = 2 / F (2 /
3) = 1 /
9. As there are six vacant sites in N , the probabilityof placing a daughter agent at the target site is ˆ P / G ( K (m) s ), is incorporated into the model in a similarway as the growth crowding function except that it is always positive, G ( K (m) s ) ∈ [0 , r = 1 for movement and r = 4 for growth. This choice is motivated by experimental images of9ell proliferation where careful examination of timelapse movies show that daughter cells areoften placed some distance from the mother cell (Simpson et al. 2010). Other choices of r canbe implemented using the software available on GitHub.
3. Continuum limit
In this section we derive the mean-field continuum limit of the discrete model. The averaged occupancy of site s , constructed from V identically-prepared realisations of the discrete model,can be written as ¯ C s = 1 V V X v =1 C ( v ) s ( t ) , (3)where C ( v ) s ( t ) ∈ { , } is the binary occupancy of site s at time t in the v th identically-preparedrealisation of the discrete model. We note that ¯ C s ∈ [0 , t , but wesuppress this dependence for notational convenience. Similarly, the averaged occupancy of N r { s } , again constructed from V identically-prepared realisations, is given by¯ K s ( r ) = 1 |N r | X s ∈N r { s } ¯ C s . (4)As we use a nearest-neighbour template, r = 1, for movement, and a larger template, r = 4,for growth, we denote the averaged occupancy of sites for potential movement events as ¯ K (m) s ,and the averaged occupancy of sites for potential growth events as ¯ K (g) s .To arrive at an approximate continuum limit description, we start by writing down anexpression for the expected change in occupancy of site s during the time interval from t to t + τ , δ (¯ C s ) = movement events into s z }| { M |N | (1 − ¯ C s ) X s ∈N { s } ¯ C s G (¯ K (m) s )1 − ¯ K (m) s − movement events out of s z }| { M ¯ C s G (¯ K (m) s )+ P |N | (1 − ¯ C s ) X s ∈N { s } H ( F (¯ K (g) s ))¯ C s F (¯ K (g) s )1 − ¯ K (g) s | {z } birth events: place new agents onto s − (1 − H ( F (¯ K (g) s )) P ¯ C s F (¯ K (g) s ) | {z } death events: remove agent from s , (5)where H is the Heaviside step function. Each term in Equation (5) has a relatively simple10hysical interpretation. The first term on the right hand side of Equation (5) representsthe change in occupancy of site s owing to the expected movement of agents in N { s } intosite s . The factor 1 / (1 − ¯ K (m) s ) accounts for the choice of the target site in N being randomlyselected from the available vacant sites. The second term on the right hand side of Equation (5)represents the change in occupancy of site s owing to the expected movement of agents outof site s . The third term on the right hand side of Equation (5) represents the change inoccupancy owing to the expected birth events of agents in N { s } that would place daughteragents onto site s , where F (¯ K (g) s ) >
0. Again, the factor 1 / (1 − ¯ K (g) s ) accounts for the choiceof the target site in N being randomly selected from the available vacant sites. The last termon the right hand side of Equation (5) represents the expected change in occupancy owing toagent death at site s , when F (¯ K (g) s ) <
0. Note that this approximate conservation statement(5) makes use of the mean-field assumption, whereby the occupancy status of lattice sites aretaken to be independent (Baker and Simpson 2010).To derive the continuum limit we replace ¯ C s with a continuous function, C ( x, y, t ), andexpand each term in Equation (5) in a truncated Taylor series about site s , and truncate termsof O (∆ ). Subsequently, we divide both sides of the resulting expression by τ and evaluatethe resulting expressions in the limit ∆ → τ → /τ heldconstant (Hughes 1995). This leads to the following nonlinear RDE, ∂C ( x, y, t ) ∂t = D ∇ · ( D ( C ) ∇ C ( x, y, t )) + λC ( x, y, t ) F ( C ) , (6)where D ( C ) = C d G ( C )d C + 1 + C − C G ( C ) , (7)and D = M ∆ ,τ → ∆ τ , λ = lim τ → Pτ . (8)Here, D is the free-agent diffusivity, D ( C ) is a nonlinear diffusivity function that relates tothe movement crowding function G ( C ), and λ is the rate coefficient associated with the sourceterm that is related to the growth crowding function F ( C ). To obtain a well-defined continuumlimit we require that P = O ( τ ) (Simpson et al. 2010). The algebraic details required to arriveat the continuum limit are outlined in Appendix B.For all simulations in this work we use ∆ = τ = 1, giving D = M/ λ = P . This is11quivalent to working in a non-dimensional framework (Simpson et al. 2010). If the modelis to be applied to a particular dimensional problem, then ∆ and τ can be re-scaled usingappropriate length and time scales. In this non-dimensional framework with τ = 1, we satisfythe requirement that P = O ( τ ) by ensuring P/M (cid:28)
1. The main focus of this work is on therole of the growth mechanism, and the question of whether the population survives or goesextinct. We therefore set G ( C ) = 1 − C , leading to D ( C ) = 1. This means that the nonlineardiffusion term in Equation (6) simplifies to a linear diffusion term, giving ∂C ( x, y, t ) ∂t = D ∇ C ( x, y, t ) + λC ( x, y, t ) F ( C ) . (9)We note that Equation (9) has been studied extensively in applications involving the spatialspread of invasive species, such as the work of Fisher (1937); Skellam (1951); Fife (1979); Lewisand Kareiva (1993) and Hastings et al. (2005). Some previous models consider a logistic-typesource term (Fisher 1937), while others consider Allee-type bistable source term (Sewalt et al.2016). Under these conditions many results have been established. For example, if we considerEquation (9) on a one-dimensional infinite domain, it is well known that this model supportstravelling wave solutions for both logistic (Fisher 1937) and bistable (Fife 1979) source terms.In this work, however, we take a different perspective by studying Equation (9) on a finitedomain and so the question of analysing travelling wave solutions is not our focus.In the rest of this work we choose F ( C ) = a (1 − C )( C − A ) , (10)since this leads to the canonical cubic source term λCF ( C ) associated with Allee kinetics. Inparticular, we set A = 2 / F ( C ) can be used to represent birth eventswhere C > / C < /
5. We further set a = 5 / F (0) = − C = 0, are always successful.In summary, our discrete model requires the specification of two crowding functions: G ( C )and F ( C ). These crowding functions are related to macroscopic quantities in the associatedRDE model. In particular, G ( C ) is related to a nonlinear diffusivity function, D ( C ), and F ( C )is related to a nonlinear source term λCF ( C ). Figure 4 shows the relationship between thesefunctions for our choice of G ( C ) and F ( C ). 12 ( C ) C (a) D ( C ) C (b) F ( C ) C . − (c) λCF ( C ) C . − . (d) Figure 4:
Specific crowding functions used in this work. (a)–(b) Setting G ( C ) = 1 − C for the movementcrowding function leads to linear diffusion, D ( C ) = 1. (c)–(d) Setting F ( C ) = 2 . − C )( C − .
4) for thegrowth crowding function leads to λCF ( C ) = 2 . C (1 − C )( C − . A = 0 .
4. Initial conditions and simulation data
In this section we consider the three initial conditions shown in Figure 1, and we introducethe corresponding continuous descriptions. In general, each of the initial conditions shown inFigure 1 can be written as C ( x, y,
0) = B, ( x, y ) ∈ H , , elsewhere , (11)where H is the region in which individuals are uniformly distributed at density B ∈ (0 , H .For the three initial conditions in Figure 1 we will report data from the stochastic modelin the following way. For simulations relating to the two-dimensional initial conditions as inFigure 1(a) we denote the averaged occupancy of site s as h C ( x, y, t ) i = 1 V V X v =1 C ( v ) ( i, j, n ) , (12)13here we note that the average denoted by the angular parenthesis is taken in the same way asthe average in Equation (3). Here, site s , indexed by i and j , are related to position, ( x, y ) viaEquation (1). Here, h C ( x, y, t ) i is a measure of the local density at location ( x, y ), and time t = nτ after the n th time step in the stochastic discrete model. The average is constructedusing V identically-prepared realisations of the stochastic model. Similarly, for simulationsrelating to the one-dimensional initial conditions as in Figure 1(b), where the initial occupancyis independent of the vertical position, we denote the averaged occupancy of any site as h C ( x, t ) i = 1 V J V X v =1 J X j =1 C ( v ) ( i, j, n ) , (13)which is a measure of the density at location x and at time t = nτ . Note that, as we willshow through simulation, the density of agents remains independent of the vertical positionfor all t > h C ( t ) i = 1 V IJ V X v =1 J X j =1 I X i =1 C ( v ) ( i, j, n ) , (14)which is a measure of the total population density at time t = nτ . As we will show throughsimulation, in this case the density of agents remains independent of position for all t > h C ( t ) i is also useful to describe simulations starting from the two-dimensional and one-dimensional initial conditions. In summary, data from the discrete modelsincludes h C ( x, y, t ) i , h C ( x, t ) i , and h C ( t ) i .To compare averaged data from the discrete model with the solution of the continuummodel we solve Equation (9) numerically to give C ( x, y, t ). Full details of the numericalmethods are presented in Appendix C. Using the numerical solution for C ( x, y, t ) we calculate C ( t ) = 1 L Z L Z L C ( x, y, t ) d x d y, which is the total number of agents divided by the total area of the L × L domain.In Figure 5, we compare data from the discrete model with numerical solutions of thecontinuum model for the two-dimensional initial condition. The initial condition in Figure 5(a)shows a square region, of size 40 ×
40 that is occupied with density B = 1. Simulations are14 a) T = 0 100100 y x (b) xy T = 0 . . h C ( x, y, T ) i (c) xy T = 0 . . C ( x, y, T )0 10001 xC ( x, , T ) , h C ( x, , T ) i (d) T C ( T ) , h C ( T ) i (e)(f) T = 0 100100 y x (g) xy T = 0 . . h C ( x, y, t ) i (h) xy T = 0 . . C ( x, y, t )0 10001 xC ( x, , T ) , h C ( x, , T ) i (i) T C ( T ) , h C ( T ) i (j) Figure 5:
Comparison of data from the discrete model with the solution of the continuum modelfor the two-dimensional initial condition.
In (a) agents are initially located in a square region of size40 ×
40 with B = 1. (b) h C ( x, y, T ) i at T = λt = 0 .
5. (c) C ( x, y, T ) at T = λt = 0 .
5. (d) h C ( x, , T ) i (dashed) and C ( x, , T ) (solid) at T = 0 . , . , .
6. (e) h C ( T ) i (dashed red) and C ( T ) (solid grey). In (f)agents are initially located at a square region of size 80 ×
80 with B = 1. (g) h C ( x, y, T ) i at T = λt = 0 . C ( x, y, T ) at time T = λt = 0 .
5. (i) h C ( x, , T ) i (dashed) and C ( x, , T ) (solid) at T = 0 . , . , .
6. (j) h C ( T ) i (dashed red) and C ( T ) (solid grey). The dashed black horizontal lines in (d), (e), (i) and (j) indicate theAllee threshold, A = 0 .
4. Arrows in (d) and (i) show the direction of increasing time. Note that we generate 40identically-prepared realisations to obtain h C ( x, y, T ) i in (b) and (g), 40 identically-prepared realisations toobtain h C ( T ) i in (e) and (j), and 4000 identically-prepared realisations to obtain h C ( x, , T ) i in (d) and (i).All discrete simulations correspond to M = 1 and P = 1 / M = 1 and P = 1 / D = 1 / λ = 1 / T = λt = 0 . T = 0 .
5. The visual comparison between the spatial arrangementand density of the profiles in Figure 5(b) and Figure 5(c) appears to match well. To makea more quantitative comparison we examine the density along the horizontal dashed linesshown in Figures 5(b)–(c) at y = 50. Figure 5(d) compares the evolution of C ( x, , T ) and h C ( x, , T ) i , and we see that the match between the solution of the continuum model andappropriately averaged data from the discrete model is excellent. Finally, in Figure 5(e) wecompare the averaged total occupancy from the discrete model, h C ( T ) i , with C ( T ) from thesolution of the continuum model. Again, we see that the discrete-continuum comparison isexcellent, and that the continuum model predicts the eventual extinction of this population.This is an interesting result given that the initial density in the centre of the domain is greaterthan the Allee threshold, yet the total population eventually becomes extinct.We now consider a second set of discrete-continuum comparisons for precisely the samemechanisms except that the spatial arrangement of the initial condition, shown in Figure 5(f) islarger and occupies the central 80 ×
80 region of the domain. As before, the match between thediscrete averaged data and numerical solutions of Equation (9) is excellent in Figures 5(g)–(j).In this case the population eventually grows to reach the maximum density.Overall, the results in Figure 5 confirm that the numerical solution of the continuum model,Equation (9), can provide an accurate way of studying the expected behaviour of the discretemodel. We are interested in the long-term outcomes of precisely the same discrete mechanismin Figures 5(a)–(e) and Figures 5(f)–(j), and we see that the population eventually becomesextinct in the former case, while surviving in the latter case. The only difference is in thespatial arrangement of the initial condition.For the one-dimensional initial condition, as shown in Figure 1(b), Equation (9) simplifiesto ∂C ( x, t ) ∂t = D ∂ C ( x, t ) ∂x + λC ( x, t ) F ( C ) , (15)where C ( x, t ) represents the column-averaged density of agents (Simpson et al. 2010). Anextensive discussion and exploration of the implications of simplifying the two-dimensionalnonlinear RDE into this simpler one-dimensional RDE is given in Simpson (2009). Given a16umerical solution of Equation (15), as outlined in Appendix C, we compute C ( t ) = 1 L Z L C ( x, t ) d x, which, again, is the total number of agents divided by the total area of the L × L domain.Results in Figure 6 give a comparison between the discrete and continuum solutions forthe one-dimensional initial condition. The initial condition in Figure 6(a) shows that thecentral strip of width 16 is occupied with density B = 1. Figures 6(b)–(c) show snapshotsfrom the discrete model as the population spreads into the domain. Figure 6(d) compares thenumerical solution of Equation (15), C ( x, T ), with averaged data from the discrete model, h C ( x, T ) i . The evolution of the total population density in the discrete model, h C ( T ) i , and inthe continuum model, C ( T ), is compared in Figure 6(e). In all cases the continuum modelaccurately captures the averaged data from the discrete model, and in this case the populationeventually becomes extinct. Again, this is an interesting result because the initial density ofagents in the central vertical strip exceeds the Allee threshold.We then consider a second set of discrete-continuum comparisons for precisely the samemechanisms except that the spatial arrangement of the initial condition, shown in Figure 6(f),is wider and occupies the central vertical strip of width 64 with density B = 1. Figures 6(g)–(h) show discrete snapshots as the population spreads. The comparisons between C ( x, T )and h C ( x, T ) i in Figure 6(i), and between C ( T ) and h C ( T ) i in Figure 6(j) are excellent. Inthis case we see that the population eventually grows to reach the maximum density. Forthe one–dimensional initial condition the same discrete mechanism again leads to differentlong-term outcomes in Figures 6(a)–(e) and Figures 6(f)–(j), where the population eventuallybecomes extinct in the former case, while surviving in the latter case. The only difference is inthe width of the initial population 17 T = 0 100100 y x (a) T = 0 . y x (b) T = 0 . y x (c) xC ( x, T ) , h C ( x, T ) i (d) T C ( T ) , h C ( T ) i (e) T = 0 100100 y x (f) T = 0 . y x (g) T = 0 . y x (h) xC ( x, T ) , h C ( x, T ) i (i) T C ( T ) , h C ( T ) i (j) Figure 6:
Comparison of data from the discrete model with the solution of the continuum modelfor the one-dimensional initial condition.
In (a), agents are initially placed within a vertical strip where x ∈ [42 , B = 1. (b)–(c) snapshots from the discrete model at T = 0 . T = 0 .
2, respectively. (d) h C ( x, T ) i (dashed) and C ( x, T ) (solid) at time T = 0 . , . , .
6. (e) h C ( t ) i (dashed red) and C ( t ) (solid grey).In (f), agents are initially placed within a vertical strip where x ∈ [18 , B = 1. (g)–(h) snapshotsfrom the discrete model at T = 0 . T = 0 .
2, respectively. (i) h C ( x, T ) i (dashed) and C ( x, T ) (solid) attime T = 0 . , . , .
6. (j) h C ( t ) i (red dashed) and C ( t ) (grey solid). The dashed black horizontal lines in (d),(e), (i) and (j) indicate the Allee threshold, A = 0 .
4. Arrows in (d) and (i) show the direction of increasingtime. Note that we generate 40 identically-prepared realisations to obtain h C ( x, T ) i in (d) and (i), and also 40identically-prepared realisations to obtain h C ( T ) i in (e) and (j). a) T = 0 100100 y x (b) T = 1 100100 y x (c) T = 2 100100 y x TC ( T ) , h C ( T ) i (d)(e) T = 0 100100 y x (f) T = 1 100100 y x (g) T = 2 100100 y x TC ( T ) , h C ( T ) i (h) Figure 7:
Comparison of data from the discrete model with the solution of the continuum model forthe zero-dimensional initial condition. (a)–(c) snapshots of discrete simulations at time T = λt = 0 , , T = 0 a fixed number of agents are randomly distributed on the lattice so that h C (0) i = 0 .
16. (d) h C ( T ) i (dashed red) and C ( T ) (solid grey). (e)–(g) snapshots of discrete simulations at time T = 0 , , h C (0) i = 0 .
64. (h) h C ( T ) i (dashed red) and C ( T ) (solid grey). The dashed black horizontal lines in (d) and (h)are the Allee threshold, A = 0 .
4. Note that we generate 40 identically-prepared realisations to obtain h C ( T ) i in(d) and (h). For the zero-dimensional initial condition, as shown in Figure 1(c), Equation (9) simplifiesto d C ( t )d t = λC ( t ) F ( C ) , (16)where C ( t ) represents the total density of agents (Simpson et al. 2010). This separable ODEcan be solved to give an implicit solution. Results in Figure 7 compare the discrete andcontinuum solutions for the zero-dimensional initial condition. In Figure 7(a), a fixed numberof agents are randomly distributed in the entire domain at T = 0, leading to h C (0) i = 0 . h C (0) i = 0 .
64 in Figure 7(e). Figures 7(f)–(g) again show discrete snapshots as thepopulation evolves, and we observe that C ( T ) approximates h C ( T ) i well in Figure 7(h). Inthis case the population survives and grows to reach the maximum density.Overall, the results in Figures 5–7 confirm that the numerical solution of the continuummodel provides a useful way of accurately studying the expected behaviour of the discrete19odel. Of interest is that the long-term fate of populations vary with the spatial arrangementof the initial conditions. Our aim now is to study these differences more carefully.
5. Role of dimensionality
Our results in Section 4 indicate that several factors are at play when we consider thelong-term fate of bistable populations. First, the initial number of individuals in the populationplays a role. Second, the spatial arrangement of the initial population also plays an importantrole. Since the initial distribution of the population is given by Equation (11), there are twodifferent ways to vary the initial population size. Either we can adjust the size or shape of theinitially occupied region H , or we can vary B to alter the initial number of individuals per unitarea within H . In the remainder of this section we investigate all these options systematically.Results in Figure 8 summarise the long-term outcome of a range of scenarios with thezero-dimensional initial condition. In this case H corresponds to the entire L × L domainand C (0) = B . We vary the initial condition by systematically varying B , as indicated inFigures 8(a)–(c), and vary the ratio P/M by holding M = 1 and varying P ∈ [1 / , / P/M = 4 λ/D , we hold D = 1 / λ in the continuous model.To systematically study the transition between population extinction to population survival,we take the ( C (0) , P/M ) space and discretise it uniformly into a square mesh, with 51 × C (0) and P/M considered, we generate 40 identically-preparedrealisations of the discrete model and we compute the survival probability, S ∈ [0 , T , which we take to be T = max(30 /P, ). Figure 8(d) summarises the outcomes of thesimulations in terms of a phase diagram. In this case the survival outcome for the continuummodel is a simple vertical line at C (0) = A . In general we see good agreement between theprediction of survival or extinction between the continuum and discrete models. There issome small discrepancy as P/M increases. This difference is consistent with the fact that thecontinuum model
P/M has to be sufficiently small otherwise the mean-field approximationis invalid and the solution of the continuum model does not necessarily provide an accuratedescription of the discrete mechanism (Baker and Simpson 2010; Simpson et al. 2010). In20 = 0 . L L (a) B = 0 . L L (b) B = 0 . L L (c)(d) . . . . . . C (0) P/M . . . . S Figure 8:
Phase diagram for extinction/survival with the zero-dimensional initial condition . (a)–(c)show how we vary the initial density with C (0) = B for this initial condition. (d) Phase diagram of a squaremesh with 51 ×
40 nodes for for C (0) ∈ [1 / , /
10] and
P/M ∈ [1 / , / M = 1. Vertical red lineindicates the survival/extinction threshold from the continuum model and the blue shading shows the survivalprobability S measured by 40 identically-prepared realisations. summary, for the zero-dimensional initial condition the long-term population survival dependssimply upon whether the initial density is above or below the Allee threshold, as expected.We now explore how the simple outcome for the zero-dimensional initial condition becomesmore complicated when we consider different initial spatial arrangements of the population.For the one-dimensional initial condition we first fix B = 1 and vary the size of H by changingthe width of the vertical strip, w . Varying the width of the strip leads to a change in theinitial density across the entire domain, C (0) = Bw /L . For example, Figures 9(a)–(c) showsthree different one-dimensional initial conditions with B = 1 and different widths. For theseinitial conditions we vary the ratio P/M = 4 λ/D by holding M = 1 and varying P in thediscrete model, and by holding D = 1 / λ in the continuous model. Again, this21 w = 20 LL (a) LL w = 40 (b) LL w = 60 (c)(d) . . . . . . C (0) P /M . . . . S Figure 9:
Phase diagram for extinction/survival with the one-dimensional initial condition . In(a)–(c) we show three different initial conditions where C (0) = Bw /L , and we vary w . (d) Phase diagram ofa square mesh with 51 ×
40 nodes for C (0) ∈ [1 / , /
10] and
P/M ∈ [1 / , / M = 1. The blackcurve indicates the survival/extinction threshold from the continuum model and the blue shading shows thesurvival probability S measured by 40 identically-prepared realisations. Vertical red line is C (0) = 0 . A = 0 . C (0) , P/M ), which we discretise into a square mesh with 51 × S , dependsupon C (0) and P/M . The boundary that separates the eventual survival and extinction inthe continuum model is shown in solid black, and the survival probability from the discretesimulations is shown in blue shading. Overall, the long-term predictions in terms of survivalor extinction are consistent between the continuum and discrete models. For completeness wealso show the red vertical line indicating the Allee threshold.It is interesting to compare the results in Figure 8(d) and Figure 9(d). For the zero-dimensional initial condition long-term survival is obtained when the initial density is greaterthan the Allee threshold. In contrast, for the one-dimensional initial condition we observesurvival even when C (0) < A . Further, in the one-dimensional case we see that the long-termsurvival is strongly dependent upon P/M whereas in the zero-dimensional case this dependenceis less pronounced. In Figure 9 we alter C (0) by fixing B = 1 and varying w . Anotherapproach is to hold w constant and vary B . Additional results in Appendix D show thatthese two approaches to varying C (0) lead to very similar outcomes.For the two-dimensional initial condition, we fix B = 1 and vary the size of H by changing w where w = w = w as shown in Figure 10(a). Varying w allows us to vary the initialdensity across the entire domain, C (0) = Bw /L . Again, we construct a phase diagram inFigure 10(d) that summarises the long-term survival outcome as a function of C (0) and P/M ,by discretising the ( C (0) , P/M ) space using a square mesh with 51 ×
40 nodes. The phasediagram in Figure 10(d) is very similar to the phase diagram in Figure 9(d). We see thatthe long-term survival strongly depends upon
P/M , and the distinction between survival andextinction predicted by the continuum limit model is a good approximation of the discretesimulation data. 23 L (a) w = 32 L L w = 60 (b) L L w = 78 (c)(d) . . . . . . C (0) P /M . . . . S Figure 10:
Phase diagram for extinction/survival with the two-dimensional initial condition . In(a)–(c) we show three different initial conditions where C (0) = Bw /L , and we vary w . (d) Phase diagram ofa square mesh with 51 ×
40 nodes for C (0) ∈ [1 / , /
10] and
P/M ∈ [1 / , / M = 1. The blackcurve indicates the survival/extinction threshold from the continuum model and the blue shading shows thesurvival probability S measured by 40 identically-prepared realisations. Vertical red line is C (0) = 0 . A = 0 . C (0), P/M and the initial arrangement of the population in a complicated manner. Whilestochasticity plays a role in the discrete model when parameters are close to the boundarythat separates the eventual survival and extinction in the continuum model. To highlight theinfluence of dimensionality on the fate of the population, we compare the outcomes from thecontinuum model in Figure 11(a) where we superimpose the boundaries that separate regionsof survival and extinction for the zero-dimensional initial condition (red), one-dimensionalinitial condition (black) and the two-dimensional initial condition (green).Superimposing these curves divides the ( C (0) , P/M ) space into five regions, R , R ,..., R , with different long-term outcomes depending on the dimensionality of the initial condition.It is interesting to note that there is no influence of dimensionality above the Allee threshold,see region R in Figure 11(a). Populations with initial conditions where C (0) > A will alwayseventually survive as the growth mechanism leads to birth events only. This suggests thatthe influence of dimensionality plays a role only when C (0) < A . To emphasise the differentlong-term outcomes we compare solutions of the continuum model with different values of C (0) and P/M in Figures 11(b)–(g). The solutions in Figures 11(b)–(g) correspond to variouskey choices of C (0) and P/M . For example, the profiles in Figure 11(b) relate to region R where all outcomes lead to extinction regardless of the dimensionality of the initial condition.In contrast, profiles in Figure 11(c) relate to region R where we see extinction for the zero-and two-dimensional initial conditions, whereas the one-dimensional initial condition leadsto survival. This is a very interesting outcome since the discrete mechanism in these threecases is identical, yet the long-term outcome is very different, and depends only upon thedimensionality of the initial condition. Similar differences are highlighted in Figures 11(d)–(g)where we see that differences in long-term survival or extinction of the population often dependon the dimensionality of the initial distribution.25 . . . . . . . . . . . R R R R R C (0) P/M (a) R R R R R . C ( T ) , C ( T ) T (b) . C ( T ) , C ( T ) T (c) . C ( T ) , C ( T ) T (d) . C ( T ) , C ( T ) T (e) . C ( T ) , C ( T ) T (f) . C ( T ) , C ( T ) T (g) Figure 11:
Role of dimensionality in long-term survival and extinction. (a) The combined phasediagrams of Figures 8–10 where the red, black and green curves highlight the boundaries between extinctionand survival for the zero-, one- and two-dimensional initial condition, respectively. Profiles in (b)–(g) showprofiles of C ( T ) and C ( T ) for six choices of P/M and C (0). The profile colours in (b)–(g) correspond to thecoloured discs shown in (a). . Conclusions and Outlook In this work we design, analyse and implement a new two-dimensional stochastic discretemodel incorporating movement, birth and death events with crowding effects to study popula-tion extinction. The continuum limit of the discrete model is a nonlinear RDE which can beused to study a wide range of macroscopic phenomena including linear diffusion, nonlineardiffusion, as well as logistic and bistable growth kinetics. Since the aim of this work is tofocus on long term survival or extinction, we choose the movement crowding function to be G ( C ) = 1 − C which corresponds to macroscopic linear diffusion. We choose the growthcrowding function to be F ( C ) = a (1 − C )( C − A ) which leads to a classical cubic bistablesource term with Allee threshold A . Using a range of initial conditions, we show that numericalsolutions of the continuum RDE compare well with appropriately averaged data from thediscrete model.The focus of our work is to use the discrete and continuum models to explore the factors thatinfluence the long-term fate of the bistable population. In particular, we explore three differentinitial conditions on a finite L × L domain with periodic boundary conditions. The zero-dimensional initial condition involves distributing agents evenly across the entire L × L domain,the one-dimensional initial condition involves distributing agents along a vertical strip withinthe L × L domain so that the initial density is independent of vertical position in the domain,and the two-dimensional initial condition involves distributing agents within a central squareregion within the L × L domain. Our results show that the dimensionality of the populationinfluences the long-term fate of the population in a way that is often overlooked. For example,we show that in many circumstances two populations that are identical with the exceptionof the dimensionality of their initial condition can lead to very different survival/extinctionoutcomes. This suggests the importance of considering the influence of spatial arrangements ofindividuals in studies of population dynamics. Standard methods of relating long-term survivalor extinction as a function of the initial population number neglects the role of dimension,whereas we show that otherwise identical populations can either survive or go extinct, withthe only difference being the initial spatial arrangement of the population.There are many avenues for extending the work presented in this study. All results presentedhere focus on a particular choice of domain size L , and the two-dimensional initial conditionsare restricted to the case where the initial population is distributed within a square region.Both of these features can be relaxed and similar numerical explorations of the long-term27urvival or extinction of the populations can be conducted using the software provided onGitHub for both the continuum and discrete models. Another feature of this work that couldbe explored is the choice of crowding functions. As we pointed out, all simulations here focuson G ( C ) = 1 − C , which gives rise to linear diffusion, and F ( C ) = a (1 − C )( C − A ) which givesrise to the classical cubic bistable term. Other choices of G ( C ) and F ( C ) can be incorporatedinto the discrete model to explore how the results presented here depend upon the precisedetails of these choices of crowding functions. We note that other choices of G ( C ) lead todifferent motility mechanisms that are associated with nonlinear diffusion mechanisms, andthat these can be important for applications where adhesion (Deroulers et al. 2009) and inertialeffects (Zhang et al. 2019) are relevant. While we have not explicitly explored these effects inthis work, our framework is sufficiently general that these mechanisms can be incorporated andexplore, if required. Another interesting extension would be to consider Allee-type dynamicswithin more complicated populations that are composed of interacting subpopulations (Simp-son et al. 2009). Under these conditions interactions among various subpopulations can alsocontribute to the eventual survival or extinction of the population (Taylor et al. 2020; Krauseand Van Gorder 2020). Acknowledgements:
This work is supported by the Australian Research Council (DP200100177,DE200100988, DP190102545). 28 ppendix A. Algorithm for discrete simulationsAlgorithm 1:
Pseudo-code for a single realisation of the stochastic model Create a two-dimensional I × J hexagonal lattice; Distribute agents withspecific initial conditions; The total number of lattice site is IJ ; Set t = 0; Calculate total agents Q ( t ); while t < t end and Q ( t ) > Q ( t ) ≤ IJ do t = t + τ ; Q ( t ) = Q ( t − τ ); B = 0; B = 0; Draw two random variables: β ∼ U [0 , β ∼ U [0 , while B < Q ( t ) do B = B + 1; Randomly choose an agent s ; if β < M then Calculate ¯ K (m) s and G (¯ K (m) s ); Draw a random variable: γ ∼ U [0 , if γ < G (¯ K (m) s ) then Randomly choose a vacant site in N ( s ) and move agent to chosensite else Nothing happens; end else Nothing happens; end end while B < N ( t ) do B = B + 1; Randomly choose an agent s ; if β < P then Calculate K (g) s and F (¯ K (g) s ); Calculate a random variable: γ ∼ U [0 , if F (¯ K (g) s ) > then if γ < F (¯ K (g) s ) then Randomly choose a vacant site in N ( s ) and place a new agenton chosen site; Q ( t ) = Q ( t ) + 1 else if F (¯ K (g) s ) < then if γ < − F (¯ K (g) s ) then Remove agent; Q ( t ) = Q ( t ) − else Nothing happens; end else Nothing happens; end end end ppendix B. Derivation of the continuum limit We recall Equation (5), that is, the expected change in occupancy of site s during the timeinterval from t to t + τ , δ (¯ C s ) = M |N | (1 − ¯ C s ) X s ∈N { s } ¯ C s G (¯ K (m) s )1 − ¯ K (m) s − M ¯ C s G (¯ K (m) s )+ P |N | (1 − ¯ C s ) X s ∈N { s } H ( F (¯ K (g) s ))¯ C s F (¯ K (g) s )1 − ¯ K (g) s − (1 − H ( F (¯ K (g) s )) P ¯ C s F (¯ K (g) s ) . (B.1)As we know that the continuum limit of the last two terms in Equation (B.1) leads to a sourceterm λCF ( C ) (Jin et al. 2016), we focus on the movement mechanism, that is, the first twoterms on the right hand side of Equation (B.1). For convenience, we will omit the overlines onnotations in the following content.It is useful to first write the general form of the Taylor series relating the occupancy ofsites ( x + a, y + b ), C x + a,y + b = C x,y + ( a ∆) ∂C x,y ∂x + ( b ∆) ∂C x,y ∂y + ( a ∆) ∂C x,y ∂x + 2 ab ∆ ∂C x,y ∂x∂y + ( b ∆) ∂C x,y ∂y + O (∆ ) . (B.2)We represent the six nearest neighbouring sites of site s located at ( x, y ) as site s with ( x − ∆ , y );site s with ( x + ∆ , y ); site s with ( x − ∆ / , y + ∆ √ / s with ( x + ∆ / , y + ∆ √ / s with ( x − ∆ / , y − ∆ √ /
2) and site s with ( x + ∆ / , y − ∆ √ / N = { s , s , s , s , s , s } . The truncated Taylor series of these sites are C s = C s − ∂C s ∂x ∆ + ∂ C s ∂x ∆ O (∆ ) , (B.3) C s = C s + ∂C s ∂x ∆ + ∂ C s ∂x ∆ O (∆ ) , (B.4) C s = C s − ∂C s ∂x ∆2 + ∂C s ∂y √ " ∂ C s ∂x + 34 ∂ C s ∂y − √ ∂ C s ∂x∂y ∆ O (∆ ) , (B.5) C s = C s + ∂C s ∂x ∆2 + ∂C s ∂y √ " ∂ C s ∂x + 34 ∂ C s ∂y + √ ∂ C s ∂x∂y ∆ O (∆ ) , (B.6) C s = C s − ∂C s ∂x ∆2 − ∂C s ∂y √ " ∂ C s ∂x + 34 ∂ C s ∂y + √ ∂ C s ∂x∂y ∆ O (∆ ) , (B.7) C s = C s + ∂C s ∂x ∆2 − ∂C s ∂y √ " ∂ C s ∂x + 34 ∂ C s ∂y − √ ∂ C s ∂x∂y ∆ O (∆ ) . (B.8)30he local density of s is obtained by summing the Taylor series of sites in N { s } , that is, K (m) s = 16 X s ∈N { s } C s = C s + ∂ C s ∂x + ∂ C s ∂y ! ∆ O (∆ ) . (B.9)Similarly, the local density of s is obtained by summing the Taylor series of sites in N { s } ,that is, K (m) s = 16 X s ∈N { s } C s = C s + ∂ C s ∂x + ∂ C s ∂y ! ∆ O (∆ ) , = C s − ∂C s ∂x ∆ + ∂ C s ∂x ∆ ∂ C s ∂x + ∂ C s ∂y ! ∆ O (∆ ) . (B.10)For simplification we rewrite Equation (B.10) as K (m) s = C s + e C s , where e C s ∼ O (∆).Subsequently, the movement crowding function at s can be expanded as G (cid:16) K (m) s (cid:17) = G (cid:16) C s + e C s (cid:17) , = G ( C s ) + d G ( C s )d C e C s + d G ( C s )d C e C s . (B.11)The expansions of G ( K (m) s ), G ( K (m) s ),..., G ( K (m) s ) have similar forms to (B.11). We then goback to the first term on the right hand side of (B.1), which gives M − C s ) X s ∈N { s } C s G ( K (m) s )1 − K (m) s . (B.12)For convenience we further drop the s notation so that C s becomes C and C s becomes C .Subsequently, (B.12) becomes M − C ) X i =1 C i G ( K (m) s i )1 − K (m) s i . (B.13)Moreover, we will use two notations A = ∂ C s ∂x + ∂ C s ∂y ! ∆ , B = (cid:18) ∂C s ∂x (cid:19) + (cid:18) ∂C s ∂y (cid:19) ! ∆ , (B.14)31n the following content. Expanding the term related to site s in (B.13) gives M − C ) (cid:16) C + e C − A (cid:17) G ( C ) + G ( C ) e C + G ( C ) e C ! − (cid:16) C + e C (cid:17) = M − C ) (cid:16) C + e C − A (cid:17) G ( C ) + G ( C ) e C + G ( C ) e C ! − C + e C (1 − C ) + e C (1 − C ) ! + O (∆ )= M (cid:16) C + e C − A (cid:17) G ( C ) + G ( C ) e C + G ( C ) e C ! e C − C + e C (1 − C ) ! + O (∆ )= M (cid:20) CG ( C ) + (cid:18) CG ( C ) + G ( C )1 − C (cid:19) e C + (cid:18) G ( C )(1 − C ) + G ( C )1 − C + CG ( C )2 (cid:19) e C − G ( C ) A (cid:21) + O (∆ ) . The terms related to other sites can be obtained in a similar way. Therefore, expanding allterms in (B.13) and neglecting terms of order O (∆ ) gives M " CG ( C ) + (cid:18) CG ( C ) + G ( C )1 − C (cid:19) X k =1 e C k + (cid:18) G ( C )(1 − C ) + G ( C )1 − C + CG ( C )2 (cid:19) X k =1 e C k − G ( C ) A . (B.15)Furthermore, since we have X k =1 e C k = 12 ∂ C∂x + ∂ C∂y ! ∆ O (∆ ) , = 12 A + O (∆ ) , (B.16)and X k =1 e C k = 12 (cid:18) ∂C∂x (cid:19) + (cid:18) ∂C∂y (cid:19) ! ∆ O (∆ ) , = 12 B + O (∆ ) , (B.17)Equation (B.15) becomes M CG ( C )+ M (cid:18) CG ( C ) − G ( C ) + 2 G ( C )1 − C (cid:19) A + M (cid:18) CG ( C ) + 2 G ( C )(1 − C ) + 2 G ( C )1 − C (cid:19) B + O (∆ ) . (B.18)Remind that the second term in (B.1) is M CG (¯ K (m) s ) = M CG ( C ) + M CG ( C ) e C, = M CG ( C ) + M CG ( C ) A + O (∆ ) . (B.19)32hen combining (B.18) and (B.19) gives δ ( C s ) = (cid:18) CG ( C ) − G ( C ) + 2 G ( C )1 − C (cid:19) M A + (cid:18) CG ( C ) + 2 G ( C )(1 − C ) + 2 G ( C )1 − C (cid:19) M B + O (∆ ) , = (cid:18) CG ( C ) + 1 + C − C G ( C ) (cid:19) M A + (cid:18) CG ( C ) + 2 G ( C )(1 − C ) + 2 G ( C )1 − C (cid:19) M B + O (∆ ) . (B.20)Dividing both sides of the resulting expression by τ , and letting ∆ → τ → /τ held constant, leads to the following nonlinear reaction-diffusion equation, ∂C∂t = D ∇ · (cid:20)(cid:18) CG ( C ) + 1 + C − C G ( C ) (cid:19) ∇ C (cid:21) + λCF ( C ) , (B.21)where D = M ∆ ,τ → ∆ τ , λ = lim τ → Pτ . (B.22)If we define D ( C ) = CG ( C ) + 1 + C − C G ( C ) , (B.23)then the continuum limit is written as ∂C∂t = D ∇ · [ D ( C ) ∇ C ] + λCF ( C ) . (B.24)33 ppendix C. Numerical methods Here, we introduce the method of lines to numerically calculate solutions of the PDE ∂C∂t = D ∇ C + R ( C ) , (C.1)on a square domain Ω = { ( x, y ) , < x < L, < y < L } . We first discretise the spatialderivative in Equation (C.1) with an ( I + 1) × ( I + 1) mesh. Nodes on the mesh are uniformlydistributed with spacing δx > x i and y j with i = 0 , , , ..., I and j =0 , , , ..., I satisfying I = L/δx . We leave the time derivative continuous and obtaind C i,j d t = D δx ( C i +1 ,j + C i − ,j + C i,j +1 + C i,j − − C i,j ) + R ( C i,j ) . (C.2)This equation is valid for interior nodes, and is modified on the boundary nodes to simulateperiodic boundary conditions. This system of I × I coupled ordinary differential equations isthen integrated through time using MATLABs function ode45 (MATLAB 2020). Followingsimilar steps, we can also calculate the numerical solution of the PDE ∂C∂t = D ∂ C∂x + R ( C ) . (C.3)34 ppendix D. Phase diagrams with B = 1 Instead of varying the size of H , we now vary C (0) by varying B . We constrain the region H as a vertical strip with width w = 64, as shown in Figure S1(a). As C (0) = Bw /L , wherewe fix w /L = 0 .
64, the initial density C (0) varies from 0 .
192 to 0 .
64 when B varies from0 . P/M = 4 λ/D by holding M = 1 andvarying P ∈ [1 / , / C (0) , P ) space into a square mesh with36 ×
40 nodes. Figure S1(d) shows a phase diagram illustrating how the survival probability, S , depends upon C (0) and P/M . The boundary that separates the eventual survival andextinction in the continuum model is shown in solid black, and the survival probability fromthe discrete simulations is shown in blue shading. The long–term predictions in terms ofsurvival or extinction are consistent between the continuum and discrete models. Furthermore,we observe a different phenomenon compared to the results in Figure 9: There is a lowerbound on C (0) for survival in Figure S1(d). This lower bound relates to B = 0 .
4, and indicatesthe Allee threshold A = 0 .
4. Unlike the solid line, C (0) = 0 .
4, which indicates a thresholdof survival in the sense of global density, the dashed line, B = 0 .
4, indicates a threshold ofsurvival in the sense of local density. 35 = 0 . L L (a) B = 0 . L L (b) B = 1 L L (c)(d) . . . . . C (0) P /M . . . . S B Figure S1:
Phase diagram for extinction/survival with the one-dimensional initial condition . In(a)–(c) we show three different initial conditions where C (0) = Bw /L , and we fix w = 64 and vary B . (d)Phase diagram on a square mesh with 36 ×
40 nodes for C (0) ∈ [0 . , .
64] and
P/M ∈ [1 / , / S measured by 40 identically-prepared realisations. Vertical solid red line is C (0) = 0 . B = 0 .
4. They both relate to the Allee threshold, A = 0 . H as a square region with width w = w = w = 80 inFigures S2(a)–(c). As C (0) = Bw /L , where we fix w /L = 0 .
64, the initial density C (0)again varies from 0 .
192 to 0 .
64 when B varies from 0 . P/M = 4 λ/D by holding M = 1 and varying P ∈ [1 / , / C (0) , P/M ) spaceinto a square mesh with 36 ×
40 nodes. With this initial condition we again construct aphase diagram summarising the long–term survival outcome as a function of C (0) and P/M in Figure S2(d), which is very similar to the phase diagram in Figure S1(d) where we see thatthe long-term survival depends on
P/M and the two red lines indicating the Allee threshold.37 = 0 . L L (a) B = 0 . L L (b) B = 1 L L (c)(d) . . . . . C (0) P /M . . . . S B Figure S2:
Phase diagram for extinction/survival with the one-dimensional initial condition . In(a)–(c) we show three different initial conditions where C (0) = Bw /L , and we fix w = 80 and vary B . (d)Phase diagram on a square mesh with 36 ×
40 nodes for C (0) ∈ [0 . , .
64] and
P/M ∈ [1 / , / S measured by 40 identically-prepared realisations. Vertical solid red line is C (0) = 0 . B = 0 .
4. They both relate to the Allee threshold, A = 0 .
38o highlight the different fates of various populations, we compare the outcomes fromthe continuum model from three phase diagrams in Figure S3(a), where we superimpose theboundaries that separate regions of survival and extinction for the zero-dimensional initialcondition (red), one-dimensional initial condition (black) and the two-dimensional initialcondition (green). Superimposing these curves divides the ( C (0) , P/M ) plane into four regionswith different long-term outcomes depending on the dimensionality of the initial conditions.To emphasise these differences we compare solutions of the continuum model with differentvalues of C (0) and P/M in S3(b)–(g). The solutions in Figures S3(b)–(g) correspond tovarious illustrative choices of C (0) and P/M . For example, the profiles in Figure S3(b) relatedto region R all lead to extinction regardless of the dimensionality of the initial condition,whereas the profiles in Figure S3(c) related to region R lead to extinction for the zero- andtwo-dimensional conditions, whereas the one-dimensional initial condition leads to survival.39 . . . . . . . . . . R R R R C (0) P/M (a) R R R R . C ( T ) , C ( T ) T (b) . C ( T ) , C ( T ) T (c) . C ( T ) , C ( T ) T (d) . C ( T ) , C ( T ) T (e) . C ( T ) , C ( T ) T (f) . C ( T ) , C ( T ) T (g) Figure S3:
Role of dimensionality in long-term survival and extinction. (a) Combined phase diagramwhere the red, black and green curves highlight the boundaries between extinction and survival for the zero-,one- and two-dimensional initial condition, respectively. Profiles in (b)–(g) show profiles of C ( T ) and C ( T )for six different choices of P/M and C (0). The profile colours in (b)–(g) correspond to the coloured discssuperimposed in (a). eferences Allee, W.C., Bowen, E.S., 1932. Studies in animal aggregations: Mass protection againstcolloidal silver among goldfishes. Journal of Experimental Zoology. 61, 185–207.Arroyo-Esquivel, J., Hastings, A., 2020. Spatial dynamics and spread of ecosystem engineers:Two patch analysis. Bulletin of Mathematical Biology. 82, 149.Baker, R.E., Simpson, M.J., 2010. Correcting mean-field approximations for birth-death-movement processes. Physical Review E. 82, 041905.B¨ottger, K., Hatzikirou, H., Voss-B¨ohme, A., Cavalcanti-Adam, E.A., Herrero, M.A., Deutsch,A., 2015. An emerging Allee effect is critical for tumor initiation and persistence. PLoSComputational Biology. 11, e1004366.Chaplain, M.A.J., Lorenzi, T., Macfarlane, F.R., 2020. Bridging the gap between individual-based and continuum models of growing cell populations. Journal of Mathematical Biology.80, 343–371.Courchamp, F., Berec, L., Gascoigne, J., 2008. Allee effects in ecology and conservation.Oxford University Press, Oxford.Courchamp, F., Clutton-Brock, T., Grenfell, B., 1999. Inverse density dependence and theAllee effect. Trends in Ecology & Evolution. 14, 405–410.Deroulers, C., Aubert, M., Badoual, M., Grammaticos, B., 2009. Modeling tumor cell migration:From microscopic to macroscopic models. Physical Review E. 79, 031917.Drake, J.M., 2004. Allee effects and the risk of biological invasion. Risk Analysis. 24, 795–802.Edelstein-Keshet, L., 2005. Mathematical models in biology. SIAM, Philadelphia.Fadai, N.T., Johnston, S.T., Simpson, M.J., 2020. Unpacking the Allee effect: determiningindividual-level mechanisms that drive global population dynamics. Proceedings of theRoyal Society A: Mathematical, Physical and Engineering Sciences. 476, 20200350.Fadai, N.T., Simpson, M.J., 2020. Population dynamics with threshold effects give rise to adiverse family of Allee effects. Bulletin of Mathematical Biology. 82, 74.41ife, P.C., 1979. Long time behavior of solutions of bistable nonlinear diffusion equationsn.Archive for Rational Mechanics and Analysis. 70, 31–36.Fisher, R.A., 1937. The wave of advance of advantageous genes. Annals of Eugenics. 7,355–369.Hastings, A., Cuddington, K., Davies, K.F., Dugaw, C.J., Elmendorf, S., Freestone, A.,Harrison, S., Holland, M., Lambrinos, J., Malvadkar, U., Melbourne, B.A., Moore, K.,Taylor, C., Thomson, D., 2005. The spatial spread of invasions: new developments in theoryand evidence. Ecology Letters. 8, 91–101.Hughes, B.D., 1995. Random walks and random environments: random walks. volume 1.Oxford University Press, Oxford.Jin, W., Penington, C.J., McCue, S.W., Simpson, M.J., 2016. Stochastic simulation tools andcontinuum models for describing two-dimensional collective cell spreading with universalgrowth functions. Physical Biology. 13, 056003.Johnston, S.T., Baker, R.E., McElwain, D.L.S., Simpson, M.J., 2017. Co-operation, competitionand crowding: A discrete framework linking Allee kinetics, nonlinear diffusion, shocks andsharp-fronted travelling waves. Scientific Reports. 7, 42134.Johnston, S.T., Simpson, M.J., Crampin, E.J., 2020. Predicting population extinction in lattice-based birth-death-movement models. Proceedings of the Royal Society A: Mathematical,Physical and Engineering Sciences. 476, 20200089.Kot, M., 2001. Elements of mathematical ecology. Cambridge University Press, Cambridge.Krause, A.L., Van Gorder, R.A., 2020. A non-local cross-diffusion model of populationdynamics II: Exact, approximate, and numerical traveling waves in single- and multi-speciespopulations. Bulletin of Mathematical Biology. 82, 113.Lewis, M.A., Kareiva, P., 1993. Allee dynamics and the spread of invading organisms.Theoretical Population Biology. 43, 141–158.Macfarlane, F.R., Lorenzi, T., Chaplain, M.A.J., 2018. Modelling the immune response tocancer: An individual-based approach accounting for the difference in movement betweeninactive and activated T cells. Bulletin of Mathematical Biology. 80, 1539–1562.42aini, P.K., McElwain, D.L.S., Leavesley, D., 2004a. Traveling wave model to interpret awound-healing cell migration assay for human peritoneal mesothelial cells. Tissue Engineering.10, 475–482.Maini, P.K., McElwain, D.L.S., Leavesley, D., 2004b. Travelling waves in a wound healingassay. Applied Mathematics Letters. 17, 575–580.MATLAB, 2020. ode45 documentation. URL:
Murray, J.D., 2002. Mathematical biology: I. An introduction. Springer, New York.Neufeld, Z., von Witt, W., Lakatos, D., Wang, J., Hegedus, B., Czirok, A., 2017. The role ofAllee effect in modelling post resection recurrence of glioblastoma. PLoS ComputationalBiology. 13, e1005818.Saltz, D., Rubenstein, D.I., 1995. Population dynamics of a reintroduced asiatic wild ass(