Global dynamics of a novel delayed logistic equation arising from cell biology
NNoname manuscript No. (will be inserted by the editor)
Global dynamics of a novel delayed logistic equationarising from cell biology
Ruth E. Baker · Gergely R¨ost
January 24, 2019
Abstract
The delayed logistic equation (also known as Hutchinson’s equation orWright’s equation) was originally introduced to explain oscillatory phenomena inecological dynamics. While it motivated the development of a large number ofmathematical tools in the study of nonlinear delay differential equations, it alsoreceived criticism from modellers because of the lack of a mechanistic biologicalderivation 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 distributeddelays. The global dynamics is completely described, and it is proven that all feasi-ble nontrivial solutions converge to the positive equilibrium. The main tools of theproof rely on persistence theory, comparison principles and an L -perturbationtechnique. Using local invariant manifolds, a unique heteroclinic orbit is con-structed that connects the unstable zero and the stable positive equilibrium, andwe show that these three complete orbits constitute the global attractor of thesystem. Despite global attractivity, the dynamics is not trivial as we can observelong-lasting transient oscillatory patterns of various shapes. We also discuss thebiological implications of these findings and their relations to other logistic typemodels of growth with delays. Keywords : cell population, logistic growth, go or grow models, time delay,global attractor, stability, long transient.
Mathematical Subject Classfication:
Ruth E. BakerMathematical Institute, University of Oxford, UKGergely R¨ostMathematical Institute, University of Oxford, UK& Bolyai Institute, University of Szeged, HungaryE-mail: [email protected] a r X i v : . [ m a t h . D S ] J a n Ruth E. Baker, Gergely R¨ost
The well known logistic differential equation N (cid:48) ( t ) = rN ( t ) (1 − N ( t ) /K ) was pro-posed 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 theintrinsic growth rate, and K as the carrying capacity. Generally, such saturatingpopulation growth can be achieved by assuming that an increase in the populationsize leads to a decrease in fertility or an increase in mortality. This is reasonablewhen resources are limited: if the population size exceeds some critical level, thehabitat cannot support further growth. The logistic equation has seen widespreaduse in modelling studies across biology with applications in single bacteria, celland animal populations, and also interacting populations, cancer, and infectiousdiseases.In the framework of the logistic ordinary differential equation, the populationalways converges to a stable steady state. However, in the context of ecology, oscil-latory behaviours have also been observed, even in the absence of external periodicforcing. For example, in Daphnia populations, oscillations are possible because thefertility of females depends not merely on the actual population density, but alsoon the past densities to which it has been exposed. This motivated Hutchinsonin 1948 [24] to propose the delayed logistic equation (also known as Hutchinson’sequation) v (cid:48) ( t ) = αv ( t ) [1 − v ( t − τ )] , (1)where we have normalized the carrying capacity to unity (by setting v ( t ) = N ( t ) /K ), and introduced the time delay τ >
0. Interestingly, independent ofHutchinson, in 1955 Wright [39], motivated by an unpublished note of Lord Cher-well about a heuristic approach to the density of prime numbers, studied in detailthe delay differential equation u (cid:48) ( t ) = − αu ( t − τ ) [1 + u ( t )] . (2)Wright’s equation is an equivalent form of the delayed logistic equation, establishedvia the change of variables u ( t ) = v ( t ) − lobal dynamics of a novel delayed logistic equation arising from cell biology 3 Since the delayed logistic equation exhibits stable periodic solutions whenever ατ > π/
2, it is very tempting to use this modification of the logistic ordinarydifferential equation to explain observed biological oscillations. One example isthe work of May fitting the delayed logistic equation to the empirical data ofoscillatory blowfly populations from Nicholson’s experiments [32]. However, theuse of Hutchinson’s delayed logistic equation received criticisms from biologicalmodellers due to the lack of a mechanistic biological derivation [14], the mainobjection being that it is not based on clearly defined birth and death terms, andthe delay was inserted on a purely phenomenological basis. Indeed, in May’s workthe parameters in the mathematical model did not directly correspond to anybiological parameters.Later, Gurney et al. [18] proposed the equation N (cid:48) ( t ) = aN ( t − τ ) e − bN ( t − τ ) − µN ( t ) , (3)which is known today as Nicholson’s blowfly equation. It is structurally differentfrom the delayed logistic equation, and the terms have clear biological interpreta-tion as the birth and mortality rates, while the delay represents the period betweenegg laying and hatching. Furthermore, in contrast to the delayed logistic equation,Nicholson’s blowfly equation can explain not only the appearance of cycling popu-lations, but also the two bursts of reproductive activity per adult population cyclethat appeared in Nicholson’s data.Incorporating a delay in the growth term, Arino et al. [1] derived an alternativeformulation of the delayed logistic equation N (cid:48) ( t ) = γµN ( t − τ ) µe µτ + κ ( e µτ − N ( t − τ ) − µN ( t ) − κN ( t ) , (4)which has recently been extended to distributed delays [29]. This model behavessimilarly to the classical logistic equation, in the sense that it cannot sustain pe- riodic oscillations, however large delays cause the extinction of the population. Inanother recent work, Lindstr¨om studied chemostat models with time lag in theconversion of the substrate into biomass [30], and obtained equations of logistictype with delays as a limiting case. Lindstr¨om’s equation shares the same dynami- Ruth E. Baker, Gergely R¨ost cal properties as the alternative delayed logistic equation from [1], both generatinga monotone semiflow, unlike Hutchinson’s equation.In addition to these biological applications, the delayed logistic equation mo-tivated the development of a large number of analytical and topological tools,including local and global Hopf-bifurcation analysis for delay differential equa-tions [5,8,22,35], asymptotic analysis [13], 3/2-type stability criteria [25,39], andthe study of slowly oscillatory solutions [28]. The most famous related problemis Wright’s conjecture from 1955, which proposes global asymptotic stability inthe delayed logistic equation for ατ ≤ π/
2. This long-standing question has beenresolved recently by combining analytical and computational tools [4,38]. More-over, many variants and generalizations of the delayed logistic equation have beenstudied in the mathematical literature, for instance equation including positiveinstantaneous feedback [19], non-autonomous terms [8], neutral terms [16], spatialdiffusion [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 bycell biology. Our equation includes both terms with and without delays, as wellas 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 inSection 4, and we provide a complete description of the global attractor. Whileall non-trivial solutions converge to the positive equilibrium, the model exhibitslong-lasting transient oscillations of various shapes; this phenomenon is discussedin Section 5. Finally, we summarize our findings and discuss their implications forfuture research in Section 6.
We derive our new delayed logistic equation by considering the mean-field limit ofa simplistic agent-based, on-lattice model that describes both cell motility and pro-liferation, and cell-cell interactions through the incorporation of volume exclusion.More specifically, we assume that agents move and proliferate on an n -dimensional lobal dynamics of a novel delayed logistic equation arising from cell biology 5 square lattice with spacing ∆ and length (in each direction) (cid:96)∆ , where (cid:96) is an in-teger describing the number of lattice sites. A great interest in current medicalresearch is to understand the roles of two particular phenotypes, namely “highproliferation-low migration” and “low proliferation-high migration”, in the pro-gression of agressive cancers [34]. Our model is designed to capture the “go orgrow” hypothesis frequently acknowledged in the cancer literature [12,15], whichproposes that cell proliferation and migration are temporally exclusive events.As such, we divide our agent population into two subpopulations, motile andproliferative. Each agent is assigned to a lattice site, from which it can move orproliferate into an adjacent site. If a motile agent attempts to move into a site thatis already occupied, the movement event is aborted. Similarly, if a proliferativeagent attempts to proliferate into a site that is already occupied, the proliferationevent is aborted. Agents convert from being motile to proliferative at constant rate r > rδt is the probability that a motile agent switchesto become proliferative in the next infinitesimally small time interval δt . Uponswitching, agents remain immobile for a fixed time τ > P s .Following arguments presented in Baker and Simpson [3], and closing the hi-erarchy of moment equations at the pair level (assuming that the occupancy ofneighbouring lattice sites is independent), we can derive a system of delay differ-ential equations that describe how the density of agents on the lattice evolves overtime: m (cid:48) ( t ) = − rm ( t ) + rm ( t − τ ) + rm ( t − τ ) (cid:18) K − p ( t ) − m ( t ) K (cid:19) , (5) p (cid:48) ( t ) = rm ( t ) − rm ( t − τ ) , (6)where m is the density of motile agents, p is the density of proliferative agents, K >
Ruth E. Baker, Gergely R¨ost as K = (cid:96) n ), and (cid:48) denotes differentiation with respect to t . The equations encodethe model assumption that motile agents become proliferative with rate r , and stayin this proliferative phase for time τ . Once the cell cycle is completed, proliferativecells switch back to being motile cells again. As they switch, they attempt to placea daughter agent into a lattice site; the probability that site is empty (in themean-field limit) is ( K − p ( t ) − m ( t )) /K .System (5)–(6) can be rewritten as a scalar equation using the natural biolog-ical relation p ( t ) = (cid:90) tt − τ rm ( s ) d s, (7)that is, the proliferative cells at time t are exactly those who entered the prolifera-tive subpopulation in the time interval [ t − τ, t ]. Using this relation in equation (5)gives m (cid:48) ( t ) = − rm ( t ) + rm ( t − τ ) + rm ( t − τ ) (cid:18) − rK (cid:90) tt − τ m ( s ) d s − m ( t ) K (cid:19) . (8)Next we rescale the density with K by letting ˜ m ( t ) = m ( t ) /K ,˜ m (cid:48) ( t ) = − r ˜ m ( t ) + r ˜ m ( t − τ ) + r ˜ m ( t − τ ) (cid:18) − r (cid:90) tt − τ ˜ m ( s ) d s − ˜ m ( t ) (cid:19) , (9)and drop the tildes and rearrange to m (cid:48) ( t ) = − rm ( t ) + rm ( t − τ ) (cid:18) − r (cid:90) tt − τ m ( s ) d s − m ( t ) (cid:19) . (10)To complete the non-dimensionalisation, we also rescale time by letting ˜ t = t/τ and x (˜ t ) = m ( t ). Then, d x/ d˜ t = τ m (cid:48) ( t ) and, moreover, m ( t − τ ) = m ( τ ˜ t − τ ) = m ( τ (˜ t − x (˜ t − (cid:90) tt − τ m ( s ) d s = τ (cid:90) ˜ t ˜ t − x (˜ s ) d˜ s, (11) we can, again, drop the tildes, to arrive at the new delayed logistic equation:d x d t = − rτ x ( t ) + rτ x ( t − (cid:18) − rτ (cid:90) tt − x ( s ) d s − x ( t ) (cid:19) . (12) lobal dynamics of a novel delayed logistic equation arising from cell biology 7 With the notation ρ = rτ >
0, we now analyse the scalar differential equation withdiscrete and distributed delays x (cid:48) ( t ) = − ρx ( t ) + ρx ( t − (cid:18) − ρ (cid:90) tt − x ( s ) d s − x ( t ) (cid:19) . (13)Let C = C ([ − , , R ) denote the Banach space of continuous real-valued functionson the interval [ − ,
0] equipped with the supremum norm. Then equation (13) isof the form x (cid:48) ( t ) = f ( x t ), where the solution segment x t ∈ C is defined by x t ( s ) := x ( t + s ) , s ∈ [ − , . (14)Then the map f : C → R is given by f ( φ ) = − ρφ (0) + ρφ ( − (cid:18) − ρ (cid:90) − φ ( s ) d s − φ (0) (cid:19) , (15)and the initial data is specified by x = φ ∈ C. (16)3.1 Well-posednessBy a solution of equations (13),(16) we mean a continuous function x ( t ) on aninterval [ − , A ) with 0 < A ≤ ∞ , which is differentiable on (0 , A ), satisfies equa-tion (13) on (0 , A ) and also satisfies equation (16). Lemma 3.1
For every φ ∈ C , there exists a unique solution of the initial value prob-lem (13) , (16) defined on an interval [ − , A ) for some < A ≤ ∞ , that dependscontinuously on the initial data.Proof Recall the standard existence and uniqueness theorem for functional differ- ential equations [37]. We shall verify that the local Lipschitz property holds for f , i.e. for any M >
0, there is an
L > | f ( φ ) − f ( ψ ) | ≤ L || φ − ψ || , whenever || φ || < M, || ψ || < M. (17) Ruth E. Baker, Gergely R¨ost
For a constant
M >
0, let || φ || < M and || ψ || < M , then we have the estimates | f ( φ ) − f ( ψ ) | ≤ ρ || φ − ψ || + 2 ρ || φ − ψ || + ρ (cid:12)(cid:12)(cid:12)(cid:12) φ ( − (cid:18) ρ (cid:90) − φ ( s ) d s + φ (0) (cid:19) − ψ ( − (cid:18) ρ (cid:90) − ψ ( s ) d s + ψ (0) (cid:19)(cid:12)(cid:12)(cid:12)(cid:12) ≤ ρ || φ − ψ || + ρ (cid:12)(cid:12)(cid:12)(cid:12) φ ( − (cid:18) ρ (cid:90) − φ ( s ) d s + φ (0) (cid:19) − ψ ( − (cid:18) ρ (cid:90) − φ ( s ) d s + φ (0) (cid:19)(cid:12)(cid:12)(cid:12)(cid:12) + ρ (cid:12)(cid:12)(cid:12)(cid:12) ψ ( − (cid:18) ρ (cid:90) − φ ( s ) d s + φ (0) (cid:19) − ψ ( − (cid:18) ρ (cid:90) − ψ ( s ) d s + ψ (0) (cid:19)(cid:12)(cid:12)(cid:12)(cid:12) ≤ (cid:16) ρ + ρ M + ρM + Mρ + ρM (cid:17) || φ − ψ || . (18)Hence L = 3 ρ + M (cid:0) ρ + 2 ρ (cid:1) is a valid Lipschitz constant. Existence, uniquenessand continuous dependence then follow from the general theory, c.f. Theorem 3.7of [37] and [27, Chapter 2]. (cid:117)(cid:116) Due to biological constraints, we are interested only in non-negative solutions, i.e. x ( s ) ≥ s ∈ [ − , x ∈ Ω , where Ω := (cid:26) φ ∈ C : φ ( s ) ≥ s ∈ [ − ,
0] and φ (0) + ρ (cid:90) − φ ( s ) d s ≤ (cid:27) . (19)This set in fact corresponds to the biologically feasible phase space since, for asolution with x t ∈ Ω , θ ( t ) := x ( t ) + ρ (cid:90) tt − x ( s ) d s (20)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 shalluse Ω ⊂ C , the set of biologically feasible states as our phase space, noting that Ω depends on the parameter ρ .3.2 Positivity and boundedness from above Lemma 3.2
The set Ω is positively invariant.Proof We claim that if θ (0) ≤ θ ( t ) ≤ t >
0. Note that θ (cid:48) ( t ) = − ρx ( t ) + ρx ( t − (cid:18) − ρ (cid:90) tt − x ( s ) d s − x ( t ) (cid:19) + ρx ( t ) − ρx ( t − , (21) lobal dynamics of a novel delayed logistic equation arising from cell biology 9 which simplifies to θ (cid:48) ( t ) = ρx ( t − − θ ( t )) . (22)For w ( t ) := 1 − θ ( t ) we have w (cid:48) ( t ) = − θ (cid:48) ( t ) = − ρx ( t − w ( t ), hence w ( t ) = w (0) exp (cid:18) − ρ (cid:90) t x ( s −
1) d s (cid:19) ≥ , (23)whenever w (0) ≥
0, or equivalently θ (0) ≤
1. We find that θ ( t ) = x ( t ) + ρ (cid:90) tt − x ( s ) d s ≤ , (24)holds if x (0) + ρ (cid:90) − x ( s ) d s ≤ . (25)The general positivity condition, i.e. that x t ≥ x ≥
0, for system(13),(16) is that x t ≥ x ( t ) = 0 implies f ( x t ) ≥ x ( t ) = 0 we have f ( x t ) = ρx ( t − − θ ( t )) ≥ x t ≥ θ ( t ) ≤
1. Since x t ∈ Ω is equivalent to x t ≥ θ ( t ) ≤
1, we obtainthat Ω is positively invariant, and additionally the estimate x ( t ) < (cid:117)(cid:116) Consequently, solutions starting from φ ∈ Ω exist globally, i.e. on the interval[ − , ∞ ), with x φt ∈ Ω for all t ≥ Theorem 3.1
Equation (13) has two steady states, x = 0 which is always unstable,and x ∗ = 1 / ( ρ + 1) which is always locally asymptotically stable.Proof Setting x (cid:48) = 0, we obtain the steady state equation0 = − ρx + ρx (2 − ρx − x ) , (26)which has the solutions x = 0, and x ∗ = 1 / ( ρ + 1) . Linearisation of equation (13) around x = 0 gives x (cid:48) ( t ) = − ρx ( t ) + 2 ρx ( t − , (27)with the characteristic equation λ + ρ = 2 ρe − λ , (28) which always has a positive real root, hence 0 is unstable.Next we look at the stability of the equilibrium x ∗ . Using the definition inequation (15), we have f ( x ∗ + φ ) = − ρ ( φ (0) + x ∗ ) + ρ ( φ ( −
1) + x ∗ ) (cid:18) − ρ (cid:90) − ( φ ( s ) + x ∗ ) d s − ( φ (0) + x ∗ ) (cid:19) = L ( φ ) + g ( φ ) , (29)which has linear part L ( φ ) = − ρ (cid:18) ρ + 1 (cid:19) φ (0) + ρφ ( − − ρ ρ + 1 (cid:90) − φ ( s ) d s, (30)and g ( φ ) = ρφ ( − (cid:18) − ρ (cid:90) − φ ( s ) d s − φ (0) (cid:19) . (31)Clearly lim || φ ||→ | g ( φ ) | / || φ || = 0, hence the linear variational equation is y (cid:48) ( t ) = − ρ (cid:18) ρ + 1 (cid:19) y ( t ) + ρy ( t − − ρ ρ + 1 (cid:90) tt − y ( s ) d s. (32)The exponential ansatz y ( t ) = e λt gives the characteristic equation λ = − ρ (cid:18) ρ + 1 (cid:19) − ρ ρ + 1 (cid:90) tt − e λ ( s − t ) d s + ρe − λ , (33)which becomes, after integration and rearranging, λ = − ρ (cid:18) ρ + 1 (cid:19) λ − ρ ρ + 1 (1 − e − λ ) + ρλe − λ , (34)provided λ (cid:54) = 0. Note that λ = 0 is not a root of equation (33), otherwise0 = − ρ (cid:18) ρ + 1 (cid:19) − ρ ρ + 1 + ρ = − ρ < , (35)which is a contradiction. Therefore we may write the characteristic equation as χ ( λ ) = P ( λ ) + Q ( λ ) e − λ = 0 , (36) with P ( λ ) = λ + ρ (cid:18) ρ + 1 (cid:19) λ + ρ ρ + 1 , (37) Q ( λ ) = − ρ ρ + 1 − ρλ. (38) lobal dynamics of a novel delayed logistic equation arising from cell biology 11 We can factor the characteristic equation as χ ( λ ) = (cid:18) λ + ρρ + 1 (cid:19) (cid:16) λ + ρ − ρe − λ (cid:17) = 0 , (39)which demonstrates that we always have the real root λ = − ρ/ ( ρ + 1) <
0. Theother roots are the roots of λ + ρ − ρe − λ = 0. This is a well-known type of charac-teristic 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 aroot (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 thanzero: assuming a root with (cid:60) λ ≥ λ (cid:54) = 0, we have ρ < | λ + ρ | = | ρe − λ | ≤ ρ, (40)which is a contradiction. Therefore the equilibrium x ∗ is locally asymptoticallystable. (cid:117)(cid:116) Theorem 4.1
Equation (13) is strongly uniformly persistent, that is there exists a δ > , independent of the initial data, such that for any solution x φ ( t ) with φ (0) > , lim inf t →∞ x ( t ) ≥ δ. (41) Proof
Let g ( ψ ) = ρψ ( − (cid:18) − ψ (0) − ρ (cid:90) − ψ ( θ ) d θ (cid:19) , (42)then equation (13) can be written as x (cid:48) ( t ) = − ρx ( t ) + g ( x t ) . (43)From the variation of constants formula, for any σ > ω > x ( σ ) = e − ρ ( σ − ω ) (cid:18) x ( ω ) + (cid:90) σω e ρ ( s − ω ) g ( x s ) d s (cid:19) . (44)We use the notation x ∞ = lim inf t →∞ x ( t ), and x φ to denote a specific solutionwith initial function φ . Assume the contrary of the statement of the theorem. Then there is a sequence φ n , such that lim n →∞ x φ n ∞ = 0. Let q ∈ (0 ,
1) suchthat q > / q > .
841 is just fine). Then there is a T n → ∞ , such that x φ n ( t ) ∈ [ qx φ n ∞ ,
1] for all t > T n −
1. Furthermore, there is a t n > T n + n such that x φ n ( t n ) ∈ [ qx φ n ∞ , q − x φ n ∞ ]. For the particular case x = x φ n , σ = t n and ω = T n , thevariation of constants formula gives x φ n ( t n ) = e − ρ ( t n − T n ) (cid:18) x φ n ( T n ) + (cid:90) t n T n e ρ ( s − T n ) g (cid:16) x φ n s (cid:17) d s (cid:19) . (45)Using the integral mean value theorem, there is an η n ∈ [ T n , t n ], such that e − ρt n (cid:90) t n T n e ρs g (cid:16) x φ n s (cid:17) d s = e − ρt n g (cid:16) x φ n η n (cid:17) (cid:90) t n T n e ρs d s = g (cid:16) x φ n η n (cid:17) ρ − (cid:16) − e ρ ( T n − t n ) (cid:17) . (46)Then we have the relation x φ n ( t n ) = e − ρ ( t n − T n ) x φ n ( T n ) + g (cid:16) x φ n η n (cid:17) ρ − (cid:16) − e ρ ( T n − t n ) (cid:17) . (47)Since x φ n ( t n ) → e − ρ ( t n − T n ) x φ n ( T n ) → (cid:16) − e T n − t n (cid:17) → n → ∞ ,necessarily g (cid:16) x φ n η n (cid:17) → n → ∞ . Since x φ n η n ∈ Ω, we have g (cid:16) x φ n η n (cid:17) ≥ ρx φ n η n ( − x φ n ( η n − → n → ∞ .Next we claim that (cid:18) − x φ n ( η n ) − ρ (cid:90) − x φ n ( η n + θ ) d θ (cid:19) → n → ∞ . (48)From the inequality x (cid:48) ( t ) ≥ − ρx ( t ), we find that x φ n ( η n + θ ) ≤ e ρ x φ n ( η n −
1) for θ ∈ [ − , x φ n ( η n ) + ρ (cid:90) − x φ n ( η n + θ ) d θ ≤ e ρ x φ n ( η n − ρ ) , (49)and the right-hand side is already shown to converge to zero. For sufficiently large n , we have the estimate g (cid:16) x φ n η n (cid:17) ≥ qρx φ n η n ( − ≥ qρqx φ n ∞ and (cid:16) − e ρ ( T n − t n ) (cid:17) > q , therefore from equation (47) we obtain q − x φ n ∞ ≥ x φ n ( t n ) ≥ q x φ n ∞ . (50)We find 1 ≥ q , which contradicts the choice of q . (cid:117)(cid:116) lobal dynamics of a novel delayed logistic equation arising from cell biology 13 Theorem 4.2
For all positive solutions, we have lim t →∞ (cid:18) x ( t ) + ρ (cid:90) tt − x ( s ) d s (cid:19) = 1 . (51) Proof
This follows from the standard comparison principle, since for large t wehave 1 > x ( t − > δ/ >
0, and by equation (22), we see that θ ( t ) is squeezedbetween θ ( t ) and ¯ θ ( t ), which are the solutions of θ (cid:48) ( t ) = ρ δ − θ ( t )) , (52)and ¯ θ (cid:48) ( t ) = ρ (1 − ¯ θ ( t )) , (53) with θ (0) = ¯ θ (0) = θ (0) >
0, each converging to unity. (cid:117)(cid:116)
Theorem 4.3
All positive solutions of equation (13) in Ω converge to the positiveequilibrium.Proof First we recall the following theorem from [20]:
Theorem A (Gy˝ori-Pituk)
Consider x (cid:48) ( t ) = ( a + a ( t )) x ( t ) + ( b + b ( t )) x ( t − τ ) , (54) where a , b ∈ R , τ > are constants and a, b : [0 , ∞ ) → R are continuous functions.Let the following assumptions be satisfied:(i) λ = a + b e − λτ has a unique root λ with largest real part, and this root λ isreal and simple;(ii) a ( t ) → , b ( t ) → as t → ∞ ;(iii) (cid:82) ∞ a ( t ) d t < ∞ , (cid:82) ∞ b ( t ) d t < ∞ ;(iv) (cid:82) ∞ (cid:12)(cid:12)(cid:12) τ a ( t ) − (cid:82) tt − τ a ( s ) d s (cid:12)(cid:12)(cid:12) d t < ∞ , (cid:82) ∞ (cid:12)(cid:12)(cid:12) τ b ( t ) − (cid:82) tt − τ b ( s ) d s (cid:12)(cid:12)(cid:12) d t < ∞ . Then, every solution x ( t ) satisfies x ( t ) = exp (cid:18)(cid:90) t λ ( s ) d s (cid:19) ( ξ + o (1)) , t → ∞ , (55) where λ ( t ) = λ + (cid:16) b τ e − λ τ (cid:17) − (cid:16) a ( t ) + e − λ τ b ( t ) (cid:17) , (56) and ξ = ξ ( x ) is a constant depending on the solution x . Fix some ψ ∈ Ω , and let x ψ ( t ) be the corresponding solution with w ψ ( t ) = 1 − x ψ ( t ) − ρ (cid:90) tt − τ x ψ ( s ) d s = w (0) exp (cid:18) − ρ (cid:90) t x ψ ( s −
1) d s (cid:19) , (57)where the last equality was given in equation (23). Now consider the equation z (cid:48) ( t ) = − ρz ( t ) + ρz ( t − (cid:16) w ψ ( t ) (cid:17) . (58) We apply Theorem A, and check all the conditions, where we have a = − ρ , b = ρ , τ = 1, a ( t ) = 0 and b ( t ) = ρw ψ ( t ). The conditions with a ( t ) are trivial, and theconditions with b ( t ) follow from the combination of persistence and equation (57)as follows, since w is exponentially convergent. (i) λ = − ρ + ρe − λ has a unique root λ with largest real part, and this root λ is real and simple (in fact λ = 0). (ii) To see that w ψ ( t ) → t → ∞ , we can use Theorem 4.1: there is a δ > T such that x ( t ) > δ for t > T −
1. From equation (57) we find w ψ ( t ) = w (0) exp (cid:32) − ρ (cid:90) T x ψ ( s −
1) d s (cid:33) exp (cid:18) − ρ (cid:90) tT x ψ ( s −
1) d s (cid:19) < w (0) exp (cid:32) − ρ (cid:90) T x ψ ( s −
1) d s (cid:33) e − ρδ ( t − T ) , (59)or w ψ ( t ) < Ke − ρδt , K = w (0) exp (cid:32) − ρ (cid:90) T x ψ ( s −
1) d s (cid:33) e rδT . (60) (iii) The estimate in equation (60) shows that (cid:90) ∞ (cid:16) w ψ ( t ) (cid:17) d t < K (cid:90) ∞ e − ρδt d t = K ρδ < ∞ . (61) lobal dynamics of a novel delayed logistic equation arising from cell biology 15 (iv) From the estimate in equation (60), we also find that (cid:90) ∞ (cid:12)(cid:12)(cid:12)(cid:12) w ψ ( t ) − (cid:90) tt − w ψ ( s ) d s (cid:12)(cid:12)(cid:12)(cid:12) d t < (cid:90) ∞ w ψ ( t ) d t + (cid:90) ∞ (cid:90) tt − w ψ ( s ) d s d t < ∞ . (62)As a result, λ ( t ) = (1 + ρ ) − ρw ψ ( t ) , (63)and all the conditions of Theorem 4.3 hold. Thus, every solution z ( t ) satisfies z ( t ) = exp (cid:18) (1 + ρ ) − ρ (cid:90) t w ψ ( s ) d s (cid:19) ( ξ + o (1)) , t → ∞ . (64)Now notice that a solution x ( t ) of equation (13) with initial function ψ is also asolution of equation (58), hence it converges to a constant. In view of Theorem4.1, this constant can only be the positive equilibrium.4.4 The global attractor Theorem 4.4
The global attractor A consists of the two equilibria and a heteroclinicorbit connecting the two.Proof First we demonstrate that the global attractor A of Φ Ω exists. From theinvariance of the bounded set Ω , this semiflow generated by the equation on Ω is point dissipative. It follows from the Arzel`a-Ascoli theorem that the solutionoperators Φ Ωt are completely continuous for t ≥
1, then by applying Theorem 3.4.8of [21], the compact global attractor exists.The two equilibria are part of the global attractor. Next we show that thereexists a unique heteroclinic orbit in Ω connecting x = 0 and the positive equi-librium, x ∗ . Recall that the characteristic equation of the linearization at x = 0is λ + ρ = 2 ρe − λ , which has a single positive root λ . The other roots form a se- quence of complex conjugate pairs ( λ j , ¯ λ j ) with (cid:60) λ j +1 < (cid:60) λ j < λ for all integers j ≥
1. We refer to Chapter XI of [6] for a complete analysis of the characteristicequation of the form z − α − βe − z , in particular Fig. XI.1. in [6] which summarizesits properties. Our equation is a special case with α = − ρ , β = 2 ρ . There is a γ > λ > γ > (cid:60) λ j for all j ≥
1, and there is a k ∈ N suchthat (cid:60) λ k > (cid:60) λ k +1 ≤
0, where the point ( − ρ, ρ ) lies between the curves C − k and C − k +1 on the ( α, β )-plane, given by C − j = (cid:110) ( α, β ) = (cid:16) ν cos ν sin ν , − ν sin ν (cid:17) : ν ∈ ((2 j − π, jπ ) (cid:111) . (65)The intersections of C − j and the half-line ( − ρ, ρ ) ( ρ >
0) are determined bycos ν = 1 /
2, thus ν = 2 jπ − π/
3, sin ν = −√ / ρ = ρ j = (2 jπ − π/ / √ . Hence, k is either zero or the largest positive integersuch that ρ > ρ k .Considering the leading real root, the corresponding eigenfunction is given by χ ( s ) := e λ s , s ∈ [ − , . The phase space C can be decomposed as C = P ⊕ P k ⊕ Q ,where the function χ ∈ C spans the linear eigenspace P := { cχ : c ∈ R } . P k is a 2 k -dimensional eigenspace corresponding to { λ j : j = 1 , . . . , k } . If λ is theonly eigenvalue with positive real part then P k = ∅ . There is an ε > (cid:60) λ k > ε . Q corresponds to the remaining part of the spectrum.There exist open neighborhoods N , N k , M in P , P k , Q , respectively, and C -maps w : N → P k ⊕ Q , w u : N ⊕ N k → Q with range in N k ⊕ M and M ,respectively, such that w (0) = 0, Dw (0) = 0, w u (0) = 0, Dw u (0) = 0; where D denotes the Fr´echet-derivative. Then, the γ -unstable set of the equilibrium x = 0,namely W (0) := (cid:110) φ ∈ N + N k + M : there is a trajectory z t , t ∈ R with z = φ,z t ∈ N + N when t ≤ z ( t ) e − γt → t → −∞ (cid:111) , (66)coincides with the graph W := { φ + w ( φ ) : φ ∈ N } . (67)Furthermore, the unstable set of the equilibrium x = 0 W u (0) := (cid:110) φ ∈ N + N k + M : there is a trajectory z t , t ∈ R with z = φ,z t ∈ N + N when t ≤ z ( t ) e − εt → t → −∞ (cid:111) , (68)coincides with the graph W u := { φ + w u ( φ ) : φ ∈ N + N k } . (69) lobal dynamics of a novel delayed logistic equation arising from cell biology 17 For the details see [10,26].For any φ ∈ N , there is a c ∈ R such that φ = cχ . We have || χ || = 1 and χ ( t ) ≥ e − λ > t ∈ [ − , Dw (0) = 0 thatlim || φ ||→ || w ( φ ) |||| φ || = 0 , φ ∈ N , (70)which means lim c → || w ( cχ ) ||| c | = 0 . (71)Therefore, there exists a c > || w ( cχ ) || / | c | < e − λ / c ∈ (0 , c ), or equivalently || w ( cχ ) || < ce − λ /
2. We may assume that c < / (2(1 + ρ )).Thenmin { cχ ( t ) + w ( cχ )( t ) : t ∈ [ − , , c ∈ (0 , c ) } ≥ ce − λ − c e − λ = c e − λ > , (72)and cχ (0) + w ( cχ )(0) + r (cid:90) − { cχ ( s ) + w ( cχ )( s ) } d s ≤ (1 + ρ )( c + || w ( cχ ) || ) < c (1 + ρ ) < , (73)thus φ c := cχ + w ( cχ ) ∈ W ∩ Ω for all c ∈ (0 , c ) . (74)The unstable set W (0) intersects Ω , and for any function φ c of this intersec-tion, there is a complete solution x ( t ) : R → R +0 such that x = φ c and x t → ∗ as t → −∞ . By the negative invariance of W (0), the trajectory is completely in Ω , hence we have proved the existence of a heteroclinic orbit connecting the twoequilibria.For uniqueness, consider the unstable manifold W u (0). Any solution on the set W u (0) \W (0) oscillates about 0 as t → −∞ , hence can not be a positive heteroclinicsolution. Denote the unique heteroclinic orbit in Ω by H . We claim that H and the two equilibria constitute the global attractor, otherwise there are further completeorbits in Ω . Assume there is a ψ ∈ Ω which is on such a complete orbit Γ ⊂ Ω ,so ψ (cid:54) = 0, ψ (cid:54) = x ∗ , and ψ / ∈ H . Then its α -limit set α ( ψ ) and ω -limit set ω ( ψ ) arenon-empty, compact and invariant. In particular, by Theorem 4.3, ω ( ψ ) = { x ∗ } .
190 192 194 196 198 200t0.0150.0200.025x 190 192 194 196 198 200t0.0050.0100.015x
190 192 194 196 198 200t0.0020.0040.0060.0080.0100.012x
Fig. 1
Representative solutions that appear periodic over short time scales. From left to right, ρ = 50 , ,
190 with φ ( t ) = 0 . t ) + 1). If α ( ψ ) = { x } , then Γ is a connecting orbit coinciding with H . Hence, we mayassume that there is a φ (cid:54) = 0 such that φ ∈ α ( ψ ). Then x φt → x ∗ as t → ∞ ,so x ∗ ∈ α ( ψ ) as well. But this contradicts the stability of x ∗ . Thus, the globalattractor does not contain any other orbit, and A = { x } ∪ { x ∗ } ∪ H . Representative numerical simulations of solutions of equation (13) are plottedin Fig. 1, with initial function φ ( t ) = 0 . t ) + 1) and parameter values ρ = 50 , , in [33] for a scalar differential equation with a constant delay. It is interesting tosee that equation (13) can exhibit rapid convergence for small ρ and oscillatorypatterns of various shapes for larger ρ , as illustrated in Fig. 3 where solutions areplotted with different initial functions. lobal dynamics of a novel delayed logistic equation arising from cell biology 19
20 22 24 26 28 30t0.0050.0100.0150.020x 220 222 224 226 228 230t0.0050.0100.0150.020x
820 822 824 826 828 830t0.0050.0100.0150.020x
Fig. 2
Top: envelope of the solution with initial function φ ( t ) = 0 . t ) + 1) and ρ = 100, showing slow convergence to the positive equilbrium. Bottom: snapshots of the samesolution over different time intervals.
190 191 192 193 194 195 t0.0040.0060.0080.0100.0120.0140.016x
Fig. 3
Left: Convergence for ρ = 10 with initial function φ ( t ) = 0 . t ) + 1). Right:Various periodic patterns for ρ = 100, emerging from initial functions φ ( t ) = 0 . t )+1)(dotted), φ ( t ) = 0 . t ) + 1) (solid), φ ( t ) = 0 . t ) + 1) (dashed). To gain some understanding of this phenomenon, we look at our model fromthe point of view of singular perturbations. In equation (13), let us divide by ρ .For large ρ , we use the notation ε = 1 /ρ , then we have ε x (cid:48) ( t ) = ε [ − x ( t ) + 2 x ( t − − x ( t − x ( t )] − x ( t − (cid:90) tt − x ( s ) d s. (75)Setting ε = 0, we obtain the singular perturbation0 = x ( t − (cid:90) tt − x ( s ) d s. (76) - - - - - - - I m - - - - - - - I m - - - - - - - I m Fig. 4
Characteristic roots of (33) on the complex plane for ρ = 1 , , ρ → ∞ ,characteristic roots converge towards the imaginary axis. Besides x ≡
0, any periodic function with unit period that is on average zero isa solution of this equation, which may be a hint as to why we see long-lastingtransient oscillations with a variety of shapes.In addition, we can look at the singular perturbation of the characteristicequation: equation (39) can be written as λ + ρ (cid:18) ρ + 1 (cid:19) λ + ρ ρ + 1 = (cid:18) ρ ρ + 1 + ρλ (cid:19) e − λ . (77)Multiplying by ( ρ + 1) /ρ we have λ ρ + 1 ρ + ρ + 2 ρ λ + 1 = (cid:18) λ ρ + 1 ρ (cid:19) e − λ , (78)which is λ ( ε + ε ) + (1 + 2 ε ) λ + 1 = (1 + λ (1 + 2 ε )) e − λ , (79)with the notation ε = 1 /ρ . Setting ε = 0, we obtain the singular perturbation λ + 1 = (1 + λ ) e − λ . (80)The roots of this equation satisfy either λ = − e − λ = 1, having the purelyimaginary roots i = 2 kπ , k ∈ Z . Indeed, as we increase ρ , we can see that allcomplex roots line up along the imaginary axis, see Fig. 4. The delayed logistic equation has received much attention in past decades in theanalysis of nonlinear delay differential equations. However, its biological validity lobal dynamics of a novel delayed logistic equation arising from cell biology 21
20 25 30 35 40 t0.0476120.0476140.0476160.0476180.047620x
10 20 30 40 t0.010.020.030.040.050.06x
Fig. 5
Illustration of the unique positive heteroclinic orbit for ρ = 20. Left: the initial phaseof the solution (solid) is plotted alongside the exponential eigenfunction (dashed) correspond-ing to the leading real eigenvalue at zero. Center: we zoomed in to see the late part of thesolution (solid), aligning nicely with the oscillatory pattern corresponding to the leading pairof complex eigenvalues at the positive equilbrium (dashed). Right: the heteroclinic solution isshown combining the two time scales. The dotted line is the positive equilbrium in all threefigures. has been questioned, despite the fact that it was introduced by Hutchinson toexplain observations of oscillatory behaviour in ecological systems. Here, we intro-duced a new logistic-type delay differential equation, derived from the go-or-growhypothesis which has been observed for some types of cancer cells. The equationnaturally includes a product of delayed and non-delayed terms, as well as bothdiscrete and distributed delays. Detailed mathematical analysis reveals that allnon-zero solutions converge to the stable positive equilibrium. For an alternativeformulation of the delayed logistic equation [1], the authors also concluded that,as opposed to Hutchinson’s equation, sustained oscillations were not possible andsolutions settled at an equilibrium. In some sense the dynamical behaviour of ourmodel is in-between these two extremes. While we have proved that all positivesolutions are attracted to the positive equilibrium, for a range of parameters theconvergence is so slow that, on a biologically realistic time scale, solutions mayappear periodic. These long-lasting transient patterns can have various shapes, asshown in Fig. 1, and these shapes also change in time (see Fig. 2).We have also fully described the global attractor as the union of the two equi-libria and a unique connecting orbit. This heteroclinic solution is depicted in Fig. ρ = 20. In this case, the leading real eigenvalue at zero is λ ≈ .
66, whilethe leading pair of eigenvalues at the positive equilibrium are λ ≈ − . ± i . Wecan numerically observe how the solution is aligned to the leading eigenspaces ofthe linearizations near the two equilibria. Our model is based on the commonly-invoked modelling assumption that pro-liferating cells abort proliferation when placement of a daughter cell is not possibledue to spatial crowding constraints. However, there are other potential models ofcrowding-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 cellbecomes available, or that, in order to enter the proliferating state, cells must beable to “reserve” a site for the daughter cell to be placed into. These, and otherpossible 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 thego-or-grow mechanism when a delay caused by the cell cycle length is explicitlyincorporated as a biological parameter. Future work will explore more accuratemodels of the effects of spatial crowding upon proliferation, using both mean-fieldand moment dynamics models. In the transformed equation (13), the parameter ρ is the product of the time delay, τ , and the proliferation rate, r . As such, increasingeither of those two parameters will have a similar effect on the dynamics (althoughon different time scales in terms of the original, dimensional equation). In modelsthat include spatial details, however, there is an additional key parameter (thatis ignored by equation (13)), namely the cell motility rate. An important futuregoal is to understand how the speed of e.g. cancer invasion depends on cell-levelparameters 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 cellcycle delays need to be developed. Travelling wave solutions of reaction-diffusion systems with delays are closely related to the heteroclinic orbits of the reactionsystems (see for example [11]). As such, the results we provide on the uniqueheteroclinic orbit connecting the zero and the positive equilibria will help shedlight on invasion speeds in these cases. lobal dynamics of a novel delayed logistic equation arising from cell biology 23 Acknowledgments
GR was supported by NKFI FK 124016 and MSCA-IF 748193. REB is a RoyalSociety Wolfson Research Merit Award holder and would like to thank the Lever-hulme Trust for a Research Fellowship.
References
1. Arino J., Wang L. and Wolkowicz G. S. (2006). An alternative formulation for a delayedlogistic equation. Journal of Theoretical Biology, 241(1):109–119.2. Baca¨er N. (2011). A short history of mathematical population dynamics. Springer Science& Business Media.3. Baker R. E. and Simpson M. J. (2010). Correcting mean-field approximations for birth-death-movement processes. Physical Review E 82(4):e041905.4. B´anhelyi B., Csendes T., Krisztin T. and Neumaier A. (2014). Global attractivity ofthe zero solution for Wright’s equation. SIAM Journal on Applied Dynamical Systems,13(1):537–563.5. Chow S.-N. and Mallet-Paret J. (1977). Integral averaging and bifurcation, Journal ofDifferential Equations 26:112–0159.6. Diekmann O., Van Gils S. A., Lunel S. M. and Walther, H. O. (2012). Delay equations:functional, complex, and nonlinear analysis (Vol. 110). Springer Science & Business Media.7. Erneux T. (2009). Applied delay differential equations (Vol. 3). Springer Science & Busi-ness Media.8. Faria T. (2006). Asymptotic stability for delayed logistic type equations. Mathematicaland computer modelling, 43(3-4), 433–445.9. Faria T., Huang W., and Wu J. (2006). Travelling waves for delayed reaction-diffusionequations with global response. Proceedings of the Royal Society of London A: Mathe-matical, Physical and Engineering Sciences 462(2065):229–261.10. Faria T., Huang W., and Wu J. (2002). Smoothness of center manifolds for maps and formaladjoints for semilinear FDEs in general Banach spaces. SIAM Journal on MathematicalAnalysis, 34(1):1730-203.11. Faria T. and Trofimchuk S. (2006). Nonmonotone travelling waves in a single speciesreaction-diffusion equation with delay. Journal of Differential Equations, 228(1):357–376.12. Farin A., Suzuki S. O., Weiker M., Goldman J. E., Bruce J. N., and Canoll P. (2006).Transplanted glioma cells migrate and proliferate on host brain vasculature: a dynamicanalysis. Glia, 53(8):799–808.13. Fowler A. C. (1982). An asymptotic analysis of the delayed logistic equation when thedelay is large. IMA Journal of Applied Mathematics, 28(1):41–49.14. Geritz S.A. and Kisdi ´E. (2012). Mathematical ecology: why mechanistic models? Journalof Mathematical Biology, 65(6):1411–1415.15. Giese A., Bjerkvig R., Berens M. E. and Westphal M. (2003). Cost of migration: invasion ofmalignant gliomas and implications for treatment. Journal of Clinical Oncology 21:1624–1636.16. Gopalsamy K., and Zhang B. G. (1988). On a neutral delayed logistic equation. Dynamicsand Stability of Systems, 2(3-4):183–195.17. Grotta-Ragazzo C., Malta C. P., and Pakdaman K. (2010). Metastable periodic patterns insingularly perturbed delayed equations. Journal of Dynamics and Differential Equations,22(2):203–252.18. Gurney W. S. C., Blythe S. P. and Nisbet, R. M. (1980). Nicholson’s blowflies revisited.Nature, 287(5777):17–21.19. Gy˝ori I., Nakata Y., R¨ost G. (2018). Unbounded and blow-up solutions for a delayedlogistic equation with positive feedback. Communications on Pure and Applied Analysis17(6):2845–2854.20. Gy˝ori I., Pituk M. (1995). L2