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

TL;DR

This paper introduces a biologically motivated delayed logistic equation with discrete and distributed delays, thoroughly analyzing its global dynamics and showing solutions tend to a positive equilibrium while exhibiting complex transient oscillations.

Contribution

The study presents a new delayed logistic model rooted in cell biology, providing a complete global dynamics analysis and novel insights into its long-term behavior.

Findings

All feasible solutions converge to the positive equilibrium.

A heteroclinic orbit connects zero and the positive equilibrium.

The system exhibits long-lasting transient oscillations.

Abstract

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.

Figures

Equations

$$v^{\prime }(t)=\alpha v(t)\left[1-v(t-\tau )\right],$$

$$u^{\prime }(t)=-\alpha u(t-\tau )\left[1+u(t)\right].$$

$$N^{\prime }(t)=aN(t-\tau )e^{-bN(t-\tau )}-\mu N(t),$$

$$N^{\prime }(t)=\frac{\gamma \mu N(t-\tau )}{\mu e^{\mu \tau }+\kappa (e^{\mu \tau }-1)N(t-\tau )}-\mu N(t)-\kappa N(t{)}^2,$$

$$m^{\prime }(t)$$

$$p^{\prime }(t)$$

$$p(t)={\int }_{t-\tau }^{t}rm(s)\text{d}s,$$

$$m^{\prime }(t)=-rm(t)+rm(t-\tau )+rm(t-\tau )\left(1-\frac{r}{K}{\int }_{t-\tau }^{t}m(s)\text{d}s-\frac{m(t)}{K}\right).$$

$${\tilde{m}}^{\prime }(t)=-r\tilde{m}(t)+r\tilde{m}(t-\tau )+r\tilde{m}(t-\tau )\left(1-r{\int }_{t-\tau }^t\tilde{m}(s)\text{d}s-\tilde{m}(t)\right),$$

$$m^{\prime }(t)=-rm(t)+rm(t-\tau )\left(2-r{\int }_{t-\tau }^{t}m(s)\text{d}s-m(t)\right).$$

$${\int }_{t-\tau }^{t}m(s)\text{d}s=\tau {\int }_{\tilde{t}-1}^{\tilde{t}}x(\tilde{s})\text{d}\tilde{s},$$

$$\frac{\text{d}x}{\text{d}t}=-r\tau x(t)+r\tau x(t-1)\left(2-r\tau {\int }_{t-1}^{t}x(s)\text{d}s-x(t)\right).$$

$$x^{\prime }(t)=-\rho x(t)+\rho x(t-1)\left(2-\rho {\int }_{t-1}^{t}x(s)\text{d}s-x(t)\right).$$

$$x_t(s):=x(t+s),s\in [-1,0].$$

$$f(\varphi )=-\rho \varphi (0)+\rho \varphi (-1)\left(2-\rho {\int }_{-1}^0\varphi (s)\text{d}s-\varphi (0)\right),$$

$$x_0=\varphi \in C.$$

$$|f(\varphi )-f(\psi )|\le L||\varphi -\psi ||,\text{ whenever }||\varphi ||<M,||\psi ||<M.$$

$$|f(\varphi )-f(\psi )|$$

$$\Omega :=\left\{\varphi \in C:\varphi (s)\ge 0\text{ for }s\in [-1,0]\text{ and }\varphi (0)+\rho {\int }_{-1}^0\varphi (s)\text{d}s\le 1\right\}.$$

$$\theta (t):=x(t)+\rho {\int }_{t-1}^{t}x(s)\text{d}s$$

$${\theta }^{\prime }(t)=-\rho x(t)+\rho x(t-1)\left(2-\rho {\int }_{t-1}^{t}x(s)\text{d}s-x(t)\right)+\rho x(t)-\rho x(t-1),$$

$${\theta }^{\prime }(t)=\rho x(t-1)(1-\theta (t)).$$

$$w(t)=w(0)\exp \left(-\rho {\int }_0^{t}x(s-1)\text{d}s\right)\ge 0,$$

$$\theta (t)=x(t)+\rho {\int }_{t-1}^{t}x(s)\text{d}s\le 1,$$

$$x(0)+\rho {\int }_{-1}^{0}x(s)\text{d}s\le 1.$$

$$0=-\rho x+\rho x\left(2-\rho x-x\right),$$

$$x^{\prime }(t)=-\rho x(t)+2\rho x(t-1),$$

$$\lambda +\rho =2\rho e^{-\lambda },$$

$$f(x_{\ast }+\varphi )$$

$$L(\varphi )=-\rho \left(1+\frac{1}{\rho +1}\right)\varphi (0)+\rho \varphi (-1)-\frac{{\rho }^2}{\rho +1}{\int }_{-1}^0\varphi (s)\text{d}s,$$

$$g(\varphi )=\rho \varphi (-1)\left(-\rho {\int }_{-1}^0\varphi (s)\text{d}s-\varphi (0)\right).$$

Full text

Keywords: cell population, logistic growth, go or grow models, time delay,

global attractor, stability, long transient.

Mathematical Subject Classfication: 92D25, 34K20.

1 Introduction

The well known logistic differential equation $N^{\prime}(t)=rN(t)\left(1-N(t)/K\right)$ was proposed by Verhulst in 1838 to resolve the Malthusian dilemma of unbounded growth bacaer . 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 hutchinson to propose the delayed logistic equation (also known as Hutchinson’s equation)

[TABLE]

where we have normalized the carrying capacity to unity (by setting $v(t)=N(t)/K$), and introduced the time delay $\tau>0$. Interestingly, independent of Hutchinson, in 1955 Wright wright , 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

[TABLE]

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 $\alpha\tau>\pi/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 may . However, the use of Hutchinson’s delayed logistic equation received criticisms from biological modellers due to the lack of a mechanistic biological derivation GeritzKisdi , 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. gurney proposed the equation

[TABLE]

which is known today as Nicholson’s blowfly equation. It is structurally different from the delayed logistic equation, and the terms have clear biological interpretation as the birth and mortality rates, while the delay represents the period between egg laying and hatching. Furthermore, in contrast to the delayed logistic equation, Nicholson’s blowfly equation can explain not only the appearance of cycling populations, but also the two bursts of reproductive activity per adult population cycle that appeared in Nicholson’s data.

Incorporating a delay in the growth term, Arino et al. arino derived an alternative formulation of the delayed logistic equation

[TABLE]

which has recently been extended to distributed delays lin . 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 lindstrom , and obtained equations of logistic type with delays as a limiting case. Lindström’s equation shares the same dynamical properties as the alternative delayed logistic equation from arino , 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 chow ; faria ; hassard ; nussbaum , asymptotic analysis fowler , 3/2-type stability criteria ivanov ; wright , and the study of slowly oscillatory solutions lessard . The most famous related problem is Wright’s conjecture from 1955, which proposes global asymptotic stability in the delayed logistic equation for $\alpha\tau\leq\pi/2$. This long-standing question has been resolved recently by combining analytical and computational tools krisztin ; vBJ . Moreover, many variants and generalizations of the delayed logistic equation have been studied in the mathematical literature, for instance equation including positive instantaneous feedback gynr , non-autonomous terms faria , neutral terms gopalsamy , spatial diffusion zou , and multiple delays yan . Further examples are discussed in kuang ; liz ; ruan .

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.

2 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 $\Delta$ and length (in each direction) $\ell\Delta$, where $\ell$ 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 levchenko . Our model is designed to capture the “go or grow” hypothesis frequently acknowledged in the cancer literature farin ; giese , 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\delta{t}$ is the probability that a motile agent switches to become proliferative in the next infinitesimally small time interval $\delta{t}$. Upon switching, agents remain immobile for a fixed time $\tau>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 $P_{s}$.

Following arguments presented in Baker and Simpson Baker:2010:CMF , 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:

[TABLE]

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=\ell^{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 $\tau$. 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 mean-field limit) is $\left(K-p(t)-m(t)\right)/K$.

System (5)–(6) can be rewritten as a scalar equation using the natural biological relation

[TABLE]

that is, the proliferative cells at time $t$ are exactly those who entered the proliferative subpopulation in the time interval $[t-\tau,t]$. Using this relation in equation (5) gives

[TABLE]

Next we rescale the density with $K$ by letting $\tilde{m}(t)=m(t)/K$,

[TABLE]

and drop the tildes and rearrange to

[TABLE]

To complete the non-dimensionalisation, we also rescale time by letting $\tilde{t}=t/\tau$ and $x(\tilde{t})=m(t)$. Then, ${\text{d}x}/{\text{d}\tilde{t}}=\tau m^{\prime}(t)$ and, moreover, $m(t-\tau)=m(\tau\tilde{t}-\tau)=m(\tau(\tilde{t}-1))=x(\tilde{t}-1)$. Noting that

[TABLE]

we can, again, drop the tildes, to arrive at the new delayed logistic equation:

[TABLE]

3 Basic mathematical properties

With the notation $\rho=r\tau>0$, we now analyse the scalar differential equation with discrete and distributed delays

[TABLE]

Let $C=C([-1,0],R)$ denote the Banach space of continuous real-valued functions on the interval $[-1,0]$ equipped with the supremum norm. Then equation (13) is of the form $x^{\prime}(t)=f(x_{t})$, where the solution segment $x_{t}\in C$ is defined by

[TABLE]

Then the map $f:C\to R$ is given by

[TABLE]

and the initial data is specified by

[TABLE]

3.1 Well-posedness

By a solution of equations (13),(16) we mean a continuous function $x(t)$ on an interval $[-1,A)$ with $0<A\leq\infty$, which is differentiable on $(0,A)$, satisfies equation (13) on $(0,A)$ and also satisfies equation (16).

Lemma 3.1

For every $\phi\in C$, there exists a unique solution of the initial value problem (13),(16) defined on an interval $[-1,A)$ for some $0<A\leq\infty$, that depends continuously on the initial data.

Proof

Recall the standard existence and uniqueness theorem for functional differential equations smith . We shall verify that the local Lipschitz property holds for $f$, i.e. for any $M>0$, there is an $L>0$ such that

[TABLE]

For a constant $M>0$, let $||\phi||<M$ and $||\psi||<M$, then we have the estimates

[TABLE]

Hence $L=3\rho+M\left(2\rho^{2}+2\rho\right)$ is a valid Lipschitz constant. Existence, uniqueness and continuous dependence then follow from the general theory, c.f. Theorem 3.7 of smith and (kuang, , Chapter 2). ∎

Due to biological constraints, we are interested only in non-negative solutions, i.e. $x_{0}(s)\geq 0$, $s\in[-1,0]$. More specifically, we consider solutions with $x_{0}\in\Omega$, where

[TABLE]

This set in fact corresponds to the biologically feasible phase space since, for a solution with $x_{t}\in\Omega$,

[TABLE]

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 $\Omega\subset C$, the set of biologically feasible states as our phase space, noting that $\Omega$ depends on the parameter $\rho$.

3.2 Positivity and boundedness from above

Lemma 3.2

The set $\Omega$ is positively invariant.

Proof

We claim that if $\theta(0)\leq 1$ then $\theta(t)\leq 1$ for $t>0$. Note that

[TABLE]

which simplifies to

[TABLE]

For $w(t):=1-\theta(t)$ we have $w^{\prime}(t)=-\theta^{\prime}(t)=-\rho x(t-1)w(t)$, hence

[TABLE]

whenever $w(0)\geq 0$, or equivalently $\theta(0)\leq 1$. We find that

[TABLE]

holds if

[TABLE]

The general positivity condition, i.e. that $x_{t}\geq 0$ whenever $x_{0}\geq 0$, for system (13),(16) is that $x_{t}\geq 0$ and $x(t)=0$ implies $f(x_{t})\geq 0$ (smith, , Theorem 3.4). This clearly holds as for $x(t)=0$ we have $f(x_{t})=\rho x(t-1)(2-\theta(t))\geq 0$ following from $x_{t}\geq 0$ and $\theta(t)\leq 1$. Since $x_{t}\in\Omega$ is equivalent to $x_{t}\geq 0$ and $\theta(t)\leq 1$, we obtain that $\Omega$ is positively invariant, and additionally the estimate $x(t)<1$ holds. ∎

Consequently, solutions starting from $\phi\in\Omega$ exist globally, i.e. on the interval $[-1,\infty)$, with $x^{\phi}_{t}\in\Omega$ for all $t\geq 0$.

3.3 Steady states and their stability

Theorem 3.1

Equation (13) has two steady states, $x_{0}=0$ which is always unstable, and $x_{*}=1/(\rho+1)$ which is always locally asymptotically stable.

Proof

Setting $x^{\prime}=0$, we obtain the steady state equation

[TABLE]

which has the solutions $x_{0}=0$, and $x_{*}=1/(\rho+1).$

Linearisation of equation (13) around $x=0$ gives

[TABLE]

with the characteristic equation

[TABLE]

which always has a positive real root, hence [math] is unstable.

Next we look at the stability of the equilibrium $x_{*}$. Using the definition in equation (15), we have

[TABLE]

which has linear part

[TABLE]

and

[TABLE]

Clearly $\lim_{||\phi||\to 0}{|g(\phi)|}/{||\phi||}=0$, hence the linear variational equation is

[TABLE]

The exponential ansatz $y(t)=e^{\lambda t}$ gives the characteristic equation

[TABLE]

which becomes, after integration and rearranging,

[TABLE]

provided $\lambda\neq 0$. Note that $\lambda=0$ is not a root of equation (33), otherwise

[TABLE]

which is a contradiction. Therefore we may write the characteristic equation as

[TABLE]

with

[TABLE]

We can factor the characteristic equation as

[TABLE]

which demonstrates that we always have the real root $\lambda=-{\rho}/({\rho+1})<0$. The other roots are the roots of $\lambda+\rho-\rho e^{-\lambda}=0$. This is a well-known type of characteristic equation of the form $\lambda=A+Be^{-\tau\lambda}$, see for example Chapter 4.5 of smith . However, we are on the stability boundary of this equation given that $\lambda=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 $\Re\lambda\geq 0$ but $\lambda\neq 0$, we have

[TABLE]

which is a contradiction. Therefore the equilibrium $x_{*}$ is locally asymptotically stable. ∎

4 Global dynamics

4.1 Persistence

Theorem 4.1

Equation (13) is strongly uniformly persistent, that is there exists a $\delta>0$, independent of the initial data, such that for any solution $x^{\phi}(t)$ with $\phi(0)>0$,

[TABLE]

Proof

Let

[TABLE]

then equation (13) can be written as

[TABLE]

From the variation of constants formula, for any $\sigma>\omega>0$ we find

[TABLE]

We use the notation $x_{\infty}=\liminf_{t\to\infty}x(t)$, and $x^{\phi}$ to denote a specific solution with initial function $\phi$. Assume the contrary of the statement of the theorem. Then there is a sequence $\phi_{n}$, such that $\lim_{n\to\infty}x_{\infty}^{\phi_{n}}=0$. Let $q\in(0,1)$ such that $q^{4}>1/2$ (any $q>0.841$ is just fine). Then there is a $T_{n}\to\infty$, such that $x^{\phi_{n}}(t)\in[qx_{\infty}^{\phi_{n}},1]$ for all $t>T_{n}-1$. Furthermore, there is a $t_{n}>T_{n}+n$ such that $x^{\phi_{n}}(t_{n})\in[qx_{\infty}^{\phi_{n}},q^{-1}x_{\infty}^{\phi_{n}}]$. For the particular case $x=x^{\phi_{n}},\sigma=t_{n}$ and $\omega=T_{n}$, the variation of constants formula gives

[TABLE]

Using the integral mean value theorem, there is an $\eta_{n}\in[T_{n},t_{n}]$, such that

[TABLE]

Then we have the relation

[TABLE]

Since $x^{\phi_{n}}(t_{n})\to 0$, $e^{-\rho(t_{n}-T_{n})}x^{\phi_{n}}(T_{n})\to 0$ and $\left(1-e^{T_{n}-t_{n}}\right)\to 1$ as $n\to\infty$, necessarily $g\left(x^{\phi_{n}}_{\eta_{n}}\right)\to 0$ as well as $n\to\infty$. Since $x^{\phi_{n}}_{\eta_{n}}\in\Omega,$ we have $g\left(x^{\phi_{n}}_{\eta_{n}}\right)\geq\rho x^{\phi_{n}}_{\eta_{n}}(-1)$, hence $x^{\phi_{n}}(\eta_{n}-1)\to 0$ as well as $n\to\infty$.

Next we claim that

[TABLE]

From the inequality $x^{\prime}(t)\geq-\rho x(t)$, we find that $x^{\phi_{n}}(\eta_{n}+\theta)\leq e^{\rho}x^{\phi_{n}}(\eta_{n}-1)$ for $\theta\in[-1,0]$. Hence,

[TABLE]

and the right-hand side is already shown to converge to zero. For sufficiently large $n$, we have the estimate $g\left(x^{\phi_{n}}_{\eta_{n}}\right)\geq 2q\rho x^{\phi_{n}}_{\eta_{n}}(-1)\geq 2q\rho qx^{\phi_{n}}_{\infty}$ and $\left(1-e^{\rho(T_{n}-t_{n})}\right)>q$, therefore from equation (47) we obtain

[TABLE]

We find $1\geq 2q^{4}$, which contradicts the choice of $q$. ∎

4.2 A crucial convergence property

Theorem 4.2

For all positive solutions, we have

[TABLE]

Proof

This follows from the standard comparison principle, since for large $t$ we have $1>x(t-1)>{\delta}/{2}>0$, and by equation (22), we see that $\theta(t)$ is squeezed between $\underline{\theta}(t)$ and $\bar{\theta}(t)$, which are the solutions of

[TABLE]

and

[TABLE]

with $\underline{\theta}(0)=\bar{\theta}(0)=\theta(0)>0$, each converging to unity. ∎

4.3 Global asymptotic stability

Theorem 4.3

All positive solutions of equation (13) in $\Omega$ converge to the positive equilibrium.

Proof

First we recall the following theorem from gyoripituk :

Theorem A** (Győri-Pituk)**

Consider

[TABLE]

where $a_{0},b_{0}\in R$, $\tau>0$ are constants and $a,b:[0,\infty)\to R$ are continuous functions. Let the following assumptions be satisfied:

1. (i)

$\lambda=a_{0}+b_{0}e^{-\lambda\tau}$* has a unique root $\lambda_{0}$ with largest real part, and this root $\lambda_{0}$ is real and simple;* 2. (ii)

$a(t)\to 0$, $b(t)\to 0$ as $t\to\infty$; 3. (iii)

$\int_{0}^{\infty}a^{2}(t)\,\mathrm{d}t<\infty$, $\int_{0}^{\infty}b^{2}(t)\,\mathrm{d}t<\infty$; 4. (iv)

$\int_{0}^{\infty}\left|\tau a(t)-\int_{t-\tau}^{t}a(s)\,\mathrm{d}s\right|\mathrm{d}t<\infty$, $\int_{0}^{\infty}\left|\tau b(t)-\int_{t-\tau}^{t}b(s)\,\mathrm{d}s\right|\mathrm{d}t<\infty$.

Then, every solution $x(t)$ satisfies

[TABLE]

where

[TABLE]

and $\xi=\xi(x)$ is a constant depending on the solution $x$.

Fix some $\psi\in\Omega$, and let $x^{\psi}(t)$ be the corresponding solution with

[TABLE]

where the last equality was given in equation (23). Now consider the equation

[TABLE]

We apply Theorem A, and check all the conditions, where we have $a_{0}=-\rho$, $b_{0}=\rho$, $\tau=1$, $a(t)=0$ and $b(t)=\rho w^{\psi}(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.

1. (i)

$\lambda=-\rho+\rho e^{-\lambda}$ has a unique root $\lambda_{0}$ with largest real part, and this root $\lambda_{0}$ is real and simple (in fact $\lambda_{0}=0$). 2. (ii)

To see that $w^{\psi}(t)\to 0$ as $t\to\infty$, we can use Theorem 4.1: there is a $\delta>0$ and a $T$ such that $x(t)>\delta$ for $t>T-1$. From equation (57) we find

[TABLE]

or

[TABLE] 3. (iii)

The estimate in equation (60) shows that

[TABLE] 4. (iv)

From the estimate in equation (60), we also find that

[TABLE]

As a result,

[TABLE]

and all the conditions of Theorem 4.3 hold. Thus, every solution $z(t)$ satisfies

[TABLE]

Now notice that a solution $x(t)$ of equation (13) with initial function $\psi$ is also a solution of equation (58), hence it converges to a constant. In view of Theorem 4.1, this constant can only be the positive equilibrium.

4.4 The global attractor

Theorem 4.4

The global attractor $\mathcal{A}$ consists of the two equilibria and a heteroclinic orbit connecting the two.

Proof

First we demonstrate that the global attractor $\mathcal{A}$ of $\Phi^{\Omega}$ exists. From the invariance of the bounded set $\Omega$, this semiflow generated by the equation on $\Omega$ is point dissipative. It follows from the Arzelà-Ascoli theorem that the solution operators $\Phi^{\Omega}_{t}$ are completely continuous for $t\geq 1$, then by applying Theorem 3.4.8 of hale88 , the compact global attractor exists.

The two equilibria are part of the global attractor. Next we show that there exists a unique heteroclinic orbit in $\Omega$ connecting $x_{0}=0$ and the positive equilibrium, $x_{*}$. Recall that the characteristic equation of the linearization at $x_{0}=0$ is $\lambda+\rho=2\rho e^{-\lambda},$ which has a single positive root $\lambda_{0}$. The other roots form a sequence of complex conjugate pairs $(\lambda_{j},\bar{\lambda}_{j})$ with $\Re\lambda_{j+1}<\Re\lambda_{j}<\lambda_{0}$ for all integers $j\geq 1$. We refer to Chapter XI of dvvw for a complete analysis of the characteristic equation of the form $z-\alpha-\beta e^{-z}$, in particular Fig. XI.1. in dvvw which summarizes its properties. Our equation is a special case with $\alpha=-\rho$, $\beta=2\rho$.

There is a $\gamma>0$ such that $\lambda_{0}>\gamma>\Re\lambda_{j}$ for all $j\geq 1$, and there is a $k\in{\mathbb{N}}$ such that $\Re\lambda_{k}>0$ and $\Re\lambda_{k+1}\leq 0$, where the point $(-\rho,2\rho)$ lies between the curves $C_{k}^{-}$ and $C_{k+1}^{-}$ on the $(\alpha,\beta)$-plane, given by

[TABLE]

The intersections of $C_{j}^{-}$ and the half-line $(-\rho,2\rho)$ ($\rho>0$) are determined by $\cos\nu=1/2$, thus $\nu=2j\pi-\pi/3$, $\sin\nu=-\sqrt{3}/2$ and the intersection is given by $\rho=\rho_{j}=(2j\pi-\pi/3)/\sqrt{3}.$ Hence, $k$ is either zero or the largest positive integer such that $\rho>\rho_{k}$.

Considering the leading real root, the corresponding eigenfunction is given by $\chi_{0}(s):=e^{\lambda_{0}s}$, $s\in[-1,0].$ The phase space $C$ can be decomposed as $C=P_{0}\oplus P_{k}\oplus Q$, where the function $\chi_{0}\in C$ spans the linear eigenspace $P_{0}:=\{c\chi_{0}\,:\,c\in{\mathbb{R}}\}$. $P_{k}$ is a $2k$-dimensional eigenspace corresponding to $\{\lambda_{j}:j=1,\ldots,k\}$. If $\lambda_{0}$ is the only eigenvalue with positive real part then $P_{k}=\emptyset$. There is an $\varepsilon>0$ such that $\Re\lambda_{k}>\varepsilon$. $Q$ corresponds to the remaining part of the spectrum.

There exist open neighborhoods $N_{0}$, $N_{k}$, $M$ in $P_{0}$, $P_{k}$, $Q$, respectively, and $C^{1}$-maps $w_{0}:N_{0}\to P_{k}\oplus Q$, $w_{u}:N_{0}\oplus N_{k}\to Q$ with range in $N_{k}\oplus M$ and $M$, respectively, such that $w_{0}(0)=0$, $Dw_{0}(0)=0$, $w_{u}(0)=0$, $Dw_{u}(0)=0$; where $D$ denotes the Fréchet-derivative. Then, the $\gamma$-unstable set of the equilibrium $x_{0}=0$, namely

[TABLE]

coincides with the graph

[TABLE]

Furthermore, the unstable set of the equilibrium $x_{0}=0$

[TABLE]

coincides with the graph

[TABLE]

For the details see fhwmanifold ; kww .

For any $\phi\in N_{0}$, there is a $c\in{\mathbb{R}}$ such that $\phi=c\chi_{0}$. We have $||\chi_{0}||=1$ and $\chi_{0}(t)\geq e^{-\lambda_{0}}>0$ for all $t\in[-1,0]$. It follows from $Dw(0)=0$ that

[TABLE]

which means

[TABLE]

Therefore, there exists a $c_{0}>0$ such that $||w(c\chi_{0})||/|c|<e^{-\lambda}/2$ whenever $c\in(0,c_{0})$, or equivalently $||w(c\chi_{0})||<ce^{-\lambda}/2$. We may assume that $c_{0}<1/(2(1+\rho))$. Then

[TABLE]

and

[TABLE]

thus

[TABLE]

The unstable set $\mathcal{W}_{0}(0)$ intersects $\Omega$, and for any function $\phi_{c}$ of this intersection, there is a complete solution $x(t):{\mathbb{R}}\to{\mathbb{R}}^{+}_{0}$ such that $x_{0}=\phi_{c}$ and $x_{t}\to 0_{*}$ as $t\to-\infty$. By the negative invariance of $\mathcal{W}_{0}(0)$, the trajectory is completely in $\Omega$, hence we have proved the existence of a heteroclinic orbit connecting the two equilibria.

For uniqueness, consider the unstable manifold $\mathcal{W}_{u}(0)$. Any solution on the set $\mathcal{W}_{u}(0)\backslash\mathcal{W}_{0}(0)$ oscillates about [math] as $t\to-\infty$, hence can not be a positive heteroclinic solution. Denote the unique heteroclinic orbit in $\Omega$ by $\mathcal{H}$. We claim that $\mathcal{H}$ and the two equilibria constitute the global attractor, otherwise there are further complete orbits in $\Omega$. Assume there is a $\psi\in\Omega$ which is on such a complete orbit $\Gamma\subset\Omega$, so $\psi\neq 0$, $\psi\neq x_{*}$, and $\psi\notin\mathcal{H}$. Then its $\alpha$-limit set $\alpha(\psi)$ and $\omega$-limit set $\omega(\psi)$ are non-empty, compact and invariant. In particular, by Theorem 4.3, $\omega(\psi)=\{x_{*}\}$. If $\alpha{(\psi)}=\{x_{0}\}$, then $\Gamma$ is a connecting orbit coinciding with $\mathcal{H}$. Hence, we may assume that there is a $\phi\neq 0$ such that $\phi\in\alpha(\psi)$. Then $x_{t}^{\phi}\to x_{*}$ as $t\to\infty$, so $x_{*}\in\alpha(\psi)$ as well. But this contradicts the stability of $x_{*}$. Thus, the global attractor does not contain any other orbit, and $\mathcal{A}=\{x_{0}\}\cup\{x_{*}\}\cup\mathcal{H}$.

5 Numerical simulations and metastability

Representative numerical simulations of solutions of equation (13) are plotted in Fig. 1, with initial function $\phi(t)=0.005(\cos(10t)+1)$ and parameter values $\rho=50,120,190$. As we can see, in a short time interval of length 10, the solutions appear to be periodic, with their periods approximately equal to the delay, despite the fact that we proved in Theorem 4.3 that all non-trivial solutions converge to the positive equilibrium. Indeed, viewing the solution on a larger time scale in Fig. 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 holmes , 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, erneux ; grotta . Long lasting transient oscillations have been reported in morozov for a scalar differential equation with a constant delay. It is interesting to see that equation (13) can exhibit rapid convergence for small $\rho$ and oscillatory patterns of various shapes for larger $\rho$, 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 $\rho^{2}$. For large $\rho$, we use the notation $\varepsilon=1/\rho$, then we have

[TABLE]

Setting $\varepsilon=0$, we obtain the singular perturbation

[TABLE]

Besides $x\equiv 0$, any periodic function with unit period that is on average zero is a solution of this equation, which may be a hint as to why we see long-lasting transient oscillations with a variety of shapes.

In addition, we can look at the singular perturbation of the characteristic equation: equation (39) can be written as

[TABLE]

Multiplying by $(\rho+1)/\rho^{2}$ we have

[TABLE]

which is

[TABLE]

with the notation $\varepsilon=1/\rho$. Setting $\varepsilon=0$, we obtain the singular perturbation

[TABLE]

The roots of this equation satisfy either $\lambda=-1$ or $e^{-\lambda}=1$, having the purely imaginary roots $i=2k\pi$, $k\in\mathbb{Z}$. Indeed, as we increase $\rho$, we can see that all complex roots line up along the imaginary axis, see Fig. 4.

6 Discussion

The delayed logistic equation has received much attention in past decades in the analysis of nonlinear delay differential equations. However, its biological validity has been questioned, despite the fact that it was introduced by Hutchinson to explain observations of oscillatory behaviour in ecological systems. Here, we introduced a new logistic-type delay differential equation, derived from the go-or-grow hypothesis which has been observed for some types of cancer cells. The equation naturally includes a product of delayed and non-delayed terms, as well as both discrete and distributed delays. Detailed mathematical analysis reveals that all non-zero solutions converge to the stable positive equilibrium. For an alternative formulation of the delayed logistic equation arino , the authors also concluded that, as opposed to Hutchinson’s equation, sustained oscillations were not possible and solutions settled at an equilibrium. In some sense the dynamical behaviour of our model is in-between these two extremes. While we have proved that all positive solutions are attracted to the positive equilibrium, for a range of parameters the convergence is so slow that, on a biologically realistic time scale, solutions may appear periodic. These long-lasting transient patterns can have various shapes, as shown 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 equilibria and a unique connecting orbit. This heteroclinic solution is depicted in Fig. 5 when $\rho=20$. In this case, the leading real eigenvalue at zero is $\lambda\approx 0.66$, while the leading pair of eigenvalues at the positive equilibrium are $\lambda\approx-0.04\pm 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 $\rho$ is the product of the time delay, $\tau$, 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 faria2 ). 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.

Acknowledgments

GR was supported by NKFI FK 124016 and MSCA-IF 748193. REB is a Royal Society Wolfson Research Merit Award holder and would like to thank the Leverhulme Trust for a Research Fellowship.