Global dynamics of a novel delayed logistic equation arising from cell biology

The delayed logistic equation (also known as Hutchinson's equation or Wright's equation) was originally introduced to explain oscillatory phenomena in ecological dynamics. While it motivated the development of a large number of mathematical tools in the study of nonlinear delay differential equations, it also received criticism from modellers because of the lack of a mechanistic biological derivation and interpretation. Here we propose a new delayed logistic equation, which has clear biological underpinning coming from cell population modelling. This nonlinear differential equation includes terms with discrete and distributed delays. The global dynamics is completely described, and it is proven that all feasible nontrivial solutions converge to the positive equilibrium. The main tools of the proof rely on persistence theory, comparison principles and an $L^2$-perturbation technique. Using local invariant manifolds, a unique heteroclinic orbit is constructed that connects the unstable zero and the stable positive equilibrium, and we show that these three complete orbits constitute the global attractor of the system. Despite global attractivity, the dynamics is not trivial as we can observe long-lasting transient oscillatory patterns of various shapes. We also discuss the biological implications of these findings and their relations to other logistic type models of growth with delays.


Introduction
The well known logistic differential equation N (t) = rN (t) (1 − N (t)/K) was proposed by Verhulst in 1838 to resolve the Malthusian dilemma of unbounded growth [2]. Here N (t) represents the population density at time t, r is referred to as the intrinsic growth rate, and K as the carrying capacity. Generally, such saturating population growth can be achieved by assuming that an increase in the population size leads to a decrease in fertility or an increase in mortality. This is reasonable when resources are limited: if the population size exceeds some critical level, the habitat cannot support further growth. The logistic equation has seen widespread use in modelling studies across biology with applications in single bacteria, cell and animal populations, and also interacting populations, cancer, and infectious diseases.
In the framework of the logistic ordinary differential equation, the population always converges to a stable steady state. However, in the context of ecology, oscillatory behaviours have also been observed, even in the absence of external periodic forcing. For example, in Daphnia populations, oscillations are possible because the fertility of females depends not merely on the actual population density, but also on the past densities to which it has been exposed. This motivated Hutchinson in 1948 [24] to propose the delayed logistic equation (also known as Hutchinson where we have normalized the carrying capacity to unity (by setting v(t) = N (t)/K), and introduced the time delay τ > 0. Interestingly, independent of Hutchinson, in 1955 Wright [39], motivated by an unpublished note of Lord Cherwell about a heuristic approach to the density of prime numbers, studied in detail the delay differential equation Wright's equation is an equivalent form of the delayed logistic equation, established via the change of variables u(t) = v(t) − 1.
Since the delayed logistic equation exhibits stable periodic solutions whenever ατ > π/2, it is very tempting to use this modification of the logistic ordinary differential equation to explain observed biological oscillations. One example is the work of May fitting the delayed logistic equation to the empirical data of oscillatory blowfly populations from Nicholson's experiments [32]. However, the use of Hutchinson's delayed logistic equation received criticisms from biological modellers due to the lack of a mechanistic biological derivation [14], the main objection being that it is not based on clearly defined birth and death terms, and the delay was inserted on a purely phenomenological basis. Indeed, in May's work the parameters in the mathematical model did not directly correspond to any biological parameters.
Later, Gurney et al. [18] proposed the equation Incorporating a delay in the growth term, Arino et al. [1] derived an alternative formulation of the delayed logistic equation which has recently been extended to distributed delays [29]. This model behaves similarly to the classical logistic equation, in the sense that it cannot sustain periodic oscillations, however large delays cause the extinction of the population. In another recent work, Lindström studied chemostat models with time lag in the conversion of the substrate into biomass [30], and obtained equations of logistic type with delays as a limiting case. Lindström's equation shares the same dynami-cal properties as the alternative delayed logistic equation from [1], both generating a monotone semiflow, unlike Hutchinson's equation.
In addition to these biological applications, the delayed logistic equation motivated the development of a large number of analytical and topological tools, including local and global Hopf-bifurcation analysis for delay differential equations [5,8,22,35], asymptotic analysis [13], 3/2-type stability criteria [25,39], and the study of slowly oscillatory solutions [28]. The most famous related problem is Wright's conjecture from 1955, which proposes global asymptotic stability in the delayed logistic equation for ατ ≤ π/2. This long-standing question has been resolved recently by combining analytical and computational tools [4,38]. Moreover, many variants and generalizations of the delayed logistic equation have been studied in the mathematical literature, for instance equation including positive instantaneous feedback [19], non-autonomous terms [8], neutral terms [16], spatial diffusion [41], and multiple delays [40]. Further examples are discussed in [27,31,36].
In this paper we derive a novel type of delayed logistic equation, inspired by cell biology. Our equation includes both terms with and without delays, as well as with distributed delays. After presenting derivation of the model in Section 2, we study its basic mathematical properties in Section 3, such as well-posedness, existence and stability of equilibria. The global dynamics is fully described in Section 4, and we provide a complete description of the global attractor. While all non-trivial solutions converge to the positive equilibrium, the model exhibits long-lasting transient oscillations of various shapes; this phenomenon is discussed in Section 5. Finally, we summarize our findings and discuss their implications for future research in Section 6.

Derivation of the new delayed logistic equation
We derive our new delayed logistic equation by considering the mean-field limit of a simplistic agent-based, on-lattice model that describes both cell motility and proliferation, and cell-cell interactions through the incorporation of volume exclusion.
More specifically, we assume that agents move and proliferate on an n-dimensional square lattice with spacing ∆ and length (in each direction) ∆, where is an integer describing the number of lattice sites. A great interest in current medical research is to understand the roles of two particular phenotypes, namely "high proliferation-low migration" and "low proliferation-high migration", in the progression of agressive cancers [34]. Our model is designed to capture the "go or grow" hypothesis frequently acknowledged in the cancer literature [12,15], which proposes that cell proliferation and migration are temporally exclusive events.
As such, we divide our agent population into two subpopulations, motile and proliferative. Each agent is assigned to a lattice site, from which it can move or proliferate into an adjacent site. If a motile agent attempts to move into a site that is already occupied, the movement event is aborted. Similarly, if a proliferative agent attempts to proliferate into a site that is already occupied, the proliferation event is aborted. Agents convert from being motile to proliferative at constant rate r > 0 per unit time. That is, rδt is the probability that a motile agent switches to become proliferative in the next infinitesimally small time interval δt. Upon switching, agents remain immobile for a fixed time τ > 0 (representing the length of time taken for cells to progress through the cell cycle to division), and they then attempt to proliferate by placing a daughter agent into one of the nearest neighbour lattice sites. After attempting proliferation, proliferative agents switch back to being motile. The initial agent distribution is, on average, spatially uniform, and is achieved by populating lattice sites uniformly at random with probability Ps.
Following arguments presented in Baker and Simpson [3], and closing the hierarchy of moment equations at the pair level (assuming that the occupancy of neighbouring lattice sites is independent), we can derive a system of delay differential equations that describe how the density of agents on the lattice evolves over time: where m is the density of motile agents, p is the density of proliferative agents, K > 0 is the number of sites on the lattice (one can think of the carrying capacity as K = n ), and denotes differentiation with respect to t. The equations encode the model assumption that motile agents become proliferative with rate r, and stay in this proliferative phase for time τ . Once the cell cycle is completed, proliferative cells switch back to being motile cells again. As they switch, they attempt to place a daughter agent into a lattice site; the probability that site is empty (in the System (5)-(6) can be rewritten as a scalar equation using the natural biological relation that is, the proliferative cells at time t are exactly those who entered the prolifera- Using this relation in equation (5) gives Next we rescale the density with K by lettingm(t) = m(t)/K, and drop the tildes and rearrange to To complete the non-dimensionalisation, we also rescale time by lettingt = t/τ and x(t) = m(t). Then, dx/dt = τ m (t) and, moreover, we can, again, drop the tildes, to arrive at the new delayed logistic equation: 3 Basic mathematical properties With the notation ρ = rτ > 0, we now analyse the scalar differential equation with discrete and distributed delays  (13) is Then the map f : C → R is given by and the initial data is specified by

Well-posedness
By a solution of equations (13), (16) we mean a continuous function x(t) on an interval [−1, A) with 0 < A ≤ ∞, which is differentiable on (0, A), satisfies equation (13) on (0, A) and also satisfies equation (16). Proof Recall the standard existence and uniqueness theorem for functional differential equations [37]. We shall verify that the local Lipschitz property holds for f , i.e. for any M > 0, there is an L > 0 such that For a constant M > 0, let ||φ|| < M and ||ψ|| < M , then we have the estimates Due to biological constraints, we are interested only in non-negative solutions, i.e.
This set in fact corresponds to the biologically feasible phase space since, for a solution with x t ∈ Ω, represents the total cell density (accounting for all mobile and proliferating cells) that should, in the non-dimensional model, lie between zero and unity. We shall use Ω ⊂ C, the set of biologically feasible states as our phase space, noting that Ω depends on the parameter ρ.

Positivity and boundedness from above
which simplifies to whenever w(0) ≥ 0, or equivalently θ(0) ≤ 1. We find that holds if The general positivity condition, i.e. that x t ≥ 0 whenever x 0 ≥ 0, for system (13), (16) is Theorem 3.4]. This clearly holds as for that Ω is positively invariant, and additionally the estimate x(t) < 1 holds.
Consequently, solutions starting from φ ∈ Ω exist globally, i.e. on the interval

Steady states and their stability
Theorem 3.1 Equation (13) has two steady states, x 0 = 0 which is always unstable, and x * = 1/(ρ + 1) which is always locally asymptotically stable.
Linearisation of equation (13) around x = 0 gives with the characteristic equation which always has a positive real root, hence 0 is unstable.
Next we look at the stability of the equilibrium x * . Using the definition in equation (15), we have which has linear part and Clearly lim ||φ||→0 |g(φ)|/||φ|| = 0, hence the linear variational equation is The exponential ansatz y(t) = e λt gives the characteristic equation which becomes, after integration and rearranging, provided λ = 0. Note that λ = 0 is not a root of equation (33), otherwise which is a contradiction. Therefore we may write the characteristic equation as with We can factor the characteristic equation as which demonstrates that we always have the real root λ = −ρ/(ρ + 1) < 0. The other roots are the roots of λ + ρ − ρe −λ = 0. This is a well-known type of characteristic equation of the form λ = A + Be −τ λ , see for example Chapter 4.5 of [37].
However, we are on the stability boundary of this equation given that λ = 0 is a root (which is, as we established, not a root of the original characteristic equation (33)). We can quickly show that all other roots have imaginary parts less than zero: assuming a root with λ ≥ 0 but λ = 0, we have which is a contradiction. Therefore the equilibrium x * is locally asymptotically stable.
4 Global dynamics

Persistence
Theorem 4.1 Equation (13) is strongly uniformly persistent, that is there exists a δ > 0, independent of the initial data, such that for any solution Proof Let then equation (13) can be written as From the variation of constants formula, for any σ > ω > 0 we find We use the notation x∞ = lim inf t→∞ x(t), and x φ to denote a specific solution with initial function φ. Assume the contrary of the statement of the theorem.
Let the following assumptions be satisfied: (i) λ = a 0 + b 0 e −λτ has a unique root λ 0 with largest real part, and this root λ 0 is real and simple; Then, every solution x(t) satisfies where and ξ = ξ(x) is a constant depending on the solution x.
Fix some ψ ∈ Ω, and let x ψ (t) be the corresponding solution with where the last equality was given in equation (23). Now consider the equation We apply Theorem A, and check all the conditions, where we have a 0 = −ρ, b 0 = ρ, τ = 1, a(t) = 0 and b(t) = ρw ψ (t). The conditions with a(t) are trivial, and the conditions with b(t) follow from the combination of persistence and equation (57) as follows, since w is exponentially convergent.
(ii) To see that w ψ (t) → 0 as t → ∞, we can use Theorem 4.1: there is a δ > 0 and a T such that x(t) > δ for t > T − 1. From equation (57) we find (iii) The estimate in equation (60) shows that (61) (iv) From the estimate in equation (60), we also find that As a result, and all the conditions of Theorem 4.3 hold. Thus, every solution z(t) satisfies Now notice that a solution x(t) of equation (13) with initial function ψ is also a solution of equation (58) There is a γ > 0 such that λ 0 > γ > λ j for all j ≥ 1, and there is a k ∈ N such that λ k > 0 and λ k+1 ≤ 0, where the point (−ρ, 2ρ) lies between the curves C − k and C − k+1 on the (α, β)-plane, given by The intersections of C − j and the half-line (−ρ, 2ρ) (ρ > 0) are determined by cos ν = 1/2, thus ν = 2jπ − π/3, sin ν = − √ 3/2 and the intersection is given by Hence, k is either zero or the largest positive integer such that ρ > ρ k .
Considering the leading real root, the corresponding eigenfunction is given by where the function χ 0 ∈ C spans the linear eigenspace P 0 := {cχ 0 : c ∈ R}. P k is a 2k-dimensional eigenspace corresponding to {λ j : j = 1, . . . , k}. If λ 0 is the only eigenvalue with positive real part then P k = ∅. There is an ε > 0 such that λ k > ε. Q corresponds to the remaining part of the spectrum.
The unstable set W 0 (0) intersects Ω, and for any function φc of this intersection, there is a complete solution x(t) : R → R + 0 such that x 0 = φc and x t → 0 * as t → −∞. By the negative invariance of W 0 (0), the trajectory is completely in Ω, hence we have proved the existence of a heteroclinic orbit connecting the two equilibria.
For uniqueness, consider the unstable manifold Wu(0). Any solution on the set Wu(0)\W 0 (0) oscillates about 0 as t → −∞, hence can not be a positive heteroclinic solution. Denote the unique heteroclinic orbit in Ω by H. We claim that H and the two equilibria constitute the global attractor, otherwise there are further complete orbits in Ω. Assume there is a ψ ∈ Ω which is on such a complete orbit Γ ⊂ Ω, so ψ = 0, ψ = x * , and ψ / ∈ H. Then its α-limit set α(ψ) and ω-limit set ω(ψ) are non-empty, compact and invariant. In particular, by Theorem 4.3, ω(ψ) = {x * }. If α(ψ) = {x 0 }, then Γ is a connecting orbit coinciding with H. Hence, we may assume that there is a φ = 0 such that φ ∈ α(ψ). Then x φ t → x * as t → ∞, so x * ∈ α(ψ) as well. But this contradicts the stability of x * . Thus, the global attractor does not contain any other orbit, and A = {x 0 } ∪ {x * } ∪ H.

Numerical simulations and metastability
Representative numerical simulations of solutions of equation (13)  2, we can see that the solution is in fact not periodic, however the convergence is very slow. This is a situation that is sometimes called metastability. Paraphrasing [23], metastability refers to a situation when something appears not to change on short time scales, while it changes after a sufficiently long period of time.
Such phenomena are common in boundary value problems of partial differential equations, and metastability has been observed for delay differential equations as well, see, for example, [7,17]. Long lasting transient oscillations have been reported in [33] for a scalar differential equation with a constant delay. It is interesting to see that equation (13) can exhibit rapid convergence for small ρ and oscillatory patterns of various shapes for larger ρ, as illustrated in Fig. 3 where solutions are plotted with different initial functions. To gain some understanding of this phenomenon, we look at our model from the point of view of singular perturbations. In equation (13), let us divide by ρ 2 .
For large ρ, we use the notation ε = 1/ρ, then we have Setting ε = 0, we obtain the singular perturbation In addition, we can look at the singular perturbation of the characteristic equation: equation (39) can be written as Multiplying by (ρ + 1)/ρ 2 we have which is with the notation ε = 1/ρ. Setting ε = 0, we obtain the singular perturbation The roots of this equation satisfy either λ = −1 or e −λ = 1, having the purely imaginary roots i = 2kπ, k ∈ Z. Indeed, as we increase ρ, we can see that all complex roots line up along the imaginary axis, see Fig. 4.

Discussion
The delayed logistic equation has received much attention in past decades in the analysis of nonlinear delay differential equations. However, its biological validity We have also fully described the global attractor as the union of the two equilibria and a unique connecting orbit. This heteroclinic solution is depicted in Fig.   5 when ρ = 20. In this case, the leading real eigenvalue at zero is λ ≈ 0.66, while the leading pair of eigenvalues at the positive equilibrium are λ ≈ −0.04 ± 6i. We can numerically observe how the solution is aligned to the leading eigenspaces of the linearizations near the two equilibria.
Our model is based on the commonly-invoked modelling assumption that proliferating cells abort proliferation when placement of a daughter cell is not possible due to spatial crowding constraints. However, there are other potential models of crowding-limited proliferation that could be encoded within the same framework.
For example, one could assume that, instead of aborting the proliferation attempt, proliferating cells instead enter a waiting state until space for the daughter cell becomes available, or that, in order to enter the proliferating state, cells must be able to "reserve" a site for the daughter cell to be placed into. These, and other possible biological hypotheses constitute a family of different logistic-type models, the behaviours of which we will systematically compare in subsequent works.
As far as we are aware, this work represents the first step to understand the go-or-grow mechanism when a delay caused by the cell cycle length is explicitly incorporated as a biological parameter. Future work will explore more accurate models of the effects of spatial crowding upon proliferation, using both mean-field and moment dynamics models. In the transformed equation (13), the parameter ρ is the product of the time delay, τ , and the proliferation rate, r. As such, increasing either of those two parameters will have a similar effect on the dynamics (although on different time scales in terms of the original, dimensional equation). In models that include spatial details, however, there is an additional key parameter (that is ignored by equation (13)), namely the cell motility rate. An important future goal is to understand how the speed of e.g. cancer invasion depends on cell-level parameters such as cell motility and proliferation rates, and the cell cycle length.
To this end, spatial models of the go-or-grow mechanism that incorporate cell cycle delays need to be developed. Travelling wave solutions of reaction-diffusion systems with delays are closely related to the heteroclinic orbits of the reaction systems (see for example [11]). As such, the results we provide on the unique heteroclinic orbit connecting the zero and the positive equilibria will help shed light on invasion speeds in these cases.