Population Dynamics with Threshold Effects Give Rise to a Diverse Family of Allee Effects.

The Allee effect describes populations that deviate from logistic growth models and arises in applications including ecology and cell biology. A common justification for incorporating Allee effects into population models is that the population in question has altered growth mechanisms at some critical density, often referred to as a threshold effect. Despite the ubiquitous nature of threshold effects arising in various biological applications, the explicit link between local threshold effects and global Allee effects has not been considered. In this work, we examine a continuum population model that incorporates threshold effects in the local growth mechanisms. We show that this model gives rise to a diverse family of Allee effects, and we provide a comprehensive analysis of which choices of local growth mechanisms give rise to specific Allee effects. Calibrating this model to a recent set of experimental data describing the growth of a population of cancer cells provides an interpretation of the threshold population density and growth mechanisms associated with the population.

Introduction

Mathematical models of population dynamics often include an Allee effect to account for dynamics that deviate from logistic growth (Stephens et al. 1999; Allee and Bowen 1932; Courchamp et al. 1999; Taylor and Hastings 2005; Courchamp et al. 2008). The logistic growth model (Table 1, Fig. 1) describes the growth rate, , as a quadratic function of density, C(t), at time . The logistic growth model has two equilibria: and , where an equilibrium is any value such that when . Since densities near will approach K, while densities near diverge away from zero (Fig. 1), we say that is a stable equilibrium, while is an unstable equilibrium. This means that the logistic growth model implicitly assumes that all densities, no matter how small, eventually thrive.

Table 1

Typical mathematical descriptions of logistic growth, the Weak Allee effect, and the Strong Allee effect

Effect	Typical mathematical description	Notes
Logistic growth	dC(t)dt=rC(t)1-C(t)K	r>0,K>0
Weak Allee	dC(t)dt=rC(t)1-C(t)K1+C(t)A	r>0,A>0,K>0
Strong Allee	dC(t)dt=rC(t)1-C(t)KC(t)B-1	r>0,0<B<K

Fig. 1

Comparison of typical logistic growth, Weak Allee, and Strong Allee models. The mathematical descriptions of the three models are shown in Table 1

Mathematical models that include an Allee effect relax the assumption that all population densities will thrive and survive, which is inherent in logistic growth models (Murray 2003; Edelstein-Keshet 2005; Stephens et al. 1999; Taylor and Hastings 2005; Courchamp et al. 2008). Consequently, populations described using Allee effect models exhibit more complicated and nuanced dynamics, including reduced growth at low densities (Neufeld et al. 2017; Johnson et al. 2006; Gerlee 2013) and extinction below a critical density threshold (Courchamp et al. 1999; Allee and Bowen 1932; Taylor and Hastings 2005; Courchamp et al. 2008). The phrase Allee effect can have many different interpretations in different parts of the literature. For instance, the Weak Allee effect (Table 1, Fig. 1) is used to describe density growth rates that deviate from logistic growth, but do not include additional equilibria (Murray 2003; Edelstein-Keshet 2005; Taylor and Hastings 2005; Stephens et al. 1999). A common mathematical description of the Weak Allee effect is shown in Table 1, where the factor represents the deviation from the classical logistic growth model. Despite the similarity between logistic growth and the Weak Allee effect, it is not possible to write down an explicit solution for Weak Allee effect in terms of C(t), like we can for logistic growth. Despite this, we are still able to examine the equilibria of the Weak Allee effect to understand its salient features. Since , the Weak Allee effect does not incorporate any additional equilibria other than and . Noting that the main feature of an Allee effect is a deviation from logistic growth, the cubic representation of the growth rate shown in Table 1 is employed predominantly for simplicity rather than explicit biological significance (Taylor and Hastings 2005; Stefan et al. 2012; Stephens et al. 1999). Therefore, in this work, we refer to the Weak Allee effect as any population density growth rate that deviates from logistic growth without incorporating additional equilibria.

Another common type of Allee effect is the Strong Allee effect (Table 1, Fig. 1), in which an additional unstable intermediate equilibrium, with , is incorporated (Murray 2003; Edelstein-Keshet 2005; Taylor and Hastings 2005; Stephens et al. 1999; Courchamp et al. 1999). In a similar fashion to the Weak Allee effect, the cubic form of the Strong Allee effect (Table 1) is chosen predominantly for simplicity (Taylor and Hastings 2005; Stefan et al. 2012; Stephens et al. 1999). Therefore, we will refer to any growth rate with two stable equilibria, and , and an additional intermediate unstable equilibrium as the Strong Allee effect. Throughout this work, we refer to growth rates that deviate from logistic growth as an Allee effect, whereas specific Allee effects (e.g. the Weak Allee effect and the Strong Allee effect) are referred to using more specific terminologies.

While Allee effects were originally used to describe population dynamics arising in ecology (Taylor and Hastings 2005; Tu et al. 2019; Courchamp et al. 1999; Johnson et al. 2006; Simberloff et al. 2013; Seebens et al. 2017; Drake 2004; Courchamp et al. 2008), there has been increasing interest in examining the potential for Allee effects in population dynamics relating to cell biology (Neufeld et al. 2017; Böttger et al. 2015; Gerlee 2013; Sarapata and de Pillis 2014; Jenner et al. 2018, 2019; Bobadilla et al. 2019; Johnston et al. 2017; Jin et al. 2017; Johnson et al. 2019; de Pillis et al. 2005; de Pillis and Radunskaya 2003). In both cell biology and ecological applications, the Allee effect provides a suitable modelling framework to describe the dynamics of well-mixed populations that exhibit non-logistic features. However, because standard models incorporating Allee effects are continuum models that describe global, population-level features of the population dynamics, the interpretation of Allee effects at the individual scale remains less clear (Johnston et al. 2017; Böttger et al. 2015).

Understanding how local, stochastic growth mechanisms give rise to global Allee effects in a population is important, since these individual-level mechanisms can ultimately determine whether a population will survive or be driven to extinction (Johnston et al. 2017; Scott et al. 2014; Colon et al. 2015; Böttger et al. 2015). Certain individual-level biological features are ubiquitous among populations displaying Allee effects, providing a unifying feature in both cell biology and ecological applications. One of these phenomena is a threshold effect (Frankham 1995; Rossignol et al. 1999; Metzger and Décamps 1997), which we also refer to as a binary switch. We define a binary switch as a local feature of a population that behaves differently when a particular biological mechanism is present or absent. Some examples of binary switches include: the go-or-grow hypothesis in cell biology (Hatzikirou et al. 2012; Vittadello et al. 2020), phenotypic plasticity (Friedl and Alexander 2011; Böttger et al. 2015), tree masting (Koenig and Knops 2005), external harvesting pressure (Courchamp et al. 1999; Kuparinen et al. 2014), density-dependent clustering (Martínez-García et al. 2015), and resource depletion (Hopf and Hopf 1985). For all of these examples, Allee effects have been proposed to potentially explain more complicated and nuanced population dynamics than are possible in a logistic growth framework. However, the link between the details of such a local binary switch and the resulting population-level Allee effect is unclear. Given that local binary switches are thought to be widely important in biology and ecology, we ask two questions: (i) how does the incorporation of a local binary switch in proliferation and death rates affect the global dynamics of a population? and (ii) how does this local binary switch relate to different forms of Allee effects?

In this work, we show that incorporating local-level binary switches in a continuum, population-level mathematical modelling framework gives rise to a surprisingly diverse family of Allee effects. Some switches in proliferation and death rates give rise to established Allee effects, whereas other binary switches lead to more generalised Allee effects that have not been previously reported. We show that incorporating local-level binary switches in proliferation and death rates leads to a diverse family of Allee effects with only a few model parameters. This model, which we refer to as the Binary Switch Model, captures key biological features, but continues to exhibit the same qualitative features as various Allee effects. We conclude by applying the Binary Switch Model to a recent cell biology data set. Interpreting this data with our modelling framework suggests that the observed growth is non-logistic and that the phenomena are best explained by a binary switch at low density.

The Binary Switch Model

We consider an individual-based model framework that incorporates individual-level growth mechanisms varying with local population density to describe the temporal evolution of the global population density. One framework incorporating these aforementioned features is the stochastic agent-based model framework that we proposed in Fadai et al. (2020), in which individuals of the same size move, die, and proliferate on a two-dimensional hexagonal lattice. This discrete model incorporates exclusion (crowding) effects, allowing the population density to saturate at a finite capacity, as well as proliferation and death rates that vary with the local population density. While local population density can be measured in many different ways, Fadai et al. (2020) take the simplest approach and use the number of nearest neighbours as a measure of local density (Fig. 2).

Fig. 2

Schematic for the Binary Switch Model. Individuals in a population a can sense nearby individuals, providing a simple measure of local density. Individuals who sense higher than a threshold density, M, are shown in blue, while more isolated individuals are shown in red. This threshold density determines the constant rates at which individuals proliferate and die. b, c The binary switch shown here occurs when individuals can sense more than neighbours

As the individual dynamics of the stochastic agent-based model are difficult to analyse mathematically, we examine the continuum limit per-capita growth rate as a means of representing the average dynamics of the spatially uniform population, noting that there is good agreement between these two modelling approaches (Fig. 3). Full details of the discrete–continuum comparison are summarised in the Supplementary Information. Since the average population dynamics obtained from the discrete stochastic individual-based model agree well with its continuum description (Fig. 3), we will only consider the features of the continuum description of the model, whose per-capita growth rate is given bywhereHere, C(t) is the population density at time t, while and are the proliferation and death rates that vary with the number of nearest neighbours, n (Fadai et al. 2020). The parameter grouping can be interpreted as the net growth mechanism for a particular local population density. Noting that represents the maximum packing density, we have . Equation (1) has a thirteen-dimensional parameter space: namely, .

Fig. 3

a When no binary switch is present, (1) reduces to logistic growth. b, c When a binary switch is incorporated in proliferation and death rates (), the continuum limit is no longer logistic. In all of these parameter regimes, the average density data determined from discrete model simulations, shown in red dashed curves in the middle column (Supplementary Information), agrees well with the continuum limit predictions (4), shown in black solid curves. Density growth rates in the right-most column show that (a) is logistic, while (b, c) are not

We incorporate a binary switch into (1) by choosingThis choice of parameters means that we have the proliferation rate when the local density is at or below the critical density M, or when the local density is above M. We refer to as the threshold density. For simplicity, we assume that the death rates are a particular fraction of the proliferation rates: and . It is useful to note that (1)–(3) relaxes to the classical logistic growth model, for any choice of by setting and (Fig. 3a).

By substituting (3) into (1), we obtain the Binary Switch Model,whereis an indicator function. The Binary Switch Model shows, for the first time, how a local binary switch in individual-level proliferation and death rates leads to a particular global density growth rate. A summary of parameters and their particular biological interpretation is shown in Table 2. In particular, we note that the Binary Switch Model reduces the thirteen-dimensional parameter space in (1) to a five-dimensional parameter space: . This reduced parameter space means that the Binary Switch Model can be used with less risk of over-fitting than (1) (Warne et al. 2019). We will discuss further merits of this reduced parameter space when calibrating the Binary Switch Model to experimental data in Sect. 3.

Table 2

Summary of parameters used in the Binary Switch Model

Parameter	Biological interpretation
r∈[0,∞)	Low-density proliferation rate
R∈[0,∞)	High-density proliferation rate
α∈[0,1]	Ratio of low-density death rate to low-density proliferation rate
β∈[0,1]	Ratio of high-density death rate to high-density proliferation rate
M∈{0,1,2,3,4,5}	Threshold density

In Fig. 3, we show how the Binary Switch Model gives rise to non-logistic growth mechanisms. When no binary switch is present (Fig. 3a), the growth mechanisms are independent of local density and assume a single proliferation and death rate, resulting in logistic growth. However, when a binary switch is incorporated into the proliferation and death rates (Fig. 3b, c), the population dynamics described by (4) deviates from the classical logistic growth model. Consequently, we now wish to examine the various kinds of Allee effects the Binary Switch Model can give rise to. The main qualitative differences between logistic growth and various Allee effects are based on the number of equilibria and their stability; therefore, we now examine the roots of (4) for various parameter values. In all parameter regimes considered in the work, the zero equilibrium, , will always be present. Additional equilibria, if present, will be denoted as , where and are ordered such that for all i. Since the right-hand side of (4) is a sixth-degree polynomial, a maximum of six equilibria can be present in (0, 1], but explicit expressions for the solutions of the polynomial cannot be determined in general. We will show that in the Binary Switch Model, a maximum of three equilibria can be present in (0, 1]. Setting and (Case 1) or and (Case 2), we will show that fewer equilibria are present in (0, 1]. In Case 3, corresponding to and , certain combinations of parameter values produce equilibria with additional qualitative features, such as double-root and triple-root equilibria. For these special equilibria, we will designate particular symbols to , which appear as required.

Case 1: and

This case corresponds to situations where individuals below the threshold density M do not proliferate or die. We will now show that in Case 1, either no equilibria are present in (0, 1], or we have one equilibrium , depending on the choice of and M. In this regime, (4) simplifies toSince appears as a linear coefficient in (6), it is easier to solve for than for . The resulting relationship between and depends on the integer value of ; however, a general solution in terms of arbitrary M is difficult to obtain. Instead, we define the family of functions, , for a particular value of M, such thatUsing , we determine the unique value of that solves for a given value of , shown in Table 3. Plotting for all and indicates that is one-to-one on . Therefore, the inverse function also has one solution, provided that . This range of is obtained by mapping the interval (0, 1] via the functions . The functions in Table 3 provide a link between and : if is known, provides the parameter value to input in the model to obtain such an equilibrium. Conversely, if is known, Table 3 indicates whether or not . Finally, we note that when , or when , only the zero equilibrium, , is present.

Table 3

Relationships between the nonzero equilibrium of the Binary Switch Model, , to and M for Case 1 when (6)

M	β=f_M(C_1)	Range of β:C_1∈(0,1]
0	C15-6C14+15C13-20C12+15C_1-5C15-6C14+15C13-20C12+15C_1-6	β∈[0,5/6)
1	4C14-19C13+35C12-30C_1+105C14-24C13+45C12-40C_1+15	β∈[0,2/3)
2	-6C13+21C12-25C_1+10-10C13+36C12-45C_1+20	β∈[0,1/2)
3	4C12-9C_1+510C12-24C_1+15	β∈[0,1/3)
4	-C_1+1-5C_1+6	β∈[0,1/6)
5	0	∅

To determine the stability of the equilibria, we consider the cases when and when separately. When , two distinct equilibria are present: and . Based on the sign of at these equilibria, is always unstable and is always stable. These features are consistent with the Weak Allee effect, whereby the density growth rate deviates from logistic growth without incorporating additional equilibria. When , or when , is the only equilibrium and it is always stable, corresponding to the qualitative features of an extinction density growth rate, where for all C(0). Both qualitative features in this parameter regime are shown in the bifurcation diagram in Fig. 4. We conclude that in Case 1, either zero or one equilibria is present in the interval (0, 1], corresponding to extinction and Weak Allee parameter regimes, respectively.

Fig. 4

Bifurcation diagram of the Binary Switch Model, shown in (6), for Case 1 when . Varying produces different qualitative features in terms of equilibria and their stability. The resulting density growth rates, , are shown as a function of C, where a stable equilibrium is represented with a black circle and an unstable equilibrium with a white circle

Case 2: and

This case corresponds to when individuals aboveM do not proliferate or die. When , we havewhich is independent of . In a similar fashion to Case 1, we consider the equilibria for various choices of and M, noting that and are always equilibria in this case. However, we will show that in Case 2, we have the possibility of a third equilibrium in (0, 1). When this additional equilibria is present, then and ; otherwise, . To determine if is the first or second nonzero equilibrium, we defineand determine the value of that solves for a given value of , shown in Table 4. Like Case 1, the family of functions provide an explicit relationship between and . Since is one-to-one on , the inverse function also has one solution, , provided . This value of provides a third equilibrium of (8); conversely, when , or when , the only two equilibria are and .

Table 4

Relation between nonzero equilibrium, , to and M for Case 2 when (8)

M	α=g_M(C_1)	Range of α:C_1∈(0,1)
0	1	∅
1	4C_1+15C_1+1	α∈(5/6,1)
2	6C12+3C_1+110C12+4C_1+1	α∈(2/3,1)
3	4C13+3C12+2C_1+110C13+6C12+3C_1+1	α∈(1/2,1)
4	C14+C13+C12+C_1+15C14+4C13+3C12+2C_1+1	α∈(1/3,1)
5	1C15+C14+C13+C12+C_1+1	α∈(1/6,1)

In the case where , examining the sign of shows that and are unstable, whereas is stable. This combination of equilibria has the opposite stability properties of the Strong Allee effect (Table 1), and so we refer to density growth rates with these stability properties as the Reverse Allee effect. In the case where , or when , stability analysis shows that is stable and is unstable, which is consistent with the qualitative features of the Weak Allee effect. Finally, when , we return to having only two equilibria, and , but the stability is the opposite of the usual Weak Allee effect. Therefore, when , for . All these qualitative features in this parameter regime are shown in the bifurcation diagram in Fig. 5. We conclude that in Case 2, either one or two equilibria are present in (0, 1], with the Extinction regime occurring when . For , a new kind of Allee effect, which we call the Reverse Allee effect, occurs if two equilibria are present in (0, 1]; otherwise, we retrieve the Weak Allee effect.

Fig. 5

Bifurcation diagram of the Binary Switch Model, shown in (8), for Case 2 when . Varying produces different qualitative features in terms of equilibria and their stability. The resulting density growth rates, , are shown as a function of C, where a stable equilibrium is represented with a black circle and an unstable equilibrium with a white circle

Case 3: and

In the most general case, the proliferation and death rates of individuals change at the threshold density M, but remain nonzero on either side of the threshold density. As a result, (4) can be written asWithout loss of generality, we assume that , since other nonzero values or r can be rescaled to unity by changing the timescale in (4), which does not affect its equilibria. Consequently, with some rearranging, we haveWe will show that in Case 3, there can be between zero and three equilibria in (0, 1], noting that is an equilibrium of (11) if and only if . When we have three distinct equilibria in (0, 1], we obtain a new type of Allee effect, referred to here as the Hyper-Allee effect (Fadai et al. 2020), in which the zero equilibrium is unstable, and an intermediate unstable equilibrium is contained between two positive, stable equilibria. However, in order for the parameter space to continuously transition from the Weak Allee effect, as in Cases 1 and 2, to the Hyper-Allee effect, there must exist a critical set of model parameters at which a double-root equilibrium occurs. Therefore, to determine what regions of parameter space exhibit Hyper-Allee effects instead of the Weak Allee effect, we focus on determining the boundary of these effects in terms of model parameters and equilibria. This boundary, defined as the Tangential Manifold, will be the focus of our analysis in this section.

In addition to determining the boundary between Weak Allee and Hyper-Allee parameter spaces, we will also show that even more Allee effects are present when . In particular, we show that in Case 3, the Extinction parameter regime continues to exist, along with the Strong Allee effect, when . We also determine an explicit relationship between and M for when the Extinction regime becomes the Strong Allee effect, which is linked to the Tangential Manifold. We now focus our attention on determining additional equilibria .

Numerical observations indicate that certain combinations of can produce up to three distinct values of satisfying . Furthermore, in parameter regimes where three distinct equilibria are present in (0, 1], stability analysis about these equilibria reveals that and are unstable equilibria, whereas and are stable equilibria. These qualitative features are consistent with the aforementioned Hyper-Allee effect, which is a higher-order effect that is very different to the usual Weak Allee and Strong Allee effects (Fig. 6).

Fig. 6

Bifurcation diagram of the Binary Switch Model for Case 3, shown in (11), with and . Pairs of parameters produce different qualitative features, in terms of equilibria and their stability. The resulting density growth rates, , are shown as a function of C, where a stable equilibrium is represented with a black circle, an unstable equilibrium with a white circle, and a semi-stable equilibrium with a half-filled circle

For solutions to continuously transition from one equilibrium in (0, 1], like the Weak Allee effect in Cases 1 and 2, to three equilibria in (0, 1], such as the Hyper-Allee effect, we must have certain values of that produce a double root for . We denote this special case of a double root equilibrium as , which can occur in either the or equilibrium position. In addition to satisfying , the double root equilibrium, , must also satisfyThe set of parameters satisfying (12) is referred to as the Tangential Manifold, where the double root equilibrium, , is a semi-stable equilibrium of (11) (Strogatz 2018). A semi-stable equilibrium has the properties that populations slightly larger than remain close to , but populations slightly smaller than diverge away from , or vice-versa. Since we have two equations with four unknowns, we parametrise the Tangential Manifold as , for particular values of and (Fig. 6). The functions and describing the Tangential Manifold are shown in Table 5.

Table 5

Relation between the semi-stable equilibrium, , to , and M for Case 3. Parameter values satisfying and are members of the Tangential Manifold. If , then is a member of the Positive Tangential Manifold; if , then is a member of the Negative Tangential Manifold. The Triple Point, , is defined implicitly via , while the Junction Point, , is defined implicitly via

M	R=F_M(C^,β)	β=H_M(C)
0	0	∅
1	(C^-1)^6C^(C^^5-6C^^4+15C^^3-20C^^2-10C^-30β+20)	2(1-C)3
2	(C^-1)^5(6C^^2+8C^+1)C^^2(6C^^5-22C^^4+21C^^3+15C^^2+10C^+60β-30)	(1-C)(1+2C)3C+2
3	(C^-1)^4(6C^^4+16C^^3+21C^^2+6C^+1)C^^3(6C^^5-8C^^4-7C^^3-6C^^2-5C^-60β+20)	(1-C)(1+2C+2C^2)3C^2+4C+3
4	(C^-1)^3(C^^6+4C^^5+10C^^4+20C^^3+10C^^2+4C^+1)C^^4(C^^5+C^^4+C^^3+C^^2+C^+30β-5)	(1-C^2)(2C^2+C+2)3(C^3+2C^2+3C+4)
5	0	∅

While the Tangential Manifold can be determined explicitly by solving (12), we observe that two forms of a semi-stable equilibrium can occur (Fig. 6). If the double root is below some critical value, , then this semi-stable equilibrium occurs between , which is unstable, and some larger equilibrium , which is stable. If , then this semi-stable equilibrium is larger than both and , which remain unstable and stable, respectively. We refer to the branch of the Tangential Manifold where as the Positive Tangential Manifold, based on the sign of the density growth rate between and (Fig. 6). In a similar fashion, we refer to the branch of the Tangential Manifold where as the Negative Tangential Manifold. When , the double root becomes a stable triple root and satisfieswhere and are chosen to ensure we remain on the Tangential Manifold. Equation (13) provides an additional constraint on the Tangential Manifold, implying that we can relate to a unique value of . We denote if (13) is satisfied, with denoting the Triple Point of (11) (Table 5).

Additionally, from Fig. 6, we note that when , the equilibria and their resulting stability change, compared to . When , the Negative Tangential Manifold is valid for a unique pair of parameters, for a particular equilibrium value, . We define this critical equilibrium value as the Junction Point, which satisfiesWe denote if (14) is satisfied (Table 5); furthermore, we determine the corresponding value of R at the Junction Point by evaluating . When and , the only equilibrium value of (11) is , which is stable. This implies that all population densities go extinct in this parameter regime. When and , (11) has three solutions: , which is stable, an intermediate-valued unstable equilibrium , and a larger-valued stable equilibrium (Fig. 6). Thus, the stability features of this density growth rate are the same as the Strong Allee effect. When , the Junction Point, , is semi-stable, while remains stable. A summary of this diverse family of Allee effects, in terms of the number and stability of the equilibria, is shown in Table 6.

Table 6

Summary of qualitative features seen in the Binary Switch Model. The stability of each equilibrium, listed in increasing order of magnitude, can be stable (S), unstable (U), or semi-stable (SS)

Effect name	Equilibria	Stability	Notes
Extinction	{0}	{S}
Logistic growth	{0,C_1}	{U,S}	r=R,α=β
Weak Allee/Triple Point	{0,C_1}	{U,S}	Triple: C_1=C
Junction Point	{0,C_1}	{S,SS}	C_1=C~
Strong Allee	{0,C_1,C_2}	{S,U,S}
Reverse Allee	{0,C_1,C_2}	{U,S,U}	C_2=1
Positive Tangential Manifold	{0,C_1,C_2}	{U,SS,S}	C_1=C^
Negative Tangential Manifold	{0,C_1,C_2}	{U,S,SS}	C_2=C^
Hyper-Allee	{0,C_1,C_2,C_3}	{U,S,U,S}

From Table 5, we note some key features of the Tangential Manifold. Firstly, when , we note that the Triple Point is for . Since the Negative Tangential Manifold must have we conclude that the Negative Tangential Manifold does not exist when , which is also observed in Fig. 7. When and , the Triple Point and the Junction Point are both , implying that no points are contained in the Tangential Manifold. Consequently, parameter pairs that result in qualitative features other than the Extinction regime or the Weak Allee effect can only occur when and , as shown in Fig. 7. Finally, we note that when or , the Tangential Manifold does not exist, since the solution of (12) requires . Therefore, the qualitative features of (11) in the entire parameter space are those seen in the Weak Allee effect when and the Extinction regime when .

Fig. 7

Bifurcation diagram of the Binary Switch Model for Case 3, shown in (11), with and varying . The qualitative forms of various effects are shown in the legend, described in further detail in Fig. 6. The parameter space exhibiting Hyper-Allee features vanishes as

To summarise, we determine that in Case 3 when , and , a diverse family of Allee effects can be found. Among these Allee effects are: the Weak Allee effect, the Extinction regime, the Strong Allee effect, and a Hyper-Allee effect parameter regime with three distinct equilibria in (0, 1]. Additional Allee effects can be observed at the boundaries of the aforementioned Allee effects, including the Tangential Manifold and Junction Point with semi-stable equilibria, as well as the Triple Point with a single stable equilibria in (0, 1]. In all of these cases, there are between zero and three equilibria in the interval (0, 1].

Interpreting Experimental Data Using the Binary Switch Model

To demonstrate how the Binary Switch Model can be used to provide biological insight, we consider population-level data sets describing the growth of populations of cancer cells. Neufeld et al. (2017) perform three experiments with U87 glioblastoma cells. Uniform monolayers of cells are grown from three different initial densities, with the data shown in Fig. 8. Here, we see that all three experiments lead to increasing population densities with time. The two experiments with the smallest initial densities lead to increasing, concave up C(t) profiles. The experiment with the largest initial density leads to an increasing C(t) profile that changes concavity at approximately h.

Fig. 8

Population density of U87 glioblastoma cells compared to the calibrated Binary Switch Model. U87 glioblastoma cells, with initial densities of , , and , are observed over the span of 120 h (black circles) (Neufeld et al. 2017). The Binary Switch Model (solid curves) is fit to minimise the combined least-square error (15), , of three experimental data sets shown in Neufeld et al. (2017). The estimates of the optimal model parameter set, for each value of M, are shown in Table 7. b A semi-log plot makes it easier to visually compare the quality of match between the data and the model

Table 7

Estimates of the Binary Switch Model parameters that minimise the combined least-squares error (15) between model predictions and experimental data from Neufeld et al. (2017). The optimal parameter set with , highlighted in bold, provides the smallest combined least-squares error for all values of M

M	r	R	α	β	C_1(0)	C_2(0)	C_3(0)	Σχ^2
0	0.0113	0.0262	0.174	2.82×10^-6	0.0250	0.0661	0.184	0.0179
1	0.0168	0.0345	0.0608	0.0692	0.0192	0.0652	0.188	0.0154
2	0.0180	0.0576	2.84×10^-5	0.139	0.0160	0.0619	0.191	0.0169
3	0.0206	0.0642	3.66×10^-9	0.0892	0.0126	0.0534	0.193	0.0268
4	0.0218	0.134	3.43×10^-9	0.0623	0.0112	0.0489	0.191	0.0366
5	0.0237	0.0110	3.73×10^-10	2.34×10^-4	0.00933	0.0420	0.183	0.0571

The density of U87 glioblastoma cells has already been rescaled by its maximum packing density in Neufeld et al. (2017), so we assume that corresponds to the maximum rescaled density. Our aim is to choose , with and as initial conditions, such that the model parameters provide the best match to all three experimental conditions simultaneously. It is important to calibrate the model to match all three data sets simultaneously, because if (4) is consistent with the experimental data, there should be a single choice of model parameters that matches the observed population dynamics, regardless of initial density (Jin et al. 2016b).

To match all experimental data sets simultaneously, we consider the combined least-squares error between model predictions and all data:Here, we treat the initial densities, as parameters to be determined; therefore, we consider the extended parameter vector, . In (15), represents all three experimental data sets obtained at times , concatenated into a single vector. While the Binary Switch Model uses the initial conditions and , we denote the experimental measurements at h as and , respectively (Fig. 8). Using fminsearch in MATLAB (MathWorks 2020), we estimate such that is minimised. Since M is discrete, while are continuous, we estimate for each value of and then choose the value of M that minimises . A MATLAB implementation of this least-squares procedure is discussed in the Supplementary Information.

In Fig. 8, we show the best match that the Binary Switch Model can provide to all three data sets from Neufeld et al. (2017) for each value of M. The optimal parameter set and minimal for each value of M are reported in Table 7. We conclude that setting a threshold of provides the best match to these data sets. While larger values of M clearly deviate from the experimental data sets at low population densities (Fig. 8b), setting or also leads to a reasonable visual match for all three experimental data sets (Fig. 8). Furthermore, it is of interest to note that the optimal model parameters associated with small values of M correspond to non-logistic growth features, since logistic growth can only be obtained when and (Table 7). The match between the experimental data and the model at has several consequences: (i) this exercise confirms that the data reported by Neufeld et al. (2017) does not follow standard logistic growth; (ii) the high-quality match between the Binary Switch Model and the data for is consistent with population dynamics similar to a Weak Allee effect, and (iii) interpreting this data using the Binary Switch Model indicates that the best way to explain the population dynamics with a relatively small threshold population density.

Applications to Ecology

Threshold effects are thought to be a common feature in biological population dynamics, both in cell biology and in ecology. In the previous section, we demonstrated that a population of U87 glioblastoma cells did not follow logistic growth and was better described using the Weak Allee effect. In a similar fashion, various populations in ecology with known threshold effects are better described using Allee effects. A common threshold effect arising in ecology is a threshold population density (Courchamp et al. 2008), whereby a particular species will go extinct below this critical density. Species that have been noted to go extinct below a threshold density include the quokka (Sinclair and Pech 1996), the woodland caribou (Wittmer et al. 2005), the red-backed vole (Morris 2002), and the gypsy moth (Tcheslavskaia et al. 2002; Liebhold and Bascompte 2003). In many of these populations, the threshold density has been measured (Courchamp et al. 2008), thereby providing an appropriate estimate of the equilibrium density employed in the Binary Switch Model. Consequently, the Binary Switch Model aligns with threshold effects known to arise in ecology, while also providing insight into the underlying individual-level mechanisms that give rise to Allee effects. Furthermore, by making use of measured threshold population densities, we are thereby able to obtain an estimate of the threshold parameter M, further reducing the parameter search space needed to calibrate the Binary Switch Model to match experimental data.

Conclusions

In this work, we examine the link between threshold effects in population growth mechanisms and Allee effects. An abrupt change in growth mechanisms, which we refer to as a binary switch, is thought to be a common feature of biological population dynamics. Despite the ubiquitous nature of local binary switches in population dynamics, an explicit connection to Allee effects has not been considered. To explore this connection in greater detail, we examine a population density growth model, in which the proliferation and death rates vary with the local density of the population. By incorporating a local binary switch in these proliferation and death rates, we greatly reduce the size of the parameter space while explicitly incorporating a biologically realistic threshold effect in the proliferation and death rates.

To provide insight into the qualitative features of population dynamics arising in the Binary Switch Model, we examine the presence and stability of the resulting equilibria. We show that when the binary switch occurs at some intermediate population density and the high-density death rate is not too large, a diverse family of Allee effects is supported by the model. Among these Allee effects are: (i) logistic growth, when no binary switch is present; (ii) the Weak Allee effect, which modifies the simpler logistic growth model without changing its equilibria or their stability; (iii) an Extinction regime, where all population densities will eventually go extinct; (iv) the Strong Allee effect, where population below a critical density will go extinct rather than grow, and (v) the Hyper-Allee effect, which has two distinct positive stable population densities. Furthermore, we show that there are additional forms of Allee effects at the boundaries in the parameter space that separate these five main classes of Allee effects.

Along with exhibiting a wide range of Allee effects, the Binary Switch Model has a restricted parameter regime, making the interpretation of the local binary switch clearer while requiring fewer parameters to identify when calibrating to experimental data. To demonstrate these advantages, we calibrate the Binary Switch Model to experimental data sets arising in cell biology. Not only can the Binary Switch Model provide a good match to all experimental data across three different initial densities, we also find that the parameters used to match the data provide a more explicit interpretation of the underlying local growth mechanisms arising in the population. Specifically, we confirm that the experimental data are inconsistent with the standard logistic model and that the phenomena is best explained by a binary switch at low density. We conclude that the Binary Switch Model is useful to theorists and experimentalists alike in providing insight into binary switches at the individual scale that produce Allee effects at the population scale.

While one of the merits of the Binary Switch Model is to show how a single local binary switch gives rise to a variety of Allee effects, further extensions of the modelling framework can be made. For instance, additional switches can be incorporated into the modelling framework, representing populations whose proliferation and death rates change at more than one density. We anticipate that this kind of extension would lead to additional forms of Allee effects in the resulting population dynamics. Another potential modification would be to generalise the notion how we measure local density. In this work, we take the simplest possible approach use the number of nearest neighbours on a hexagonal lattice to represent the local density. Several generalisations, such as working with next nearest neighbours or working with a weighted average of nearest neighbours, could be incorporated into our modelling framework (Fadai et al. 2020; Jin et al. 2016a). Again, we expect that such extensions would lead to an even richer family of population dynamics models. We leave these extensions for future considerations.

Below is the link to the electronic supplementary material.

Supplementary material 1 (pdf 272 KB)