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 non-trivial solutions converge to the positive equilibrium. The main tools of the proof rely on persistence theory, comparison principles and an L2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$L^2$$\end{document}-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) = r N(t) (1 − N (t)/K ) was proposed by Verhulst in 1838 to resolve the Malthusian dilemma of unbounded growth (Bacaër 2011). 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) 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 1955Wright (1955, 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 (May 1973). However, the use of Hutchinson's delayed logistic equation received criticisms from biological modellers due to the lack of a mechanistic biological derivation (Geritz and Kisdi 2012), 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. (1980) proposed the equation 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. (2006) derived an alternative formulation of the delayed logistic equation which has recently been extended to distributed delays (Lin et al. 2018). 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 (Lindström 2017) 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 et al. (2006), 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 and Mallet-Paret 1977;Faria 2006;Hassard et al. 1981;Nussbaum 1975), asymptotic analysis Fowler (1982), 3/2-type stability criteria (Ivanov et al. 2002;Wright 1955) and the study of slowly oscillatory solutions (Lessard 2010). 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 (Bánhelyi et al. 2014;van den Berg and Jaquette 2018). Moreover, many variants and generalizations of the delayed logistic equation have been studied in the mathematical literature, for instance equation including positive instantaneous feedback (Győri et al. 2018), nonautonomous terms (Faria 2006), neutral terms (Gopalsamy and Zhang 1988), spatial diffusion (Zou 2002) and multiple delays (Yan and Shi 2017). Further examples are discussed in Kuang (1993), Liz (2014), Ruan (2006).
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 Sect. 2, we study its basic mathematical properties in Sect. 3, such as well-posedness, existence and stability of equilibria. The global dynamics is fully described in Sect. 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 Sect. 5. Finally, we summarize our findings and discuss their implications for future research in Sect. 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 aggressive cancers (Noren et al. 2016). Our model is designed to capture the "go or grow" hypothesis frequently acknowledged in the cancer literature (Farin et al. 2006;Giese et al. 2003), 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. We assume that daughter agents are initially motile. After attempting proliferation, proliferative agents switch back to being motile. We assume that agents do not die, reflecting that malignant cells tend to evade apoptosis. The initial agent distribution is, on average, spatially uniform and is achieved by populating lattice sites uniformly at a given probability.
Following arguments presented in Baker and Simpson (2010) 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. This mean-field model can be obtained directly as well by the following biological modelling arguments. The equations encode the model assumption that motile agents become proliferative with per capita rate r and stay in this proliferative phase for time τ , this is reflected by (6). Once the cell cycle is completed, proliferative cells switch back to being motile cells again, that is the middle term rm(t − τ ) in equation (5). The last term of Eq. (5) is obtained as follows. As the time needed for cell division elapsed, a daughter agent is placed into a lattice site, but only if that site is empty, otherwise the proliferation event is cancelled. The probability that site is empty at time t (in the mean-field limit) is (K − p(t) − m(t)) /K , which is the number of empty sites (K − p(t) − m(t)) divided by the total size K of the cell space. At time t, such attempts of placing daughter cells occur with rate rm(t − τ ), since this is the rate cells converted from motile to proliferative τ time ago. Hence, successfully completed cell proliferations occur at rate , which is the increase in the m population due to proliferation, since we assumed that daughter cells are initially motile. Let us note that since we assumed that there is no cell death in the model, there are no mortality terms in (5)-(6), and the term representing switching back from p to m is exactly rm(t − τ ).
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 proliferative subpopulation in the time interval [t − τ, t]. 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-dimensionalization, we also rescale time by lettingt = t/τ and we can, again, drop the tildes, to arrive at the new delayed logistic equation:

Basic Mathematical Properties
With the notation ρ = r τ > 0, we now analyze the scalar differential equation with discrete and distributed delays 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 (9) is of the form Then, the map f : C → R is given by and the initial data are specified by
Lemma 3.1 For every φ ∈ C, there exists a unique solution of the initial value problem (9), (11) defined on an interval [−1, A) for some 0 < A ≤ ∞ that depends continuously on the initial data.
Proof Recall the standard existence and uniqueness theorem for functional differential equations (Smith 2011). 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 Hence, L = 3ρ + M 2ρ 2 + 2ρ is a valid Lipschitz constant. Existence, uniqueness and continuous dependence then follow from the general theory, c.f. Theorem 3.7 of Smith (2011) and (Kuang 1993, Chapter 2).
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 For , and hence, 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 (9), (11) is that x t ≥ 0 and x(t) = 0 implies f (x t ) ≥ 0 (Smith 2011, Theorem 3.4). This clearly holds as for obtain 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 (9) has two steady states,x = 0 which is always unstable, and x * = 1/(ρ + 1) which is always locally asymptotically stable.
Linearization of Eq. (9) around x = 0 gives with the characteristic equation which always has a positive real root, and hence, 0 is unstable.
Next, we look at the stability of the equilibrium x * . Using the definition in Eq. (10), we have which has linear part Clearly, lim ||φ||→0 |g(φ)|/||φ|| = 0, and 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 Eq. (15), 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 Smith (2011). 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 (15)).
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.

Persistence
Theorem 4.1 Equation (9) is strongly uniformly persistent; that is, there exists a δ > 0, independent of the initial data, such that for any solution then Eq. (9) can be written as From the variation in constants formula, for any σ > ω > 0 we find We use the notation x ∞ = lim inf t→∞ x(t), and x φ denotes a specific solution with 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) such that q 4 > 1/2 (any q > 0.841 is just fine). Then, there is a T n → ∞, such that For the particular case x = x φ n , σ = t n and ω = T n , the variation in constants formula gives Using the integral mean value theorem, there is an η n ∈ [T n , t n ], such that Then, we have the relation Since x φ n (t n ) → 0, e −ρ(t n −T n ) x φ n (T n ) → 0 and 1 − e T n −t n → 1 as n → ∞, necessarily g x φ n η n → 0 as well as n → ∞. Since x φ n η n ∈ Ω, we have g x φ n η n ≥ ρx φ n η n (−1), and hence, x φ n (η n − 1) → 0 as well as n → ∞. Next, we claim that From the inequality x (t) ≥ −ρx(t), we find that x φ n (η n + s) ≤ e ρ x φ n (η n − 1) for s ∈ [−1, 0]. Hence, and the right-hand side is already shown to converge to zero, and thus, (20) holds. For sufficiently large n, we have the estimate g x φ n η n ≥ 2qρx φ n η n (−1) ≥ 2qρqx φ n ∞ (where the first inequality follows from (18) and (20)) and 1 − e ρ(T n −t n ) > q, therefore, from equation (19) we obtain We find 1 ≥ 2q 4 , which contradicts the choice of q. Proof This follows from the standard comparison principle, since for large t we have 1 > x(t − 1) > δ/2 > 0, and by Eq. (13), we see that θ(t) is squeezed between θ (t) andθ(t), which are the solutions of

Global Asymptotic Stability
where a 0 , b 0 ∈ R, τ > 0 are constants and a, b : [0, ∞) → R are continuous functions. Let the following assumptions be satisfied: (i) λ = a 0 + b 0 e −λτ has a unique root λ 0 with the largest real part, and this root λ 0 is real and simple; Then, every solution x(t) satisfies

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 Eq. (14). 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 Eq. (21) as follows, since w is exponentially convergent.
(i) λ = −ρ + ρe −λ has a unique root λ 0 with the largest real part, and this root λ 0 is real and simple (in fact λ 0 = 0). (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 (21), we find or (iii) The estimate in Eq. (24) shows that (iv) From the estimate in Eq. (24), we also find that As a result, and all the conditions of Theorem A hold. Thus, every solution z(t) satisfies Now notice that a solution x(t) of Eq. (9) with initial function ψ is also a solution of Eq. (22), and hence, it converges to a constant. In view of Theorem 4.1, this constant can only be the positive equilibrium.

The Global Attractor
Denote by Φ Ω the semiflow on Ω generated by (9). We say that a subset of Ω is the global attractor of Φ Ω , if it is compact, invariant and attracts each bounded subset of Ω.

Theorem 4.4 The global attractor A consists of the two equilibria and a heteroclinic orbit connecting the two.
Proof First, we demonstrate that the global attractor A of Φ Ω exists. From the positive invariance of the bounded set Ω, this semiflow generated by the equation on Ω is point dissipative. It follows from the Arzelà-Ascoli theorem that the solution operators Φ Ω t are completely continuous for t ≥ 1; then, by applying Theorem 3.4.8 of Hale (1988), 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 Ω connectingx = 0 and the positive equilibrium, x * . Recall that the characteristic equation of the linearization atx = 0 is λ + ρ = 2ρe −λ , which has a single positive root λ 0 . The other roots form a sequence of complex conjugate pairs (λ j ,λ j ) with λ j+1 < λ j < λ 0 for all integers j ≥ 1. We refer to Chapter XI of Diekmann et al. (2012) for a complete analysis of the characteristic equation of the form z − α − βe −z , in particular Fig. XI.1. in Diekmann et al. (2012) which summarizes its properties. Our equation is a special case with α = −ρ, β = 2ρ.
Considering the leading real root, the corresponding eigenfunction is given by χ 0 (s) := e λ 0 s , s ∈ [−1, 0]. The phase space C can be decomposed as C = P 0 ⊕P k ⊕Q, 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.
For uniqueness, consider the unstable manifold W u (0). Any solution on the set W u (0)\W 0 (0) oscillates about 0 as t → −∞, and hence, it cannot 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}, 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} ∪ {x * } ∪ H.

Numerical Simulations and Metastability
Representative numerical simulations of solutions of Eq. (9) are plotted in Fig. 1, with initial function φ(t) = 0.005(cos(10t) + 1) and parameter values ρ = 50, 120, 190. To avoid numerical issues with representing the integral term, for the numerical simulations we used the system version of the model that includes only constant delay. As we can see in Fig. 1, 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 2012), 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 (2009), Grotta-Ragazzo et al. (2010. Long-lasting transient oscillations have been reported in Morozov et al. (2016) for a scalar differential equation with a constant delay. It is interesting to see that Eq. (9) 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. Nevertheless, for such large ρ, these oscillations in x will be very small compared to the total cell density which converges to one.
Since the long-lasting transient oscillations emerge for large ρ, to gain some understanding of this phenomenon, we look at our model from the point of view of singular perturbations, letting ρ → ∞. In Eq. (9), let us divide by ρ 2 . For large ρ, we use the notation ε = 1/ρ, and then, we have x(s) ds. (28) Setting ε = 0, we obtain the singular perturbation Besides x ≡ 0, any periodic function with unit period that is on average zero is a solution of this equation. While the singular limit is an equation from a different class, and its periodic solutions do not have non-negativity that we require for the differential equation, the behaviour of this singular limit may be a hint as to why we see long-lasting transient oscillations with a variety of shapes for large ρ.
In addition, we can look at the singular perturbation of the characteristic equation: Eq. (17) can be written as Multiplying by (ρ + 1)/ρ 2 , we have which is with the notation ε = 1/ρ. Setting ε = 0, we obtain the singular perturbation

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 nonzero solutions converge to the stable positive equilibrium. For an alternative formulation of the delayed logistic equation (Arino et al. 2006), 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 ρ = 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-orgrow 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 (9), 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 Eq. (9)), namely the cell motility rate. An important future goal is to understand how the speed of, for example, 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, Faria and Trofimchuk 2006). 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.