Unpacking the Allee effect: determining individual-level mechanisms that drive global population dynamics

(a) The individual-based model

We perform non-dimensional simulations where the hexagonal lattice spacing is Δ. These simulations can be rescaled to match any particular application by rescaling Δ [23,28]. Each lattice site has position

with i = 1, …, I, j = 1, …, J and Δ = 1.

Crowding effects are important in both cell biology and ecology [10,23,24,28–30], so potential motility and proliferation events that would result in more than one agent per site are aborted. Agents attempt to undergo nearest neighbour motility events at rate m_n ≥ 0, proliferation at rate p_n ≥ 0, and death events at rate d_n ≥ 0 (figure 2c,d). Here, $n\in \{0,1,\dots ,6\}$ is the number of occupied nearest neighbour sites, providing a simple measure of the local density. While individuals attempt to undergo motility and proliferation events at a constant rate, the actual rate of successful events is density dependent, with the net rates being decreasing functions of density. This is the key feature of the model that gives rise to non-logistic phenomena.

For simplicity, we assume that all agents, regardless of their local density, have the same motility rate (m_n ≡ m). Furthermore, we choose m such that p_n/m ≪ 1 and d_n/m ≪ 1 for all n, since the characteristic timescale for motility is assumed to be much shorter than the characteristic timescale for proliferation and death [31]. Agents are initially seeded on the lattice with a constant probability, representing spatially uniform initial conditions. Furthermore, we impose reflecting boundary conditions and, using a Gillespie approach [32], we simulate the number of agents as a function of time and space (electronic supplementary material, Algorithm 1).

To compare data from the IBM with the global population description, we average data from the IBM using

Here, Q_ℓ(t) is the total number of agents on the lattice at time t, in the ℓth identically prepared realization of the IBM. The total number of identically prepared realizations is L; we choose L = 100 for the results presented in figures 3 and 4. A description of the numerical algorithm (electronic supplementary material) and a Matlab implementation of this algorithm are available (Code Availability).

• Download figure
• Open in new tab
• Download PowerPoint

• Download figure
• Open in new tab
• Download PowerPoint

(b) Continuum limit

While the IBM allows us to visualize realistic-looking individual simulations of population dynamics, as well as to explicitly specify the individual-level behaviour of agents, it is convenient to derive a simpler mathematical description of the average behaviour of the IBM, called the continuum limit description [10,28]. The continuum limit description gives us the ability to study global, deterministic, features of the IBM when the number of lattice sites is large, as well as to understand how individual-level differences translate into global outcomes.

Since the IBM employs spatially uniform initial conditions, the net flux of agents entering and leaving each lattice site due to motility events is, on average, zero. Therefore, spatial derivatives in the continuum limit will vanish, meaning that the continuum description of the average agent density, 0 ≤ C ≤ 1, is a function of time alone. Furthermore, we assume that the occupancy status of lattice sites is independent. This assumption, called the mean-field approximation [10,23,28,33], is mathematically convenient and is consistent with setting m ≫ p_n and m ≫ d_n [23,33]. Our results confirm that the mean-field approximation is very accurate for populations that migrate faster than they proliferate or die. Populations with slow migration rates, such as populations considered in plant ecology, are not expected to agree with the mean-field assumption continuum limit, as local clustering will develop and thereby violate the assumption of spatial uniformity [23,33].

Using the mean-field approximation, we describe the local density of agents in terms of the average agent density. Specifically, for an agent to have a local density corresponding to n nearest neighbours, we require: (i) an agent to be present at the particular site; (ii) n nearest neighbour sites are occupied and; (iii) the remaining (6 − n) nearest neighbour sites are vacant. These conditions give the density of agents with n nearest neighbours, 0 ≤ I_n ≤ 1, as

hence, I_n follows a binomial distribution. To determine how the global agent density evolves, we examine the time rate of change in expected agent density due to proliferation and death events, giving where n(i) is the number of nearest neighbours is each neighbouring agent i. The binomial coefficient, $\left(\begin{array}{c}6\\ j\end{array}\right)=6!/[j!(6-j)!]$, accounts for all possible configurations of an agent’s j nearest neighbours. Furthermore, by assuming that the agent configurations I_n follow a binomial distribution, we can write We omit the saturated configuration, I₆ (black hexagon in figure 2), from this sum, since proliferation events require an empty neighbouring lattice site to produce a daughter agent.

By combining equations (2.4)–(2.5), we obtain the continuum limit description (i.e. the global population description) of C(t):

For convenience, we quantify the death rate as a fraction of the proliferation rate, i.e.d_n = α_n p_n, where α_n ≥ 0. With this substitution and some rearranging of terms in equation (2.6), we can write the evolution of the global agent density in terms of the per capita growth rate, (1/C)(dC/dt): where

We deliberately refer to the per capita growth rate associated with the continuum limit as g(C), and the per capita growth rate prescribed via the global population behaviour as f(C). Later, we will show how to choose the IBM parameters such that f(C) ≡ g(C). Additionally, we note that equation (2.7) groups the 14 parameters (p_n, α_n) as 7 linearly independent parameters, γ_n.

(c) Allee effects arising from the individual-based model: the forward problem

We now demonstrate that the IBM framework gives rise to a rich variety of Allee effects. In the first instance, we set α_n = 0 for all n. Equation (2.8) implies that γ_n = p_n for 0 ≤ n ≤ 5 and γ₆ = 0, such that agents only proliferate and move, and do not die. While this parameter regime is by no means a complete account of all possible parameter combinations, it serves to highlight that g(C), from equation (2.7), can give rise to a suite of Allee effects.

We consider three different choices of p_n:

We note that in cases 2 and 3, the range of proliferation rates varies by a factor of 10 as n varies from zero to five, demonstrating a significant range of proliferation rates within a single IBM realization. Equation (2.7) gives,

The global population density, C, can be determined by solving equation (2.7) with a specified initial condition, C(0).

While we can retrieve logistic growth from the continuum limit of the IBM (case 1; figure 3), we can also obtain more complicated, nuanced, per capita growth rates, including the Weak Allee effect (case 2; figure 3) and completely novel Weak Allee-like per capita growth rates never previously described (case 3; figure 3). This analysis of three simple IBM parameter regimes suffices to show that this IBM framework is related to a large class of global population descriptions. Furthermore, in figure 3, as well as videos presented in the electronic supplementary material, we show that averaged density data from the IBM, 〈C〉, agrees very well with C, confirming that the average agent density from the simulated IBM data is faithfully captured by the global population description. With these results, we now turn to the inverse problem to determine which individual-level rates describe global Allee effects.

(d) Choosing individual-based model rates to match Allee effect models: the inverse problem

A standard practice for the application of an Allee-type model involves matching population-level experimental data to a particular continuum model, without any regard for the underlying individual-level mechanisms. Therefore, it is natural to ask, for a given per capita growth rate, f(C), how do we choose the IBM parameters so that the per capita growth rate determined from the continuum limit of the IBM, g(C) in equation (2.7), is identically f(C)?

Mathematically, we seek the parameters γ_n in equation (2.7) such that g(C) ≡ f(C). To do this, we first note that equation (2.7) is linear in γ_n. Therefore, we can evaluate equation (2.7) at seven distinct points $C_i\in \{C_0,C_1,\dots ,C_6\}$ and obtain a corresponding linear system in γ_n. By denoting

we obtain the linear system where the 7 × 7 matrix $\mathbb{M}$ has entries [M_ij], $\gamma \gamma$ is the vector of parameter values [γ_j] and f is the vector of function values [f(C_i)]. We note that the polynomials appearing in (2.7), Cⁿ(1 − C)⁶⁻ⁿ, are linearly independent, so $\mathbb{M}$ is never singular and the solution of this system is $\gamma \gamma ={\mathbb{M}}^{-1}\mathbf{f}$. Furthermore, if f(C) is a polynomial of degree 6 or less and f ≠ 0, this unique solution of $\gamma \gamma$ provides the IBM parameters such that g(C) ≡ f(C). We discuss a method of determining the IBM parameters for other forms of f(C) in the electronic supplementary material.

To prevent populations from becoming infinite, we require that f(C) will be non-positive when C = 1, so we impose f(1) ≤ 0. Noting that g(1) = −γ₆ = −p₆α₆ ≤ 0, this implies that γ₆ = −f(1) ≥ 0, which is consistent with the restriction that p₆, α₆ ≥ 0. Common per capita rates and the corresponding IBM parameters are tabulated in the electronic supplementary material.

Upon determining suitable parameters γ_n that match a specific f(C), we then consider choosing p_n and α_n from equation (2.8). However, since the birth and death rates (p_n, α_n) need to be determined solely from γ_n, there are infinitely many different combinations of (p_n, α_n) that satisfy equation (2.8). If there is further information about either p_n or α_n, such as additional experimental data, a unique estimate of (p_n, α_n) can be determined. However, in the absence of such data, we make a straightforward choice of parameter combinations

Here, R = max _1≤n≤6 p_n and only appears in the case when f(0) < 0. This choice of R provides a balance between minimizing the relative death rate α₀ > 1 while preventing p₀ from dominating other proliferation rates. These parameters can be interpreted as proliferation-dominated behaviour when γ_n > 0 and death-dominated behaviour when γ_n < 0. Furthermore, having γ₆ > 0 only occurs when f(1) < 0, implying that the death rate d₆ = p₆α₆, shown in black in figure 2, must be strictly positive.

(e) Mechanistic interpretation of experimental data

Per capita growth rates describing population growth and extinction can match the global trends in experimental data [2–7], but fail to provide any insight at the individual-level scale. By contrast, the mathematical tools presented in this work provide the missing link that connects various global per capita rates to specific individual-level mechanisms, via solving equation (2.12). Therefore, we can provide insight into individual-level behaviour from experimental data with the following approach. First, we choose a per capita growth rate to match the global features of the experimental data. Second, the associated global model parameters are then fit to the data: for example, by minimizing the least-square error between the population model and the experimental data. Last, we solve the inverse problem and determine the associated IBM parameters that give rise to the experimentally observed global behaviour.

To highlight the insight possible from this approach, we consider two population-level datasets and provide previously hidden detail about individual-level behaviour [6,7]. In Johnson et al. [6], BT-474 breast cancer cells are seeded at three initial densities in a 96-well plate, and cell proliferation is observed for 328 h. Similar epithelial cell lines $20–30\hspace{0.17em}\text{μ}\text{m}$ in diameter have motility rates of about $4–9\hspace{0.17em}{\text{h}}^{-1}$, while their proliferation rates are approximately $0.03–0.06\hspace{0.17em}{\text{h}}^{-1}$ [31]. Therefore, we assume that these BT-474 breast cancer cells have motility rates that are sufficiently large compared to both proliferation and death rates. As shown in the electronic supplementary material, a motility rate that is ten times larger than the proliferation rate is sufficient for the well-mixed assumption of the IBM continuum limit to hold.

Model selection analysis in Johnson et al. [6] suggests that an Allee effect is required to describe this data. As there is no evidence of population extinction, we consider the Weak Allee effect with per capita rate

where r is the per capita growth rate at C = 0 and K is the carrying capacity with units cells/image. We match this model to three experimental datasets simultaneously by minimizing the total least-squares error (figure 5) and we determine the underlying individual-level mechanisms associated with equation (2.14) using equation (2.12). Specifically, the global Weak Allee parameters are determined to be A = 148 cells/image, K = 315 cells/image and r = 0.00757 h⁻¹. In figure 5b, we show that the IBM proliferation rates associated with this Weak Allee effect increase linearly with local density, as we observed when discussing the forward problem (figure 3). This is to be expected, since the Weak Allee effect also features a reduced growth rate at low densities. Furthermore, the individual-level rates increase by a factor of about three, which is significantly different to simpler models such as the logistic growth model, where the proliferation rates are independent of local density (figure 4).

• Download figure
• Open in new tab
• Download PowerPoint

As carrying capacity densities are often reported in units of cells/well, we note that K = 315 cells/image is equivalent to K = 33 500 cells/well. To convert K from units of cells/image to units of cells/well, we make use of two key measurements: the size of the well, and the size of the image containing cells. The Trueline 96-well plate [6] has a well diameter of 6400 μm, while the largest region of an image shown in [6] that contains cells is approximately a 550 μm × 550 μm square. Therefore, we have that

It is important to note that K = 33 500 cells/well is simply the predicted carrying capacity density based on these low-cell number experiments at early time. Without any late time experimental data, where the density approaches the carrying capacity, it is impossible to predict the carrying capacity [35]. Nevertheless, we infer the carrying capacity that minimizes the least-squares error between model predictions and early-time experimental data. We conclude that the experimental observations presented in Johnson et al. [6] can be explained by allowing the attempted proliferation rates of BT-474 breast cancer cells to linearly increase as a function of local density. Without posing and solving the inverse problem described here, this individual-level insight is not possible.

To compare continuum descriptions with experimental results that arise in ecology, we consider the datasets in Melica et al. [7], where the dynamics of the population density of Aurelia aurita polyps growing on oyster shells is measured. New Aurelia aurita polyps form via free-swimming propagules, which are released from the existing polyp into the surrounding water and land elsewhere on the oyster shell [36]. Consequently, these polyps experience a faster motility rate, relative to proliferation and death rates, over the surface of the oyster shell. As such, the mean-field assumption inherent in our modelling framework is satisfied. Melica et al. [7] observe that a series of experiments where polyps are initially distributed at low density leads the population density evolve to some particular carrying capacity. Interestingly, when the same experiments are performed at a much higher initial density, the population decays to a different, higher, carrying capacity. This observation is explained in Melica et al. [7] by supposing the population dynamics are logistic, with the requirement that there are two different carrying capacities, despite the fact that the classic logistic model specifies a single carrying capacity only. Instead, we consider the ‘Hyper-Allee effect’ model, figure 4d, with per capita rate

Equation (2.15) has the advantage that two stable equilibria, A and K, exist. While other cubic per capita growth rates have been proposed [37], this work focuses on the underlying individual-level mechanisms giving rise to such a model. We fit this model to two experimental datasets simultaneously by minimizing the combined least-squares error of both datasets. As seen in figure 6, the Hyper-Allee effect model agrees with both experimental datasets with only a change in the initial condition.

• Download figure
• Open in new tab
• Download PowerPoint

To provide further insight about the population dynamics of this experiment, we determine the corresponding IBM parameters that give rise to this Hyper-Allee effect. We note that if f(C) is rescaled such that the largest recorded density point, the high-density initial condition, corresponds to C = 1 (electronic supplementary material), then we must have f(1) < 0. Equation (2.12) implies that γ₆ > 0, producing non-zero attempted proliferation and death rates when an individual has n = 6 nearest neighbours, shown in black bars in figure 6b,c. However, the presence of death rates at various local densities (figure 6c) drives low-density populations to the smaller carrying capacity A, while high density populations are driven to the larger carrying capacity K. We conclude that the experimental observations in Melica et al. [7], displaying the coexistence of two stable population densities, can be explained by varying the death rates with local density. In both examples, there is a clear need to use global population models more nuanced than logistic growth to capture the experimentally observed features.