A novel mathematical model of heterogeneous cell proliferation
Sean T. Vittadello, Scott W. McCue, Gency Gunasingh, Nikolas K. Haass, Matthew J. Simpson
AA novel mathematical model of heterogeneous cell proliferation
Sean T. Vittadello † , Scott W. McCue , Gency Gunasingh , Nikolas K. Haass , Matthew J. Simpson Abstract
We present a novel mathematical model of heterogeneous cell proliferation where the totalpopulation consists of a subpopulation of slow-proliferating cells and a subpopulation of fast-proliferating cells. The model incorporates two cellular processes, asymmetric cell divisionand induced switching between proliferative states, which are important determinants for theheterogeneity of a cell population. As motivation for our model we provide experimental datathat illustrate the induced-switching process. Our model consists of a system of two coupleddelay differential equations with distributed time delays and the cell densities as functionsof time. The distributed delays are bounded and allow for the choice of delay kernel. Weanalyse the model and prove the non-negativity and boundedness of solutions, the existenceand uniqueness of solutions, and the local stability characteristics of the equilibrium points.We find that the parameters for induced switching are bifurcation parameters and thereforedetermine the long-term behaviour of the model. Numerical simulations illustrate and supportthe theoretical findings, and demonstrate the primary importance of transient dynamics forunderstanding the evolution of many experimental cell populations.
Key words and phrases : Delay differential equation, Integro-differential equation, Distributed delay, Asymmetriccell division, Induced switching, Transient dynamics School of Mathematical Sciences, Queensland University of Technology, Brisbane, Queensland, Australia. Current address:
Melbourne Integrative Genomics, School of BioSciences, The University of Melbourne, Parkville,Victoria, Australia. † Corresponding author: [email protected] The University of Queensland Diamantina Institute, The University of Queensland, Brisbane, Queensland, Australia.
July 27, 2020 a r X i v : . [ q - b i o . CB ] J u l Introduction
Cell proliferation is the fundamental function of the cell cycle [Matson and Cook, 2017], which is a complexprocess regulated by both intracellular signals and the extracellular environment [Zhu and Thompson,2019]. Such complexity necessitates that mathematical models of cell proliferation are often restricted todetails that are most pertinent to the experimental situation under consideration. The main requirementis that the model must account for progression through the cell cycle in a manner relevant to the cellpopulation and the surrounding environment. Despite all of the underlying complexity the cell cycle hastwo basic fates, either progression or arrest [Matson and Cook, 2017]. These two cellular fates form thebasis of many mathematical models of cell proliferation in the literature, typically based on exponentialgrowth [Lebowitz and Rubinow, 1974, Webb, 1986, Swanson et al., 2003, Sarapata and de Pillis, 2014] orlogistic growth [Sherratt and Murray, 1990, Maini et al., 2004, Cai et al., 2007, Byrne and Drasdo, 2009,Scott et al., 2013, Sarapata and de Pillis, 2014]. Exponential growth explicitly accounts for progressiononly, while logistic growth accounts for progression and density-dependent arrest, which can result fromcontact inhibition [Pavel et al., 2018].An important detail of the cell cycle not explicitly accounted for in exponential and logistic growthmodels is the duration of the cell cycle, which is always nonzero and exhibits considerable variationbetween different cell types and different extracellular environments [Weber et al., 2014, Chao et al.,2019, Vittadello et al., 2019, Vittadello et al., 2020]. From a modelling perspective the cell cycle durationis a positive time delay between two sequential cell proliferation events. There are two main typesof models which incorporate time delays: one involves functional differential equations [Mackey andRudnicki, 1994, Byrne, 1997, Baker et al., 1997, Baker et al., 1998, Villasana and Radunskaya, 2003, Gettoand Waurick, 2016,Getto et al., 2019,Cassidy and Humphries, 2020], of which delay differential equationsare a specific type; and multi-stage models [Yates et al., 2017, Simpson et al., 2018, Vittadello et al.,2018, Vittadello et al., 2019, Gavagnin et al., 2019]. Models incorporating time delays are consistent withthe kinetics of cell proliferation, and can result in a better qualitative and quantitative fit of the modelto experimental data [Baker et al., 1998]. The inclusion of a time delay must be based on whether theimproved model fit outweighs the increase in model complexity arising from additional parameters and,for functional differential equations, an infinite-dimensional state space. Models with time delays areparticularly relevant when the transient dynamics of a cell population are of interest, especially whenmodelling slow-proliferating cells. We briefly note that age-structured models [Arino, 1995, Gabriel et al.,2012, Billy et al., 2014, Clairambault and Fercoq, 2016, Cassidy et al., 2019], which are related to delaydifferential equations, provide another approach to incorporating realistic cell cycle durations into models1f cell population growth.In this article we introduce a delay differential equation model for cell proliferation in which the cellpopulation consists of a slow-proliferating subpopulation and a fast-proliferating subpopulation. The cellscan switch between the proliferative states of slow and fast proliferation through two cellular processes:asymmetric cell division and induced switching of proliferative states by surrounding cells. Our modelis motivated by the proliferative heterogeneity , with respect to cell cycle duration, of tumours, which areoften composed of a large proportion of fast-proliferating cells and a small proportion of slow-proliferatingcells which can repopulate the fast-proliferating subpopulation [Perego et al., 2018, Vallette et al., 2019].The slow-proliferating subpopulation is sometimes considered to be quiescent, or arrested, however it ispossible that this subpopulation is actually in a very-slow-proliferating state [Moore and Lyle, 2011, Ahnet al., 2017]. Experimental studies have found slow-proliferating cells with cell cycle durations greaterthan four weeks, whereas the predominant fast-proliferating cells have cell cycle durations of around 48hours [Roesch et al., 2010].In the literature there are various mathematical models that consider proliferative heterogeneity. Somemodels account for one proliferating subpopulation [Cassidy and Humphries, 2020], which may undergoasymmetric division [Arino and Kimmel, 1989, Greene et al., 2015], while the other subpopulations arequiescent or differentiated. Other models consider subpopulations with different proliferative states with-out any cells switching between the subpopulations [Jin et al., 2018]. Our model, and our mathematicalanalysis of the model, are novel in several ways: (1) for each subpopulation we model a distribution ofcell cycle durations using a distributed delay with an arbitrary delay kernel on a bounded interval, whichallows us to freely choose an appropriate proliferative state for each subpopulation; (2) cells can switch be-tween the slow- and fast-proliferating subpopulations through two important processes, either during celldivision or induced by surrounding cells; (3) we provide formal proofs of existence, uniqueness, non-neg-ativity, and boundedness of the solutions for our model under appropriate initial conditions; (4) the localstability of all equilibrium points is characterised and bifurcation parameters identified, involving theanalysis of an interesting transcendental characteristic equation; (5) numerical simulations are providedwhich illustrate and support the theoretical results, and demonstrate the importance of considering thetransient dynamics of experimental cell populations.The remainder of this article is organised as follows. In Section 2 we discuss the biological andmathematical motivations for our model, which we then present in Section 3. Our main analyticalresults are in Section 4 in the form of three theorems: Theorem 2 for non-negativity and boundedness ofsolutions, Theorem 4 for the existence and uniqueness of solutions, and Theorem 5 for the local stabilityof the equilibrium points. Some examples of numerical simulations of our model are provided in Section 5,2llustrating the long-term dynamics described in Theorem 5, and demonstrating the importance of thetransient dynamics. Finally, in Section 6 we summarise our results, discuss the utility of our model todescribe experimental scenarios, and note some possibilities for future work.
In this section we discuss the biological and mathematical considerations that motivate the developmentof our model.
The eukaryotic cell cycle (Figure 1) is a sequence of four phases, namely gap 1 (G1), synthesis (S),gap 2 (G2) and mitosis (M). The primary function of the cell cycle is the replication of cellular DNA M G1 SG2M
Figure 1:
Schematic of the eukaryotic cell cycle, indicating the colour of flu-orescent ubiquitination-based cell cycle indicator (FUCCI), see [Sakaue-Sawanoet al., 2008], in each phase. During very early G1 phase there is no fluorescence asboth FUCCI reporters are downregulated. In G1 phase the red FUCCI reporteris upregulated and red fluorescence is observed. During the transition from G1to S phase, called early S, both the red and green FUCCI reporters are upreg-ulated producing yellow. Through S/G2/M phase the red FUCCI reporter isdownregulated and only the green FUCCI reporter is upregulated so that greenfluorescence is observed during S phase, followed by the division of the replicated chromosomes and cytoplasm into two daughtercells during M phase [Vermeulen et al., 2003]. Progression through the cell cycle is tightly regulated innormal cells, which are subject to density-dependent contact inhibition producing reversible cell-cyclearrest [McClatchey and Yap, 2012, Puliafito et al., 2012]. In cancer cells, however, cell cycle regulation isgenerally lost [Hanahan and Weinberg, 2011] resulting in cell populations with proliferative heterogeneity,as exemplified by tumours of solid cancers [Roesch et al., 2010, Perego et al., 2018]. In particular a smallsubpopulation of slow-proliferating cells is often present in tumours, and this subpopulation tends tosurvive anticancer drug treatment and can maintain the tumour by repopulating the fast-proliferatingsubpopulation [Perego et al., 2018, Vallette et al., 2019].Experimental studies have revealed the highly dynamic nature of intratumoural heterogeneity, whichcan cause adverse outcomes from cancer therapy, notably drug resistance [Haass et al., 2014, Beaumontet al., 2016, Ahmed and Haass, 2018, Gallaher et al., 2019]. Therefore the nonequilibrium, or transient,state of a tumour tends to be of greater relevance than the equilibrium states. The main purpose of ourmodel is to provide insight into the transient dynamics of intratumoural heterogeneity, specifically withregard to cells switching their proliferative states through cellular mechanisms.3he range of mechanisms leading to proliferative heterogeneity in cancer cell populations are notcompletely understood, although asymmetric cell division is known to be partly responsible [Bajaj et al.,2015, Dey-Guha et al., 2015]. Asymmetric cell division is a normal process of stem cell proliferation, re-quired for development and the maintenance of tissue homeostasis, whereby a stem cell divides to produceone daughter stem cell, called self renewal , and a second daughter cell that will undergo differentiation.In contrast, symmetric division of a stem cell produces either two daughter stem cells or two daughtercells that will both undergo differentiation [Bajaj et al., 2015]. It is known that cancer cells can utilise thepathway of asymmetric cell division to produce heterogeneous populations of cancer cells that supportsurvival of the cancer [Smalley and Herlyn, 2009, Dey-Guha et al., 2011, Dey-Guha et al., 2015].Another important mechanism contributing to proliferative heterogeneity in cancer cell populationsis cell-induced switching between the slow- and fast-proliferating states, which occurs through cell–cellsignalling and direct contact between cells [Nelson and Chen, 2002, West and Newton, 2019]. The pos-sibility that cells can switch their proliferative state, either through asymmetric cell division or inducedswitching by surrounding cells, means that the growth rate of a cell population is highly dependent onthe influence of each of these two processes.To illustrate the concept of induced switching we show in Figure 2 a series of our experimental imagesfrom a two-dimensional proliferation assay using FUCCI-C8161 melanoma cells [Haass et al., 2014,Spoerriet al., 2017]. See Electronic Supplementary Material for further details. While this experiment is not (a) t = 0 h t = 3 h t = 7 h t = 9 h3 2 1 3 2 13 2 1 (f) t = 14 h t = 17 h t = 19 h (b) (c) (d) (e) (g) (h)
Figure 2:
Experimental images from a proliferation assay using FUCCI-C8161 melanoma cells, illustrating thepossibility of induced switching between proliferative states, as discussed in the text. explicitly designed to study slow- and fast-proliferating subpopulations, it is possible that switching of4ell cycle speeds occurs and can be observed by careful inspection of the time-series images. For example,consider the cells labelled 1, 2 and 3 throughout the images. Cells 1 and 2 are daughter cells from the sameparent cell and are in an early stage of G1, while cell 3 is a daughter cell from a different parent cell and isin a later stage of G1 at time 0 hours (Figure 2(a)). At time 3 hours we observe cell 3 interact closely withcell 2 and not cell 1 (Figure 2(b)). At times 7 and 9 hours the three cells continue to progress throughthe cell cycle with no close interaction between cell 3 and cells 1 and 2 (Figure 2(c)–(d)). At 11 hours,cell 3 interacts closely with cell 2 again; cell 2 is in S/G2/M phase, which is further through the cell cyclethan cell 1 which is in eS phase (Figure 2(e)). At time 14 hours cell 3, which is still close to cell 2, is in Mphase and is undergoing mitotic rounding in preparation for cell division (Figure 2(f)). At 17 hours, cell3 undergoes division to produce two daughter cells, cell 2 is undergoing mitotic rounding in preparationfor division, whereas cell 1 is in an earlier stage of the cell cycle (Figure 2(g)). At time 19 hours cell 2has divided to produce two daughter cells, whereas cell 1 has only just undergone mitotic rounding inpreparation for division (Figure 2(h)). These experimental observations illustrate the possibility that cellscan switch between states of slow and fast proliferation, induced by surrounding cells. In summary, itseems plausible that cell 2 progresses through the cell cycle faster than cell 1 because of the interactionswith cell 3. Indeed, for a given cell line, G1 phase tends to have the most variable duration of thecell cycle phases [Chao et al., 2019], and cells 1 and 2 appear to progress through G1 at a similar rate(Figure 2(a)–(c)), so it is possible that asymmetric division does not account for the overall variation incell cycle duration between cells 1 and 2.
Delay differential equations are often used when the evolution of the process to be modelled depends onthe history of the process, represented as a time delay which may be discrete [Lu, 1991,Engelborghs et al.,2000,Sun, 2006], distributed [McCluskey, 2010,Khasawneh and Mann, 2011,Huang et al., 2016,Kaslik andNeamtu, 2018, Cassidy and Humphries, 2020], or, more generally, state-dependent [Getto and Waurick,2016]. We employ a system of two coupled nonlinear delay differential equations to model the transientdynamics of cell proliferation in a population consisting of slow- and fast-proliferating cells. The timedelays are distributed so that only cells of a certain age can proliferate according to an appropriateprobability distribution, or delay kernel , and our system is therefore of integro-differential type. Analternative to modelling cell proliferation with delay differential equations is to use a multi-stage model,however this is not suitable for modelling the cell proliferation scenario that we consider here. Indeed,cell cycle durations in the multi-stage model are hypoexponentially distributed, and we want to allowfor more general distributions of cell cycle duration. Further, the multi-stage model can be difficult to5arameterise due to the large number of stages and hence parameters required to represent the cell cycleas stages.
Distributed delays
Standard deterministic mathematical models of cell proliferation, such as expo-nential and logistic growth models, are based on cell cycle durations with an exponential distribution,which allows for a relatively large probability of arbitrarily small cell cycle durations. Experimental in-vestigations, however, suggest that the duration of the cell cycle, and in particular each cell cycle phase,is not exponentially distributed, rather the hypoexponential distribution is often found to be a reasonableapproximation [Weber et al., 2014, Yates et al., 2017, Vittadello et al., 2019, Chao et al., 2019, Gavagninet al., 2019]. Our experimental data support these observations (Figure S1, Electronic SupplementaryMaterial). Ordinary differential equations therefore tend to overestimate the cell population growth rate,and may not qualitatively and quantitatively represent the transient growth dynamics of cell populations,particularly for slow-proliferating cells.The distributions of cell cycle durations for all cell lines are naturally bounded, so unbounded distri-butions such as the hypoexponential distribution are unrealistic since they theoretically have a nonzeroprobability of arbitrarily large values. For this reason we consider only bounded distributions for theirgreater biological realism, and we otherwise allow complete generality for the distributions. Unboundeddistributions can be left- or right-truncated to form bounded distributions, and while our model assumesa left bound of zero, alternative left bounds are easily incorporated. Our distributed delays have thefollowing form. Let X denote either the slow-proliferating cells S or fast-proliferating cells F , then thedistributed delay X ( t ) is defined by X ( t ) = Z U X X ( t − z ) g X ( z ) d z, (1)where the upper limit of integration U X ∈ (0 , ∞ ) corresponds to the maximum possible duration of thecell cycle, and g X is the delay kernel which is normalised so that R U X g X ( z ) d z = 1. The delay kernel isthe probability density for the distribution of the cell cycle durations. X ( t ) is a weighted average overthe past population densities X ( t − z ), and corresponds to the subpopulation of cells at time t that areready to divide. Contact inhibition
Normal cells are subject to contact inhibition [Pavel et al., 2018], so we assumethat the growth rate of each subpopulation at time t has a logistic density dependence given by (cid:0) − P ( t ) /K (cid:1) , where P is the total population density and K is the carrying capacity density. For cancer cells,however, contact inhibition may be lost [Hanahan and Weinberg, 2011] resulting in density-independent6rowth, so in this case we could set the carrying capacity to infinity in the logistic-growth terms. Wecan keep the carrying capacity finite in the induced-switching terms, as increasing the cell density beyonda finite carrying capacity would, realistically, increase the probability of cell-induced switching betweenproliferative states due to surrounding cells. Proliferation, switching between proliferative states, and apoptosis
Given the cells that areable to undergo division based on the time delays and density constraints, the intrinsic growth rates0 < r S ≤ r F for slow- and fast-proliferating cells, respectively, determine the cells that are parent cellsand divide at time t . Parent cells can divide symmetrically, where the daughter cells have the sameproliferative state as the parent cell, or asymmetrically, where a daughter cell can have a proliferativestate different from the parent cell. When a subpopulation of slow-proliferating cells divides to producetwice as many daughter cells, the parameter α S determines the proportion of these daughter cells thatare also slow-proliferating cells, so the proportion 1 − α S of the daughter cells are fast-proliferating cells.Similarly, the parameter α F determines the proportion of daughter cells from fast-proliferating cells thatare also fast-proliferating cells, so the proportion 1 − α F of the daughter cells are slow-proliferating cells.We only consider asymmetric cell division that is self-renewing, that is α S , α F ∈ [ , α S , α F ∈ [0 , α S , α F ∈ [ , β S corresponds tothe per capita interaction strength of fast-proliferating cells to induce slow-proliferating cells to switchto fast proliferation. Similarly, the parameter β F corresponds to the per capita interaction strength ofslow-proliferating cells to induce fast-proliferating cells to switch to slow proliferation.Because we focus on growing tumours, for which proliferation outcompetes cell death, and mechanismsof switching between proliferative states, we simplify our model by not incorporating apoptosis.7 Mathematical model
The total cell population consists of the two subpopulations of slow-proliferating cells and fast-proliferatingcells, where P ( t ) ≥ S ( t ) ≥
0, and F ( t ) ≥ P ( t ) = S ( t ) + F ( t ).The model is d S ( t )d t = (2 α S − r S Z U S S ( t − z ) g S ( z ) d z | {z } slow-proliferating cells from asymmetricdivision of slow-proliferating cells (cid:18) − (cid:0) S ( t ) + F ( t ) (cid:1) K (cid:19)| {z } contact inhibitionof proliferation + 2(1 − α F ) r F Z U F F ( t − z ) g F ( z ) d z | {z } slow-proliferating cells from asymmetricdivision of fast-proliferating cells (cid:18) − (cid:0) S ( t ) + F ( t ) (cid:1) K (cid:19)| {z } contact inhibitionof proliferation − β S S ( t ) F ( t ) K | {z } induced switchingof slow-proliferating cellsto fast-proliferating cells + β F F ( t ) S ( t ) K | {z } induced switchingof fast-proliferating cellsto slow-proliferating cells , (2)d F ( t )d t = 2(1 − α S ) r S Z U S S ( t − z ) g S ( z ) d z | {z } fast-proliferating cells from asymmetricdivision of slow-proliferating cells (cid:18) − (cid:0) S ( t ) + F ( t ) (cid:1) K (cid:19)| {z } contact inhibitionof proliferation + (2 α F − r F Z U F F ( t − z ) g F ( z ) d z | {z } fast-proliferating cells from asymmetricdivision of fast-proliferating cells (cid:18) − (cid:0) S ( t ) + F ( t ) (cid:1) K (cid:19)| {z } contact inhibitionof proliferation + β S S ( t ) F ( t ) K | {z } induced switchingof slow-proliferating cellsto fast-proliferating cells − β F F ( t ) S ( t ) K | {z } induced switchingof fast-proliferating cellsto slow-proliferating cells , (3)where the parameters satisfyIntrinsic growth rates: r S , r F ∈ (0 , ∞ ) with r S ≤ r F , (4)Proportion of symmetric divisions: α S , α F ∈ [ , U S , U F ∈ (0 , ∞ ), (6)Induced switching rates: β S , β F ∈ [0 , ∞ ), (7)and the mean values of the delay kernels g S and g F satisfy Z U S z g S ( z ) d z ≥ Z U F z g F ( z ) d z, (8)8o that the mean cell cycle duration for the slow-proliferating cells is not smaller than the mean cell cycleduration for the fast-proliferating cells.As functional differential equations depend on the solution and perhaps derivatives of the solutionat past times it is necessary to specify a function for the initial condition, called the history function .Defining b U = max { U S , U F } , the history function φ = ( φ S , φ F ) for our model satisfies φ ∈ C ([ − b U , , R > ) , (9) φ S + φ F ∈ C ([ − b U , , (0 , K )) , (10)where C ([ − b U , , R > ) is the space of continuous functions on [ − b U ,
0] into R > and C ([ − b U , , (0 , K ))is the space of continuous functions on [ − b U ,
0] into (0 , K ). Note that state space is therefore aninfinite-dimensional function space. For bounded delays the state space is typically the Banach space C ([ − h, , R n ), for some h ∈ (0 , ∞ ), of continuous functions χ : [ − h, → R n on the closed interval [ − h, k·k defined by k χ k = sup { k χ ( t ) k | t ∈ [ − h, } , where k·k is the Euclideannorm on R n . In our case, state space C is defined by C = C ([ − b U , , R ) . (11)Finally, we note that adding Equations (2) and (3) givesd P ( t )d t = (cid:18) r S Z U S S ( t − z ) g S ( z ) d z + r F Z U F F ( t − z ) g F ( z ) d z (cid:19)(cid:18) − P ( t ) K (cid:19) , (12)so, from the perspective of the whole population, asymmetric division and induced switching have no neteffects. Moreover, if we consider the total population as composed of cells in the same proliferative statewith r S = r F , and g S ( z ) = g F ( z ) = δ ( z ) is the Dirac kernel for zero delay, then Equation (12) reduces tothe logistic growth model (Electronic Supplementary Material).9 Main results
We now discuss our analysis of Equations (2) and (3), namely non-negativity and boundedness of solutions,existence and uniqueness of solutions, and local stability analysis of the equilibrium points. A solutionfor Equations (2) and (3) means the following [Kuang, 1993, Smith, 2011]:
Definition . Consider the system of delay differential equations (2) and (3) with parameters,delay kernels, and history functions that satisfy (4)–(10). A solution for the system is a function (
S, F ) ∈ C ([ − b U , u ) , R ≥ ) for some u ∈ [0 , ∞ ] such that: • S and F are differentiable on (0 , u ) and right-differentiable at 0; • ( S, F ) satisfies Equations (2) and (3) for t ∈ [0 , u ).Additionally, ( S, F ) is a solution with initial condition φ = ( φ S , φ F ) ∈ C ([ − b U , , R > ) if • ( S, F ) | [ − b U, = φ . Since the dependent variables S ( t ) and F ( t ) in Equations (2) and (3) represent cell densities they mustassume non-negative values at all times. Further, the densities of normal cells are bounded above bythe carrying capacity density K arising from contact inhibition. Cancer cells typically have unregulatedgrowth due to the loss of contact inhibition, however continued growth depends on environmental con-ditions such as nutrient availability, so it is reasonable to assume that the density of cancer cells is alsobounded by a carrying capacity. Theorem 2 (Non-negativity and boundedness of solutions) . Let ( S, F ) be a solution for the system of delay differential equations (2) and (3) with parameters, delaykernels, and history functions that satisfy (4) – (10) . Then S ( t ) , F ( t ) ∈ [0 , K ] for all t > , therefore thesolutions S and F are non-negative and bounded. We give an elementary proof of this theorem as it facilitates understanding of the non-negativity andboundedness of the solutions. We first require a lemma.
Lemma 3.
Let ( S, F ) be a solution for the system of delay differential equations (2) and (3) with pa-rameters, delay kernels, and history functions that satisfy (4) – (10) . If there exists T > such that thedistributed delays satisfy the relations S ( t ) ≥ and F ( t ) ≥ for all t ∈ (0 , T ] then P ( t ) ∈ [0 , K ] for all t ∈ (0 , T ] . roof. If P ( t ) > K for some t ∈ (0 , T ] then, since P (0) < K by Equation (10), we may assumewithout loss of generality that d P ( t ) / d t | t = t >
0. Since Equation (12) gives d P ( t ) / d t | t = t ≤
0, we havea contradiction. It follows that P ( t ) ≤ K for all t ∈ (0 , T ]. Similarly, if P ( t ) < t ∈ (0 , T ]then, since P (0) > P ( t ) / d t | t = t < P ( t ) / d t | t = t ≥
0, we have a contradiction. It follows that P ( t ) ≥ t ∈ (0 , T ]. Theorem 2.
It suffices to prove that S ( t ) ≥ F ( t ) ≥ t >
0, for then it follows fromLemma 3 that S ( t ) ≤ K and F ( t ) ≤ K for all t >
0. Define t and t by t = inf { t > | S ( t ) < } and t = inf { t > | F ( t ) < } . We consider the infima in the extended real numbers, so that t , t ∈ (0 , + ∞ ],where either infimum is equal to + ∞ if the corresponding set is empty. The proof consists of four separatecases which are proved similarly. We demonstrate one case here, and provide the complete proof in theElectronic Supplementary Material.Case 1: Let β S − β F ≥ S ( t ) < t > t ∈ R . If t < t then choose t ∈ ( t , t ) such that S ( t ) <
0, d S ( t ) / d t | t = t <
0, and the delayssatisfy S ( t ) ≥ F ( t ) ≥
0. Then, since P ( t ) ≤ K by Lemma 3 and since F ( t ) ≥
0, Equation (2)gives d S ( t ) / d t | t = t ≥
0, a contradiction.If t ≥ t then, since F (0) > F ( t ) = 0, there exists t ∈ (0 , t ) such that d F ( t ) / d t | t = t <
0. Then,since P ( t ) ≤ K by Lemma 3, since the delays satisfy S ( t ) ≥ F ( t ) ≥
0, and since S ( t ), F ( t ) ≥ F ( t ) / d t | t = t ≥
0, a contradiction. We conclude that S ( t ) ≥ t > We begin by introducing some simplifying notation. For ρ in the state space C we denote the componentfunctions by ρ S and ρ F so that ρ = ( ρ S , ρ F ), and then define ρ S and ρ F by ρ S = Z U S ρ S ( − z ) g S ( z ) d z and ρ F = Z U F ρ F ( − z ) g F ( z ) d z . (13)Now we define f : C → R by f ( ρ ) = (cid:16) (2 α S − r S ρ S + 2(1 − α F ) r F ρ F (cid:17)(cid:18) − ( ρ S + ρ F )(0) K (cid:19) − ( β S − β F ) K ρ S (0) ρ F (0) (cid:16) − α S ) r S ρ S + (2 α F − r F ρ F (cid:17)(cid:18) − ( ρ S + ρ F )(0) K (cid:19) + ( β S − β F ) K ρ S (0) ρ F (0) . (14)11ote that f is continuous. If ( S, F ) is a solution for Equations (2) and (3), t ≥
0, and we define( S t , F t ) ∈ C ([ − b U , , R ) by ( S t ( r ) , F t ( r )) = ( S ( t + r ) , F ( t + r )) for r ∈ [ − b U , f (( S t , F t )) =[d S ( t ) / d t, d F ( t ) / d t ] T from (2) and (3). Theorem 4 (Existence and uniqueness of solutions) . Consider the system of delay differential equations (2) and (3) with parameters, delay kernels, and historyfunctions which satisfy (4) – (10) . Then there exists a unique solution ( S, F ) ∈ C ([ − b U , ∞ ) , R ≥ ) of (2) and (3) .Proof. Here we give an outline of the proof. The complete proof is provided in the Electronic Supplemen-tary Material. We first show that f defined in Equation (14) satisfies the following Lipschitz conditionon every bounded subset of C : for all M >
L > ρ , ψ ∈ C ([ − b U , , R )with k ρ k , k ψ k ≤ M we have k f ( ρ ) − f ( ψ ) k ≤ L k ρ − ψ k .To further simplify the notation in Equation (14) we define κ = (2 α S − r S , κ = 2(1 − α F ) r F , κ = 2(1 − α S ) r S , κ = (2 α F − r F , and κ = ( β S − β F ) /K . Now, f ( ρ ) − f ( ψ ) = (cid:0) κ ( ρ S − ψ S ) + κ ( ρ F − ψ F ) (cid:1)(cid:18) − ( ρ S + ρ F )(0) K (cid:19) + (cid:16) κ ψ S + κ ψ F (cid:17)(cid:18) ( ψ S − ρ S )(0) + ( ψ F − ρ F )(0) K (cid:19) + κ ( ψ S − ρ S )(0) ψ F (0) + κ ( ψ F − ρ F )(0) ρ S (0) (cid:0) κ ( ρ S − ψ S ) + κ ( ρ F − ψ F ) (cid:1)(cid:18) − ( ρ S + ρ F )(0) K (cid:19) + (cid:16) κ ψ S + κ ψ F (cid:17)(cid:18) ( ψ S − ρ S )(0) + ( ψ F − ρ F )(0) K (cid:19) + κ ( ρ S − ψ S )(0) ψ F (0) + κ ( ρ F − ψ F )(0) ρ S (0) so, using the triangle inequality, we obtain k f ( ρ ) − f ( ψ ) k ≤ (cid:18)q κ + κ + q κ + κ (cid:19)(cid:18) MK (cid:19) + 2 √ | κ | M ! k ρ − ψ k , so we can set L to be L = (cid:18)q κ + κ + q κ + κ (cid:19)(cid:18) MK (cid:19) + 2 √ | κ | M and then f satisfies the Lipschitz condition. Then [Smith, 2011, Page 32, Theorem 3.7] provides localexistence and uniqueness of solutions for the system (2) and (3). Since our solutions of interest arebounded by Theorem 2, it follows from [Smith, 2011, Page 37, Proposition 3.10] that the solutions are12ontinuable to all positive time. Here we consider the local stability analysis of the equilibrium points for the system in (2) and (3), andshow that β S and β F are bifurcation parameters with bifurcation point when β S = β F . We will prove thefollowing theorem. Theorem 5 (Local stability) . Consider the system of delay differential equations (2) and (3) with parameters, delay kernels, and historyfunctions that satisfy (4) – (10) . • When β S = β F the system has the three equilibrium points (0 , , (0 , K ) , and ( K, with the followingproperties: – (0 , is locally unstable. – If β S > β F then ( K, is locally unstable and (0 , K ) is locally stable. – If β S < β F then ( K, is locally stable and (0 , K ) is locally unstable. • When β S = β F the system has infinitely many equilibrium points corresponding to the line segmentjoining ( K, and (0 , K ) , all of which are locally stable.The parameters β S and β F are therefore bifurcation parameters. Note that, since equilibrium points for delay differential equations are functions in a Banach space, wehave (0 , , K ), ( K, ∈ C ([ − b U , ∞ ) , R ≥ ).The proof of Theorem 5 follows immediately from Propositions 6, 8, 9 and 10. We begin by non-dimensionalising Equations (2) and (3), and then linearising the non-dimensional system about the equi-librium points. Denoting the dimensionless variables with a caret, we define ˆ t = r F t , ˆ S (ˆ t ) = S ( t ) /K and ˆ F (ˆ t ) = F ( t ) /K . We also define the dimensionless parameters r = r S /r F and β = ( β S − β F ) /r F .Equations (2) and (3) then become, dropping the caret notation for simplicity,d S ( t )d t = (cid:18) (2 α S − r Z U S S ( t − r F z ) g S ( z ) d z + 2(1 − α F ) Z U F F ( t − r F z ) g F ( z ) d z (cid:19) × (cid:0) − S ( t ) − F ( t ) (cid:1) − βS ( t ) F ( t ) , (15)d F ( t )d t = (cid:18) − α S ) r Z U S S ( t − r F z ) g S ( z ) d z + (2 α F − Z U F F ( t − r F z ) g F ( z ) d z (cid:19) × (cid:0) − S ( t ) − F ( t ) (cid:1) + βS ( t ) F ( t ) . (16)13ince S and F are cell densities, hence non-negative, we only consider equilibrium points ( S * , F * ) ∈ C ([ − b U , ∞ ) , R ≥ ). To find the equilibrium points we substitute S = S * and F = F * into (15) and (16) togive 0 = (cid:0) (2 α S − rS * + 2(1 − α F ) F * (cid:1)(cid:0) − S * − F * (cid:1) − βS * F * , (17)0 = (cid:0) − α S ) rS * + (2 α F − F * (cid:1)(cid:0) − S * − F * (cid:1) + βS * F * , (18)hence ( S * , F * ) = (0 , ,
0) or (0 ,
1) when β = 0. When β = 0 the equilibrium points consist of the twolines ( S * , F * ) = ( u, − u ) for all u ∈ R and ( S * , F * ) = ( u, − ru ) for all u ∈ R , for which the non-negativepoints are ( S * , F * ) = ( u, − u ) for all u ∈ [0 ,
1] and ( S * , F * ) = (0 , S * , F * ) we linearise the system in Equa-tions (15) and (16) about each point. Defining x ( t ) = S ( t ) − S * and y ( t ) = F ( t ) − F * we obtain thelinearised system:d x ( t )d t = (cid:18) (2 α S − r (cid:16) Z U S x ( t − r F z ) g S ( z ) d z + S * (cid:17) + 2(1 − α F ) (cid:16) Z U F y ( t − r F z ) g F ( z ) d z + F * (cid:17)(cid:19) × (cid:0) − x ( t ) − y ( t ) − S * − F * (cid:1) − β (cid:0) x ( t ) + S * (cid:1)(cid:0) y ( t ) + F * (cid:1) ∼ (cid:18) (2 α S − r Z U S x ( t − r F z ) g S ( z ) d z + 2(1 − α F ) Z U F y ( t − r F z ) g F ( z ) d z (cid:19) × (cid:0) − S * − F * (cid:1) + (cid:16) (2 α S − rS * + 2(1 − α F ) F * (cid:17)(cid:0) − x ( t ) − y ( t ) − S * − F * (cid:1) − β (cid:0) x ( t ) F * + y ( t ) S * + S * F * (cid:1) , as x ( t ), y ( t ) → , (19)d y ( t )d t = (cid:18) − α S ) r (cid:16) Z U S x ( t − r F z ) g S ( z ) d z + S * (cid:17) + (2 α F − (cid:16) Z U F y ( t − r F z ) g F ( z ) d z + F * (cid:17)(cid:19) × (cid:0) − x ( t ) − y ( t ) − S * − F * (cid:1) + β (cid:0) x ( t ) + S * (cid:1)(cid:0) y ( t ) + F * (cid:1) ∼ (cid:18) − α S ) r Z U S x ( t − r F z ) g S ( z ) d z + (2 α F − Z U F y ( t − r F z ) g F ( z ) d z (cid:19) × (cid:0) − S * − F * (cid:1) + (cid:16) − α S ) rS * + (2 α F − F * (cid:17)(cid:0) − x ( t ) − y ( t ) − S * − F * (cid:1) + β (cid:0) x ( t ) F * + y ( t ) S * + S * F * (cid:1) , as x ( t ), y ( t ) → . (20)14y the Principle of Linearised Stability [Diekmann et al., 1995, Page 240, Theorem 6.8] it suffices toconsider the stability of the equilibrium points for the linearisation in Equations (19) and (20). Proposition 6 (Equilibrium point (0 , . Consider the system of delay differential equations (2) and (3) with parameters, delay kernels, and historyfunctions which satisfy (4) – (10) . Then (0 , is locally unstable. Proposition 6 follows immediately from Proposition 7, in which we analyse the transcendental charac-teristic equation associated with (0 ,
0) of the linearised system (19) and (20) to show that the characteristicequation has at least one zero in C with positive real part. While there are alternative methods for prov-ing Proposition 6 (Electronic Supplementary Material), we consider the direct approach of analysing theassociated transcendental characteristic equation to have mathematical relevance for other studies involv-ing delay differential equations. Indeed, transcendental characteristic equations are generally difficult toanalyse [Diekmann et al., 1995, Chapter XI], so new analysis of such equations is of mathematical interest.For ( S * , F * ) = (0 , x ( t )d t = (2 α S − r Z U S x ( t − r F z ) g S ( z ) d z + 2(1 − α F ) Z U F y ( t − r F z ) g F ( z ) d z, (21)d y ( t )d t = 2(1 − α S ) r Z U S x ( t − r F z ) g S ( z ) d z + (2 α F − Z U F y ( t − r F z ) g F ( z ) d z. (22)Equations (21) and (22) have a solution of the form x ( t ) y ( t ) = c c e λt , where c c ∈ C is nonzero and λ ∈ C , (23)so substitution gives e λt λ λ c c = (2 α S − r Z U S e − λr F z g S ( z ) d z − α F ) R U F e − λr F z g F ( z ) d z − α S ) r Z U S e − λr F z g S ( z ) d z (2 α F − R U F e − λr F z g F ( z ) d z c c e λt . (24)15o ensure that ( c , c ) | = 0 we must have the characteristic equation G ( λ ) = (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) (2 α S − r Z U S e − λr F z g S ( z ) d z − λ − α F ) R U F e − λr F z g F ( z ) d z − α S ) r Z U S e − λr F z g S ( z ) d z (2 α F − R U F e − λr F z g F ( z ) d z − λ (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) = λ − λ (cid:18) (2 α S − r Z U S e − λr F z g S ( z ) d z + (2 α F − Z U F e − λr F z g F ( z ) d z (cid:19) + (2 α S + 2 α F − r Z U S Z U F e − λr F ( z + v ) g F ( v ) g S ( z ) d v d z (25)equal to zero. The zeros of the transcendental equation G ( λ ) are the eigenvalues. Proposition 7 showsthat G ( λ ) has at least one zero with positive real part, so the equilibrium point (0 ,
0) is locally unstable.In the proof of Proposition 7 we consider three cases for G ( λ ) in Equation (25), depending on whether2 α S + 2 α F − α S + 2 α F − α S + 2 α F − α S − − − α F ). Referring to Equation (15), (2 α S − − − α F )is the difference between the proportion of slow-proliferating parent cells that produce slow-proliferatingdaughter cells beyond self renewal, and the proportion of fast-proliferating parent cells that produce slow-proliferating daughter cells, at a given time. A similar interpretation follows by referring to Equation (16)and noting that 2 α S + 2 α F − α F − − − α S ). When 2 α S + 2 α F − G ( λ ) has a zero in R which is positive. When 2 α S + 2 α F − h be a meromorphic functionon Ω with no zeros or poles on Γ. Then Cauchy’s argument principle is [Ahlfors, 1979, Page 152, Theorem18] 12 πi I Γ h ( λ ) h ( λ ) d λ = Z − P , (26)where Z is the number of zeros of h inside Γ and P is the number of poles of h inside Γ, includingmultiplicities.Now, let Γ be a piecewise differentiable closed curve in C with positive orientation that does not passthrough the point z . Then the index of z with respect to Γ, denoted Ind Γ ( z ), is defined by [Ahlfors,1979, Page 115] Ind Γ ( z ) = 12 πi I Γ d zz − z . (27)Note that Ind Γ ( z ) is also referred to as the winding number of Γ with respect to z . By substituting16 = h ( λ ) into Equation (26) and using Equation (27) we arrive at the standard observation Z − P = 12 πi I h (Γ) d zz = Ind h (Γ) (0) , (28)where the term on the right of the last equality is the winding number of the closed curve h (Γ) withrespect to the origin. Therefore, the number of zeros minus the number of poles of h inside Γ, includingmultiplicities, can be determined by calculating the winding number of the image h (Γ) with respect tothe origin.Note that G ( λ ) in Equation (25) is a holomorphic function, hence meromorphic, on C , and has nopoles. Therefore, the number of zeros of G ( λ ) inside a contour Γ which satisifes the conditions forEquations (26) and (27) is equal to the winding number of G (Γ) with respect to the origin. We will beconsidering rectangular contours with positive orientation in the right half-plane, and we need to knowthat G ( λ ) is not identically zero in the region bounded by the closed contour. Our contour will bound aninterval of the positive real axis arbitrarily close to, but excluding, the origin, and with arbitrary upperbound. For a contour which bounds sufficiently large positive real numbers, and by considering Re( G ( λ ))on the positive real axis, we can see that G ( λ ) is not identically zero in the region bounded by the contour.Furthermore, by [Rudin, 1986, Page 208, Theorem 10.18] it follows that the zeros of G ( λ ) are isolatedand countable, so we can always choose an appropriate rectangular contour which does not pass througha zero of G ( λ ). In the Electronic Supplementary Material we graphically illustrate our application ofCauchy’s argument principle. We now state and prove the required proposition. Proposition 7.
The transcendental characteristic equation G ( λ ) in Equation (25) has at least one zeroin C with positive real part.Proof. We consider three cases according to whether 2 α S + 2 α F − α S + 2 α F − < G (0) = (2 α S + 2 α F − r < λ →∞ G | R ( λ ) =lim λ →∞ λ = ∞ , it follows that G has a positive real zero by the intermediate value theorem.Case 2: If 2 α S + 2 α F − G factors as G ( λ ) = λH ( λ ), where H ( λ ) = λ − (cid:18) (2 α S − r Z U S e − λr F z g S ( z ) d z + (2 α F − Z U F e − λr F z g F ( z ) d z (cid:19) . (29)Since 2 α S + 2 α F − α S − < α F >
1, so 2 α S − ≥
0. Similarly, 2 α F − ≥ α S − α F − α S + 2 α F − α S + 2 α F − H (0) = − (cid:0) (2 α S − r + (2 α F − (cid:1) <
0. Since lim λ →∞ H | R ( λ ) = lim λ →∞ λ = ∞ ,17t follows that H has a positive real zero by the intermediate value theorem.Case 3: Suppose now that 2 α S + 2 α F − >
0. We employ Cauchy’s argument principle to show thatthe holomorphic function G ( λ ) has a zero with positive real part. Let Γ be the simple closed contour withpositive orientation in the right half-plane of the complex plane defined piecewise as follows:Γ : (cid:0) m (1 − t ) + ( N/ t (cid:1) − iN, ≤ t ≤ , (30)Γ : (cid:0) ( N/ − t ) + (3 N/ t (cid:1) − iN, ≤ t ≤ , (31)Γ : (3 N/
2) + i (cid:0) ( − N )(1 − t ) + N t (cid:1) , ≤ t ≤ , (32)Γ : (cid:0) (3 N/ − t ) + ( N/ t (cid:1) + iN, ≤ t ≤ , (33)Γ : (cid:0) ( N/ − t ) + mt (cid:1) + iN, ≤ t ≤ , (34)Γ : m + i (cid:0) N (1 − t ) + ( − N ) t (cid:1) , ≤ t ≤ , (35)where we fix m > N > S j =1 Γ j is rectangular, with vertices at m − iN , 3 N/ − iN , 3 N/ iN , and m + iN . Since the zeros of aholomorphic function that is not identically zero are isolated and countable, we can choose m arbitrarilysmall and N arbitrarily large while ensuring G ( λ ) is nonzero on Γ. Therefore, since G ( λ ) has no poles,the number of zeros of G ( λ ) inside Γ is equal to the index of the image contour G (Γ) with respect to theorigin, Ind G (Γ) (0). We now begin our application of the argument principle. Since we only need to showthe existence of one zero with a positive real part, it suffices to prove that Ind G (Γ) (0) ≥
1. Specifically,we show that the image contour G (Γ) crosses the positive real axis at least once in a counter-clockwisedirection while encircling the origin, and doesn’t cross the positive real axis in a clockwise direction.We will traverse Γ for one cycle in a counter-clockwise direction beginning with Γ , and determinewhen G (Γ) crosses the positive real axis. For this it is helpful to consider the real and imaginary parts of G ( λ ), so evaluating G ( λ ) at the arbitrary complex number λ = a + ib we haveRe( G ( λ )) = a − b − a (2 α S − r Z U S e − ar F z cos( br F z ) g S ( z ) d z − a (2 α F − Z U F e − ar F z cos( br F z ) g F ( z ) d z − b (2 α S − r Z U S e − ar F z sin( br F z ) g S ( z ) d z − b (2 α F − Z U F e − ar F z sin( br F z ) g F ( z ) d z + (2 α S + 2 α F − r Z U S Z U F e − ar F ( z + v ) cos( br F ( z + v )) g F ( v ) g S ( z ) d v d z , (36)18nd Im( G ( λ )) = 2 ab + a (2 α S − r Z U S e − ar F z sin( br F z ) g S ( z ) d z + a (2 α F − Z U F e − ar F z sin( br F z ) g F ( z ) d z − b (2 α S − r Z U S e − ar F z cos( br F z ) g S ( z ) d z − b (2 α F − Z U F e − ar F z cos( br F z ) g F ( z ) d z − (2 α S + 2 α F − r Z U S Z U F e − ar F ( z + v ) sin( br F ( z + v )) g F ( v ) g S ( z ) d v d z . (37)In the following argument we shall generally use that m is arbitrarily small and N is arbitrarily largewithout further comment.Consider G ( λ ) along Γ , where b = − N and a increases from m to N/
2. For sufficiently large N andfor all a ∈ [ m, N/ G ( λ )) < G ( λ )) is dominated by N . At the end of Γ and for sufficientlylarge N , Im( G ( λ )) < G ( λ )) is dominated by N . So G (Γ ) starts in the left half-plane and endsin the third quadrant.Consider G ( λ ) along Γ , where b = − N and a increases from N/ N/
2. At the end of Γ and forsufficiently large N , Re( G ( λ )) > G ( λ )) is dominated by N . For sufficiently large N and for all a ∈ [ N/ , N/ G ( λ )) < G ( λ )) is dominated by N . So G (Γ ) starts in the third quadrantand ends in the fourth quadrant.Consider G ( λ ) along Γ , where a = 3 N/ b increases from − N to N . For sufficiently large N , Re( G ( λ )) > G ( λ )) is dominated by N . At the end of Γ and for sufficiently large N ,Im( G ( λ )) > G ( λ )) is dominated by N . So G (Γ ) starts in the fourth quadrant and ends inthe first quadrant. The image contour G (Γ) has now crossed the positive real axis in a counter-clockwisedirection.Consider G ( λ ) along Γ , where b = N and a decreases from 3 N/ N/
2. At the end of Γ and forsufficiently large N , Re( G ( λ )) < G ( λ )) is dominated by N . For sufficiently large N and for all a ∈ [3 N/ , N/ G ( λ )) > G ( λ )) is dominated by N . So G (Γ ) starts in the first quadrantand ends in the second quadrant.Consider G ( λ ) along Γ , where b = N and a decreases from N/ m . At the end of Γ and forsufficiently large N , Re( G ( λ )) < G ( λ )) is dominated by N . Im( G ( λ )) could be positive ornegative. So G (Γ ) starts in the second quadrant and ends in the left half-plane.Consider G ( λ ) along Γ , where a = m and b decreases from N to − N , which completes one circuitaround Γ in a counter-clockwise direction. If we fix N to be as large as required then we can choose m the Equations (36) and (37) are approximated arbitrarily closely bythe equations Re( G ( b )) = − b − b (2 α S − r Z U S sin( br F z ) g S ( z ) d z − b (2 α F − Z U F sin( br F z ) g F ( z ) d z + (2 α S + 2 α F − r Z U S Z U F cos( br F ( z + v )) g F ( v ) g S ( z ) d v d z , (38)Im( G ( b )) = − b (2 α S − r Z U S cos( br F z ) g S ( z ) d z − b (2 α F − Z U F cos( br F z ) g F ( z ) d z − (2 α S + 2 α F − r Z U S Z U F sin( br F ( z + v )) g F ( v ) g S ( z ) d v d z . (39)For notational simplicity, define the functions f , f , f , g , g , and g by f ( b ) = − b (2 α S − r Z U S sin( br F z ) g S ( z ) d z , (40) f ( b ) = − b (2 α F − Z U F sin( br F z ) g F ( z ) d z , (41) f ( b ) = (2 α S + 2 α F − r Z U S Z U F cos( br F ( z + v )) g F ( v ) g S ( z ) d v d z , (42) g ( b ) = − b (2 α S − r Z U S cos( br F z ) g S ( z ) d z , (43) g ( b ) = − b (2 α F − Z U F cos( br F z ) g F ( z ) d z , (44) g ( b ) = − (2 α S + 2 α F − r Z U S Z U F sin( br F ( z + v )) g F ( v ) g S ( z ) d v d z , (45)so that Equations (38) and (39) becomeRe( G ( b )) = − b + f ( b ) + f ( b ) + f ( b ) , (46)Im( G ( b )) = g ( b ) + g ( b ) + g ( b ) . (47)Now, consider decreasing b from N to − N . The curves f ( b ) + ig ( b ) and f ( b ) + ig ( b ) have the sameorientation as the spiral − b sin( b ) − ib cos( b ), which is traversed counter-clockwise as b decreases. Similarly,the curve f ( b ) + ig ( b ) has the same orientation as the circle cos( b ) − i sin( b ), which is also counter-clockwise. Note that for discrete delays the curves f ( b ) + ig ( b ) and f ( b ) + ig ( b ) are spirals and the20urve f ( b ) + ig ( b ) is a circle. The sum of these three curves, P k =1 f k ( b ) + ig k ( b ), has counter-clockwiseorientation, and the − b term in Re( G ( b )) translates these curves along the negative real axis. It followsthat if G ( b ) encircles the origin as b decreases from N to − N then it does so in a counter-clockwisedirection. In particular, G ( b ) does not encircle the origin in a clockwise direction. So, for sufficientlysmall m , G (Γ ) can only encircle the origin in a counter-clockwise direction.To ensure that the image contour G (Γ) encircles the origin at least once, note that G (Γ ) crosses thepositive real axis in a counter-clockwise direction at approximately the point (2 α S + 2 α F − r > m , so it follows that G (Γ ) must cross the negative real axis at a point closer to the startof Γ . Therefore, the symmetry of Im( G ( λ )) with respect to the real axis implies that G (Γ) completesa cycle around the origin before the end of G (Γ ). We conclude that Ind G (Γ) (0) ≥
1, and our proof iscomplete.
Proposition 8 (Equilibrium point (1 ,
0) when β S − β F = 0) . Consider the system of delay differential equations (2) and (3) with parameters, delay kernels, and historyfunctions which satisfy (4) – (10) . For all β S , β F ∈ R , (1 , is locally stable when β S − β F < and locallyunstable when β S − β F > .Proof. For ( S * , F * ) = (1 , x ( t )d t = (2 α S − r (cid:0) − x ( t ) − y ( t ) (cid:1) − βy ( t ) , (48)d y ( t )d t = 2(1 − α S ) r (cid:0) − x ( t ) − y ( t ) (cid:1) + βy ( t ) . (49)Equations (48) and (49) have a solution of the form in Equation (23), so substitution gives e λt λ λ c c = − (2 α S − r − (2 α S − r − β − − α S ) r − − α S ) r + β c c e λt . (50)To ensure that ( c , c ) | = 0 we must have the characteristic equation G ( λ ) = (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) − (2 α S − r − λ − (2 α S − r − β − − α S ) r − − α S ) r + β − λ (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) = λ + λ ( r − β ) − βr (51)equal to zero. The zeros of G ( λ ) are the eigenvalues, given by λ = −
12 ( r − β ) ± | r + β | = − r , β, (52)21nd the result follows. Proposition 9 (Equilibrium point (0 ,
1) when β S − β F = 0) . Consider the system of delay differential equations (2) and (3) with parameters, delay kernels, and historyfunctions which satisfy (4) – (10) . For all β S , β F ∈ R , (0 , is locally stable when β S − β F > and locallyunstable when β S − β F < .Proof. For ( S * , F * ) = (0 , x ( t )d t = 2(1 − α F ) (cid:0) − x ( t ) − y ( t ) (cid:1) − βx ( t ) , (53)d y ( t )d t = (2 α F − (cid:0) − x ( t ) − y ( t ) (cid:1) + βx ( t ) . (54)Equations (53) and (54) have a solution of the form in Equation (23), so substitution gives e λt λ λ c c = − − α F ) − β − − α F ) − (2 α F −
1) + β − (2 α F − c c e λt . (55)To ensure that ( c , c ) | = 0 we must have the characteristic equation G ( λ ) = (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) − − α F ) − β − λ − − α F ) − (2 α F −
1) + β − (2 α F − − λ (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) = λ + λ (1 + β ) + β (56)equal to zero. The zeros of G ( λ ) are the eigenvalues, given by λ = −
12 (1 + β ) ± | − β | = − − β, (57)and the result follows. Proposition 10 (Equilibrium points ( u, − u ) for u ∈ [0 ,
1] when β S = β F ) . Consider the system of delay differential equations (2) and (3) with parameters, delay kernels, and historyfunctions which satisfy (4) – (10) . For all β S , β F ∈ R such that β S = β F and for all u ∈ [0 , the equilibriumpoint ( u, − u ) is locally stable.Proof. For ( S * , F * ) = ( u, − u ) with u ∈ R and β = 0, Equations (19) and (20) becomed x ( t )d t = (cid:0) (2 α S − ru + 2(1 − α F )(1 − u ) (cid:1)(cid:0) − x ( t ) − y ( t ) (cid:1) , (58)d y ( t )d t = (cid:0) − α S ) ru + (2 α F − − u ) (cid:1)(cid:0) − x ( t ) − y ( t ) (cid:1) . (59)22quations (58) and (59) have a solution of the form in Equation (23), so substitution gives e λt λ λ c c = − (2 α S − ru − (2 α S − ru − − α F )(1 − u ) − − α F )(1 − u ) − − α S ) ru − − α S ) ru − (2 α F − − u ) − (2 α F − − u ) c c e λt . (60)To ensure that ( c , c ) | = 0 we must have the characteristic equation G ( λ ) = (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) − (2 α S − ru − λ − (2 α S − ru − − α F )(1 − u ) − − α F )(1 − u ) − − α S ) ru − − α S ) ru − λ − (2 α F − − u ) − (2 α F − − u ) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) = λ − λ (cid:0) u (1 − r ) − (cid:1) (61)equal to zero. The zeros of G ( λ ) are the eigenvalues, given by λ = 0, u (1 − r ) −
1. (62)Since r ∈ (0 , u ∈ [0 ,
1] we have u (1 − r ) − <
0, therefore ( u, − u ) is locally stable. We obtain numerical solutions of Equations (2) and (3) using the forward Euler method, for whichthe temporal domain, [0 , t = 0 . S ( t ) and F ( t ), which are obtained by interpolating between the previously estimated values for S ( t ) and F ( t ). The interpolation is achieved using piecewise cubic Hermite interpolating polynomials, which areshape preserving. The sizes of the time step and the integration subintervals ensure grid-independencefor our results. Examples of the simulations are shown in Figure 3.The delay kernels in our model are set as probability density functions of a right-truncated Erlangdistribution (Electronic Supplementary Material), shown in Figure 3(a). For slow-proliferating cells theErlang density has shape k = 12000 and rate λ = 20 h − with mean 600 h, and is truncated at U S = 700h. For fast-proliferating cells the Erlang density has shape k = 20 and rate λ = 1 h − with mean 20 h,23 a) (b) S ( t ) , F ( t ) S ( t ) F ( t ) t (h)α S = 1, α F = 1, β S = 0 h -1 , β F = 0 h -1 (c) S ( t ) , F ( t ) S ( t ) F ( t ) t (h)α S = 1, α F = 0.7, β S = 0 h -1 , β F = 0 h -1 (d) S ( t ) , F ( t ) S ( t ) F ( t ) t (h)α S = 1, α F = 0.5, β S = 0 h -1 , β F = 0 h -1 (e) S ( t ) , F ( t ) S ( t ) F ( t ) t (h)α S = 1, α F = 1, β S = 0 h -1 , β F = 0.008 h -1 (f) S ( t ) , F ( t ) S ( t ) F ( t ) t (h)α S = 1, α F = 0.7, β S = 0 h -1 , β F = 0.008 h -1 (g) S ( t ) , F ( t ) S ( t ) F ( t ) t (h)α S = 1, α F = 1, β S = 0.008 h -1 , β F = 0 h -1 (h) S ( t ) , F ( t ) S ( t ) F ( t ) t (h)α S = 1, α F = 0.5, β S = 0.008 h -1 , β F = 0 h -1 (i) S ( t ) , F ( t ) t (h)α S = 1, α F = 1, β S = 0 h -1 , β F = 0.024 h -1 S ( t ) F ( t ) (j) S ( t ) , F ( t ) t (h) S ( t ) F ( t ) α S = 1, α F = 1, β S = 0.016 h -1 , β F = 0 h -1 Slow-proliferating cellsFast-proliferating cells
100 700
Figure 3:
Numerical simulations of our model in Equations (2) and (3). (a) Each delay kernel is the probabilitydensity function (PDF) of a right-truncated Erlang distribution (Electronic Supplementary Material). For slow-proliferating cells the Erlang density function has shape k = 12000 and rate λ = 20 h − with mean 600 h, and istruncated at U S = 700 h. For fast-proliferating cells the Erlang density function has shape k = 20 and rate λ = 1h − with mean 20 h, and is truncated at U F = 100 h. (b)–(j) The simulations all use the parameters K = 500and r S = r F = 0 . − . For (b)–(h) the history functions are φ S ( t ) = 100 e r S t and φ F ( t ) = 100 e r F t , for (i) thehistory functions are φ S ( t ) = 10 − e r S t and φ F ( t ) = 100 e r F t , and for (j) the history functions are φ S ( t ) = 100 e r S t and φ F ( t ) = 10 − e r F t . Parameters specific to each simulation, namely α S , α F , β S , and β F , are indicated on thefigure U F = 100 h. All simulations use the parameters K = 500 and r S = r F = 0 . − . Theparameters α S , α F , β S , and β F are varied for the different simulations, as indicated in Figure 3(b)–(h).There are many options for the functional form of the history functions. One simple option is touse constant functions, however it is reasonable to assume that the cells grew exponentially in the past,so we use exponential functions with growth rates equal to the intrinsic growth rates of the slow- andfast-proliferating cells. The history functions are φ S ( t ) = 100 e r S t for t ∈ [ − ,
0] and φ F ( t ) = 100 e r F t for t ∈ [ − , φ S ( t ) = 10 − e r S t for t ∈ [ − ,
0] and φ F ( t ) = 100 e r F t for t ∈ [ − , φ S ( t ) = 100 e r S t for t ∈ [ − ,
0] and φ F ( t ) = 10 − e r F t for t ∈ [ − , C in Equation (11), choosing a different historyfunction in C results in a different solution. When β S = β F , different history functions may producedifferent transient dynamics, whereas the solutions have the same long-term behaviour. When β S = β F ,however, there are infinitely many equilibrium points so different history functions can produce solutionswith different long-term behaviour.Figure 3(b)–(d) shows simulations with β S = β F = 0 h − , so no induced switching between slow andfast proliferation. By Theorem 5 there are infinitely many locally-stable equilibrium points correspondingto the line segment between ( K,
0) and (0 , K ). The different equilibrium states are obtained by varyingthe levels of asymmetric division through α S and α F , or using different history functions.In Figure 3(e)–(f) we show simulations with β S = 0 h − and β F = 0 .
008 h − , therefore inducedswitching only from fast to slow proliferation. By Theorem 5 the equilibrium point ( K,
0) is locallystable and the equilibrium point (0 , K ) is locally unstable. In Figure 3(g)–(h) we show simulations with β S = 0 .
008 h − and β F = 0 h − , so induced switching only from slow to fast proliferation. By Theorem 5the equilibrium point (0 , K ) is locally stable and the equilibrium point ( K,
0) is locally unstable. Thesesimulations illustrate that induced switching determines the long-term behaviour of the solutions, whileasymmetric division only influences the transient dynamics.In Figure 3(i)–(j) we set one of the history functions φ S or φ F close to zero over its domain, illustratinghow a very small subpopulation can become the main subpopulation through induced switching, possiblyrequiring a long time period. It is particularly interesting that, in Figure 3(i), the density of the fast-proliferating cells is very close to carrying capacity and appears to be at equilibrium for a long time,however through induced switching the slow-proliferating cells eventually become the main subpopulation.25 Discussion and outlook
Proliferative heterogeneity in cancer cell populations constitutes a crucial challenge for cancer therapy,as slow-proliferating cells tend to be highly aggressive, have increased resistance to cytotoxic drugs, andcan replenish the fast-proliferating subpopulation [Ahn et al., 2017,Ahmed and Haass, 2018,Perego et al.,2018, Vallette et al., 2019]. The dynamics underlying tumour heterogeneity are not well understood, soimproving cancer therapy depends on furthering this understanding [Haass et al., 2014,Haass, 2015,Ahmedand Haass, 2018]. Theoretical approaches are well-placed to assist in elucidating the transient dynamicsof intratumoural heterogeneity.In this article we present a delay differential equation model for heterogeneous cell proliferation inwhich the total population consists of a subpopulation of slow-proliferating cells and a subpopulationof fast-proliferating cells. Our model incorporates the two cellular processes of asymmetric cell divisionand induced switching between proliferative states, which are important contributors to the dynamicheterogeneity of a cancer cell population [Nelson and Chen, 2002,Bajaj et al., 2015,Dey-Guha et al., 2015,West and Newton, 2019]. The model is designed for investigating the transient dynamics of intratumouralheterogeneity with respect to cell proliferation. We employ delay differential equations in our model ratherthan ordinary differential equations in order to obtain transient dynamics consistent with the dynamicsin a tumour. While the equilibrium states for our model are the same as those for the correspondingordinary differential equations, the transient dynamics are very different, and model parameterisationwith biologically-realistic values requires a model that incorporates realistic cell cycle durations for theslow- and fast-proliferating subpopulations.Because the transient dynamics of a tumour are of primary interest, and the local stability analysis ofour model provides only long-term behaviour near to the equilibrium, we must numerically simulate ourmodel to explore the transient dynamics. We provide some examples of numerical simulations in Figure 3,where we specify the delay kernels to be right-truncated Erlang distributions (Section 2 of ElectronicSupplementary Material), and vary the parameters to demonstrate some of the possible dynamics within atumour cell population. To exemplify some of the experimental scenarios to which our model is applicable,we consider a tumour that is treated with a cytotoxic drug which may induce cellular stress, causing thefast-proliferating cells to switch to the drug-resistant slow-proliferating phenotype [Ahmed and Haass,2018] through the mechanisms of asymmetric cell division and induced switching.We show simulations where there is no induced switching between the slow- and fast-proliferatingsubpopulations in Figure 3(b)–(d). If tumour cells experience no microenvironmental stress, then all celldivisions may be symmetric, corresponding to the situation in Figure 3(b) where the fast-proliferating cells26apidly populate the tumour until the total cell density reaches carrying capacity. If a drug is introducedinto the tumour microenvironment then the fast-proliferating cells may experience cellular stress, inducingthe fast-proliferating cells into asymmetric cell division as a survival strategy, as the slow-proliferatingphenotype is drug resistant. Figure 3(c)–(d) illustrates this behaviour, first for an intermediate level ofasymmetric division in (c), and then for the maximum level of asymmetric division in (d).Figure 3(e) shows a simulation where fast-proliferating cells are under stress due to the presence of adrug, and are induced to switch to slow proliferation through signals from slow-proliferating cells as a sur-vival strategy. Alternatively, the induced signalling may arise from the highly invasive slow-proliferatingcells [Chapman et al., 2014, Ahmed and Haass, 2018] influencing the less invasive fast-proliferating cellsto switch to the more invasive slow-proliferating phenotype. In our simulation the total cell populationappears to reach the carrying capacity density at around 200 hours, however the dynamics of inducedswitching of cells from fast to slow proliferation continues until all cells are slow proliferating. This isimportant because, while tumour growth has effectively ceased, the tumour is becoming increasingly drugresistant and invasive over time until the whole tumour is composed of drug resistant and invasive cells.Therefore, effective early treatment of the tumour is required in order to prevent the tumour from becom-ing more aggressive and treatment resistant. If the per capita interaction strength of slow-proliferatingcells to induce fast-proliferating cells to switch to slow proliferation is obtained experimentally, then ourmodel could be used to predict the transient change in the proportion of slow-proliferating cells in thepopulation, and therefore the changes in invasiveness and drug resistance of the tumour.Now consider a tumour composed of mostly fast-proliferating cells and a very small proportion ofslow-proliferating cells, as in Figure 3(i). The fast-proliferating cells undergo induced switching to slowproliferation, perhaps due to stress from an introduced drug or to increase invasiveness. Initially thefast-proliferating subpopulation grows to reach near the carry capacity density, and the system appearsto be in equilibrium for an extended period of time. As the tumour is almost completely composed offast-proliferating cells, it is in the least invasive and most drug sensitive state. Without knowledge of thepresence of induced switching, experimental investigations may not reveal that the tumour is in the processof becoming highly invasive and drug resistant. Indeed, once the density of the slow-proliferating cells hasreached a sufficient but still very low level, the tumour rapidly becomes populated by slow-proliferatingcells through induced switching of the fast-proliferating cells.Finally, consider a small tumour comprised mostly of slow-proliferating cells and a very small pro-portion of fast-proliferating cells, as in Figure 3(j). Signals from the fast-proliferating cells induce theslow-proliferating cells to switch to fast proliferation. Experimentally, this could correspond to a tumourthat has been treated with a drug which caused the death of most of the fast-proliferating cells. For27n extended period of time the tumour grows very little, until the density of fast-proliferating cells ishigh enough that the induced switching from slow- to fast-proliferation rapidly grows the tumour to themaximum sustainable size. Our model is therefore able to provide an estimate of tumour growth overtime following drug treatment, when the cells can undergo induced switching.There are numerous possibilities for future work. Induced switching between proliferative statescould take many forms. In tumours the slow-proliferating state may continually arise and disappear[Roesch et al., 2010], so it would be interesting to accommodate time-dependent induced switching intothe model, which could be either periodic or aperiodic. We could also consider the induced switchingto have an explicit dependence on density, so that no switching occurs from a particular proliferativestate when the density of cells from the other proliferative state is above a certain value. A similarexplicit density dependence could be implemented for asymmetric cell division, which occurs at constantproportions in our current model. These explicit dependences on density may be relevant for slow-proliferating subpopulations in tumours that appear to persist over time and maintain the relative size ofthe subpopulation [Perego et al., 2018, Vallette et al., 2019]. Our model could also be extended to includethe additional process of spontaneous switching between proliferative states, which is independent of othercells and may be stochastic.While our model has implicit spatial structure, since the dependent variables are cell densities, wecould include spatial structure explicitly. We could then explicitly model cell migration with a diffusiveterm [Vittadello et al., 2018]. Further, induced switching could be modelled as a more localised processwhere the rate of a cell switching proliferative states is determined by the density of cells in a differentproliferative state within a given radius of the cell, where the interaction strength decreases with distancefrom the cell.We could also extend our model to more than two dependent variables. For example, we couldconsider fast-, slow-, and very-slow-proliferating subpopulations. Another possible extension is to includeapoptosis. Much of our analysis in this article is likely to be easily generalised to an extended versionof our model. The more challenging aspect could be the analysis of the corresponding transcendentalcharacteristic equations, however taking a more abstract approach for an extended model with an arbitrary n dependent variables could simplify the problem. 28 ode availability The code for the algorithm to replicate the numerical simulations in this work is available on GitHub athttps://github.com/DrSeanTVittadello/Vittadello2020.
Acknowledgements
NKH is a Cameron fellow of the Melanoma and Skin Cancer Research Institute, and is supported by theNational Health and Medical Research Council of Australia (APP1084893). MJS is supported by theAustralian Research Council (DP170100474).
Conflict of interest
The authors declare that they have no conflict of interest.
References [Ahlfors, 1979] Ahlfors, L. V. (1979).
Complex Analysis . McGraw-Hill, third edition.[Ahmed and Haass, 2018] Ahmed, F. and Haass, N. K. (2018). Microenvironment-driven dynamic het-erogeneity and phenotypic plasticity as a mechanism of melanoma therapy resistance.
Frontiers inOncology , 8:173.[Ahn et al., 2017] Ahn, A., Chatterjee, A., and Eccles, M. R. (2017). The slow cycling phenotype: Agrowing problem for treatment resistance in melanoma.
Molecular Cancer Therapeutics , 16:1002–1009.[Arino, 1995] Arino, O. (1995). A survey of structured cell population dynamics.
Acta Biotheoretica ,43:3–25.[Arino and Kimmel, 1989] Arino, O. and Kimmel, M. (1989). Asymptotic behavior of a nonlinearfunctional-integral equation of cell kinetics with unequal division.
Journal of Mathematical Biology ,27:341–354.[Bajaj et al., 2015] Bajaj, J., Zimdahl, B., and Reya, T. (2015). Fearful symmetry: subversion of asym-metric division in cancer development and progression.
Cancer Research , 75:792–797.[Baker et al., 1997] Baker, C. T. H., Bocharov, G. A., and Paul, C. A. H. (1997). Mathematical modellingof the interleukin-2 T-cell system: a comparative study of approaches based on ordinary and delaydifferential equation.
Journal of Theoretical Medicine , 1:117–128.29Baker et al., 1998] Baker, C. T. H., Bocharov, G. A., Paul, C. A. H., and Rihan, F. A. (1998). Modellingand analysis of time-lags in some basic patterns of cell proliferation.
Journal of Mathematical Biology ,37:341–371.[Beaumont et al., 2016] Beaumont, K. A., Hill, D. S., Daignault, S. M., Lui, G. Y., Sharp, D. M.,Gabrielli, B., Weninger, W., and Haass, N. K. (2016). Cell cycle phase-specific drug resistance asan escape mechanism of melanoma cells.
Journal of Investigative Dermatology , 136:1479–1489.[Billy et al., 2014] Billy, F., Clairambaultt, J., Fercoq, O., Gaubertt, S., Lepoutre, T., Ouillon, T., andSaito, S. (2014). Synchronisation and control of proliferation in cycling cell population models withage structure.
Mathematics and Computers in Simulation , 96:66–94.[Byrne and Drasdo, 2009] Byrne, H. and Drasdo, D. (2009). Individual-based and continuum models ofgrowing cell populations: a comparison.
Journal of Mathematical Biology , 58:657–687.[Byrne, 1997] Byrne, H. M. (1997). The effect of time delays on the dynamics of avascular tumor growth.
Mathematical Biosciences , 144:83–117.[Cai et al., 2007] Cai, A. Q., Landman, K. A., and Hughes, B. D. (2007). Multi-scale modeling of awound-healing cell migration assay.
Journal of Theoretical Biology , 245:576–594.[Cassidy et al., 2019] Cassidy, T., Craig, M., and Humphries, A. R. (2019). Equivalences between agestructured models and state dependent distributed delay differential equations.
Mathematical Bio-sciences and Engineering , 16:5419–5450.[Cassidy and Humphries, 2020] Cassidy, T. and Humphries, A. R. (2020). A mathematical model of viraloncology as an immuno-oncology instigator.
Mathematical Medicine and Biology , 37:117–151.[Chao et al., 2019] Chao, H. X., Fakhreddin, R. I., Shimerov, H. K., Kedziora, K. M., Kumar, R. J.,Perez, J., Limas, J. C., Grant, G. D., Cook, J. G., Gupta, G. P., and Purvis, J. E. (2019). Evidencethat the human cell cycle is a series of uncoupled, memoryless phases.
Molecular Systems Biology ,15:e8604.[Chapman et al., 2014] Chapman, A., Fernandez del Ama, L., Ferguson, J., Kamarashev, J., Wellbrock,C., and Hurlstone, A. (2014). Heterogeneous tumor subpopulations cooperate to drive invasion.
CellReports , 8:688–695.[Clairambault and Fercoq, 2016] Clairambault, J. and Fercoq, O. (2016). Physiologically structured cellpopulation dynamic models with applications to combined drug delivery optimisation in oncology.
Mathematical Modelling of Natural Phenomena , 11:45–70.30Dey-Guha et al., 2015] Dey-Guha, I., Alves, C. P., Yeh, A. C., Salony, Sole, X., Darp, R., and Ra-maswamy, S. (2015). A mechanism for asymmetric cell division resulting in proliferative asynchronicity.
Molecular Cancer Research , 13:223–230.[Dey-Guha et al., 2011] Dey-Guha, I., Wolfer, A., Yeh, A. C., Albeck, J. G., Darp, R., Leon, E.,Wulfkuhle, J., Petricoin, E. F., Wittner, B. S., and Ramaswamy, S. (2011). Asymmetric cancer celldivision regulated by AKT.
Proceedings of the National Academy of Sciences of the United States ofAmerica , 108:12845–12850.[Diekmann et al., 1995] Diekmann, O., van Gils, S. A., Verduyn Lunel, S. M., and Walther, H.-O. (1995).
Delay Equations . Springer New York.[Engelborghs et al., 2000] Engelborghs, K., Luzyanina, T., and Roose, D. (2000). Numerical bifurcationanalysis of delay differential equations.
Journal of Computational and Applied Mathematics , 125:265–275.[Gabriel et al., 2012] Gabriel, P., Garbett, S. P., Quaranta, V., Tyson, D. R., and Webb, G. F. (2012).The contribution of age structure to cell population responses to targeted therapeutics.
Journal ofTheoretical Biology , 311:19–27.[Gallaher et al., 2019] Gallaher, J. A., Brown, J. S., and Anderson, A. R. A. (2019). The impact ofproliferation-migration tradeoffs on phenotypic evolution in cancer.
Scientific Reports , 9:2425.[Gavagnin et al., 2019] Gavagnin, E., Ford, M. J., Mort, R. L., Rogers, T., and Yates, C. A. (2019). Theinvasion speed of cell migration models with realistic cell cycle time distributions.
Journal of TheoreticalBiology , 481:91–99.[Getto et al., 2019] Getto, P., Gyllenberg, M., Nakata, Y., and Scarabel, F. (2019). Stability analysis ofa state-dependent delay differential equation for cell maturation: analytical and numerical methods.
Journal of Mathematical Biology , 79:281–328.[Getto and Waurick, 2016] Getto, P. and Waurick, M. (2016). A differential equation with state-dependent delay from cell population biology.
Journal of Differential Equations , 260:6176–6200.[Greene et al., 2015] Greene, J. M., Levy, D., Fung, K. L., Souza, P. S., Gottesman, M. M., and Lavi,O. (2015). Modeling intrinsic heterogeneity and growth of cancer cells.
Journal of Theoretical Biology ,367:262–277.[Haass, 2015] Haass, N. K. (2015). Dynamic tumor heterogeneity in melanoma therapy: how do weaddress this in a novel model system?
Melanoma Management , 2:93–95.31Haass et al., 2014] Haass, N. K., Beaumont, K. A., Hill, D. S., Anfosso, A., Mrass, P., Munoz, M. A.,Kinjyo, I., and Weninger, W. (2014). Real-time cell cycle imaging during melanoma growth, invasion,and drug response.
Pigment Cell & Melanoma Research , 27:764–776.[Hanahan and Weinberg, 2011] Hanahan, D. and Weinberg, R. A. (2011). Hallmarks of cancer: The nextgeneration.
Cell , 144:646–674.[Huang et al., 2016] Huang, C., Cao, J., Wen, F., and Yang, X. (2016). Stability analysis of SIR modelwith distributed delay on complex networks.
PLOS ONE , 11:e0158813.[Jin et al., 2018] Jin, W., McCue, S. W., and Simpson, M. J. (2018). Extended logistic growth model forheterogeneous populations.
Journal of Theoretical Biology , 445:51–61.[Kaslik and Neamtu, 2018] Kaslik, E. and Neamtu, M. (2018). Stability and hopf bifurcation analysisfor the hypothalamic-pituitary-adrenal axis model with memory.
Mathematical Medicine and Biology ,35:49–78.[Khasawneh and Mann, 2011] Khasawneh, F. A. and Mann, B. P. (2011). Stability of delay integro-differential equations using a spectral element method.
Mathematical and Computer Modelling , 54:2493–2503.[Kuang, 1993] Kuang, Y. (1993).
Delay Differential Equations: With Applications in Population Dynam-ics . Academic Press.[Lebowitz and Rubinow, 1974] Lebowitz, J. L. and Rubinow, S. I. (1974). A theory for the age andgeneration time distribution of a microbial population.
Journal of Mathematical Biology , 1:17–36.[Lu, 1991] Lu, L. (1991). Numerical stability of the θ -methods for systems of differential equations withseveral delay terms. Journal of Computational and Applied Mathematics , 34:291–304.[Mackey and Rudnicki, 1994] Mackey, M. C. and Rudnicki, R. (1994). Global stability in a delayed partialdifferential equation describing cellular replication.
Journal of Mathematical Biology , 33:89–109.[Maini et al., 2004] Maini, P. K., McElwain, D. L. S., and Leavesley, D. I. (2004). Traveling wave modelto interpret a wound-healing cell migration assay for human peritoneal mesothelial cells.
Tissue Engi-neering , 10:475–482.[Matson and Cook, 2017] Matson, J. P. and Cook, J. G. (2017). Cell cycle proliferation decisions: theimpact of single cell analyses.
The FEBS Journal , 284:362–375.32McClatchey and Yap, 2012] McClatchey, A. I. and Yap, A. S. (2012). Contact inhibition (of proliferation)redux.
Current Opinion in Cell Biology , 24:685–694.[McCluskey, 2010] McCluskey, C. C. (2010). Global stability of an sir epidemic model with delay andgeneral nonlinear incidence.
Mathematical Biosciences and Engineering , 7:837–850.[Moore and Lyle, 2011] Moore, N. and Lyle, S. (2011). Quiescent, slow-cycling stem cell populations incancer: A review of the evidence and discussion of significance.
Journal of Oncology , 2011:396076.[Nelson and Chen, 2002] Nelson, C. M. and Chen, C. S. (2002). Cell–cell signaling by direct contactincreases cell proliferation via a PI3K-dependent signal.
FEBS Letters , 514:238–242.[Pavel et al., 2018] Pavel, M., Renna, M., Park, S. J., Menzies, F. M., Ricketts, T., Fllgrabe, J., Ashke-nazi, A., Frake, R. A., Lombarte, A. C., Bento, C. F., Franze, K., and Rubinsztein, D. C. (2018).Contact inhibition controls cell survival and proliferation via YAP/TAZ-autophagy axis.
Nature Com-munications , 9:2961.[Perego et al., 2018] Perego, M., Maurer, M., Wang, J. X., Shaffer, S., M¨uller, A. C., Parapatics, K., Li,L., Hristova, D., Shin, S., Keeney, F., Liu, S., Xu, X., Raj, A., Jensen, J. K., Bennett, K. L., Wagner,S. N., Somasundaram, R., and Herlyn, M. (2018). A slow-cycling subpopulation of melanoma cells withhighly invasive properties.
Oncogene , 37:302–312.[Puliafito et al., 2012] Puliafito, A., Hufnagel, L., Neveu, P., Streichan, S., Sigal, A., Fygenson, D. K., andShraiman, B. I. (2012). Collective and single cell behavior in epithelial contact inhibition.
Proceedingsof the National Academy of Sciences of the United States of America , 109:739–744.[Roesch et al., 2010] Roesch, A., Fukunaga-Kalabis, M., Schmidt, E. C., Zabierowski, S. E., Brafford,P. A., Vultur, A., Basu, D., Gimotty, P., Vogt, T., and Herlyn, M. (2010). A temporarily distinctsubpopulation of slow-cycling melanoma cells is required for continuous tumor growth.
Cell , 141:583–594.[Rudin, 1986] Rudin, W. (1986).
Real and Complex Analysis . McGraw-Hill Education, third edition.[Sakaue-Sawano et al., 2008] Sakaue-Sawano, A., Kurokawa, H., Morimura, T., Hanyu, A., Hama, H.,Osawa, H., Kashiwagi, S., Fukami, K., Miyata, T., Miyoshi, H., Imamura, T., Ogawa, M., Masai, H.,and Miyawaki, A. (2008). Visualizing spatiotemporal dynamics of multicellular cell-cycle progression.
Cell , 132:487–498.[Sarapata and de Pillis, 2014] Sarapata, E. A. and de Pillis, L. G. (2014). A comparison and catalog ofintrinsic tumor growth models.
Bulletin of Mathematical Biology , 76:2010–2024.33Scott et al., 2013] Scott, J. G., Basanta, D., Anderson, A. R. A., and Gerlee, P. (2013). A mathematicalmodel of tumour self-seeding reveals secondary metastatic deposits as drivers of primary tumour growth.
Journal of the Royal Society Interface , 10:20130011.[Sherratt and Murray, 1990] Sherratt, J. A. and Murray, J. D. (1990). Models of epidermal wound healing.
Proceedings of the Royal Society B , 241:29–36.[Simpson et al., 2018] Simpson, M. J., Jin, W., Vittadello, S. T., Tambyah, T. A., Ryan, J. M., Gunas-ingh, G., Haass, N. K., and McCue, S. W. (2018). Stochastic models of cell invasion with fluorescentcell cycle indicators.
Physica A: Statistical Mechanics and its Applications , 510:375–386.[Smalley and Herlyn, 2009] Smalley, K. S. M. and Herlyn, M. (2009). Integrating tumor-initiating cellsinto the paradigm for melanoma targeted therapy.
International Journal of Cancer , 124:1245–1250.[Smith, 2011] Smith, H. (2011).
An Introduction to Delay Differential Equations with Applications to theLife Sciences . Springer-Verlag GmbH.[Spoerri et al., 2017] Spoerri, L., Beaumont, K. A., Anfosso, A., and Haass, N. K. (2017). Real-time cellcycle imaging in a 3d cell culture model of melanoma.
Methods in Molecular Biology , 1612:401–416.[Sun, 2006] Sun, L. (2006). Stability analysis for delay differential equations with multidelays and nu-merical examples.
Mathematics of Computation , 75:151–165.[Swanson et al., 2003] Swanson, K. R., Bridge, C., Murray, J., and Alvord, E. C. (2003). Virtual andreal brain tumors: using mathematical modeling to quantify glioma growth and invasion.
Journal ofthe Neurological Sciences , 216:1–10.[Vallette et al., 2019] Vallette, F. M., Olivier, C., L´ezot, F., Oliver, L., Cochonneau, D., Lalier, L.,Cartron, P.-F., and Heymann, D. (2019). Dormant, quiescent, tolerant and persister cells: Four syn-onyms for the same target in cancer.
Biochemical Pharmacology , 162:169–176.[Vermeulen et al., 2003] Vermeulen, K., Van Bockstaele, D. R., and Berneman, Z. N. (2003). The cellcycle: a review of regulation, deregulation and therapeutic targets in cancer.
Cell Proliferation , 36:131–149.[Villasana and Radunskaya, 2003] Villasana, M. and Radunskaya, A. (2003). A delay differential equationmodel for tumor growth.
Journal of Mathematical Biology , 47:270–294.34Vittadello et al., 2018] Vittadello, S. T., McCue, S. W., Gunasingh, G., Haass, N. K., and Simpson,M. J. (2018). Mathematical models for cell migration with real-time cell cycle dynamics.
BiophysicalJournal , 114:1241–1253.[Vittadello et al., 2019] Vittadello, S. T., McCue, S. W., Gunasingh, G., Haass, N. K., and Simpson,M. J. (2019). Mathematical models incorporating a multi-stage cell cycle replicate normally-hiddeninherent synchronization in cell proliferation.
Journal of the Royal Society Interface , 16:20190382.[Vittadello et al., 2020] Vittadello, S. T., McCue, S. W., Gunasingh, G., Haass, N. K., and Simpson,M. J. (2020). Examining go-or-grow using fluorescent cell-cycle indicators and cell-cycle-inhibitingdrugs.
Biophysical Journal , 118:1243–1247.[Webb, 1986] Webb, G. F. (1986). A model of proliferating cell populations with inherited cycle length.
Journal of Mathematical Biology , 23:269–282.[Weber et al., 2014] Weber, T. S., Jaehnert, I., Schichor, C., Or-Guil, M., and Carneiro, J. (2014). Quan-tifying the length and variance of the eukaryotic cell cycle phases by a stochastic model and dualnucleoside pulse labelling.
PLoS Computational Biology , 10:e1003616.[West and Newton, 2019] West, J. and Newton, P. K. (2019). Cellular interactions constrain tumorgrowth.
Proceedings of the National Academy of Sciences of the United States of America , 116:1918–1923.[Yates et al., 2017] Yates, C. A., Ford, M. J., and Mort, R. L. (2017). A multi-stage representation ofcell proliferation as a markov process.
Bulletin of Mathematical Biology , 79:2905–2928.[Zhu and Thompson, 2019] Zhu, J. and Thompson, C. B. (2019). Metabolic regulation of cell growth andproliferation.
Nature Reviews. Molecular Cell Biology , 20:436–450.35 upplementary Material for:
A novel mathematical model of heterogeneous cell proliferation
S. T. Vittadello, S. W. McCue, G. Gunasingh, N. K. Haass, M. J. Simpson
July 27, 2020 a r X i v : . [ q - b i o . CB ] J u l ontents Experimental
We briefly describe the materials and methods employed to obtain our experimental data for the durationsof cells in G1 and S/G2/M phases, shown in Figure S1. Further experimental details are given in [1].Each histogram is constructed using data from 50 individual cells.Our experimental data consist of microscopy time-series images of two-dimensional cell proliferationassays using the three human melanoma cell lines C8161 (kindly provided by Mary Hendrix, Chicago,IL, USA), WM983C and 1205Lu (both kindly provided by Meenhard Herlyn, Philadelphia, PA, USA),which have cell cycle durations of approximately 21, 23 and 37 h, respectively [2]. The cell lines weregenotypically characterised [3–6], grown as described in [7], and authenticated by STR fingerprinting(QIMR Berghofer Medical Research Institute, Herston, Australia).We maintain the cell cultures to prevent any induced synchronisation from cell cycle arrest in G1phase, by passaging the cells every three days, and on the day prior to setting up an experiment, tomaintain a subconfluent cell density and a fresh growth medium.Experimental investigation of the progression of the cell cycle is visually enabled with fluorescentubiquitination-based cell cycle indicator (FUCCI) technology [8]. FUCCI consists of two genetically-encoded reporters that enable visualisation of the cell cycle of individual live cells: when the cell isin G1 phase the nucleus fluoresces red, and when the cell is in S/G2/M phase the nucleus fluorescesgreen (Figure 1 in the main document). During the transition from G1 to S phase, called early S, bothreporters fluoresce and the nucleus appears yellow. FUCCI is utilised in experimental studies of thecycling dynamics of cells in tumours [2], and reveals the differential cycling of the cell population.2 .350 Duration (h)C8161 in G1 phase P r o b a b ili t y d e n s i t y (a) P r o b a b ili t y d e n s i t y (b) P r o b a b ili t y d e n s i t y (c) P r o b a b ili t y d e n s i t y (d) P r o b a b ili t y d e n s i t y (e) P r o b a b ili t y d e n s i t y (f)
10 20 30 400 10 20 30 40010 20 30 400 10 20 30 40010 20 30 400
S/G2/M-duration dataRight-truncatedErlang distributionRight-truncated Erlang distributionExponential distributionMean of distributions P r o b a b ili t y d e n s i t y (g) G1-duration dataRight-truncatedErlang distribution
Figure S1:
Histograms of G1- and S/G2/M-phase duration data with best fits of the right-truncated Erlangdistribution in Equation (S4), and a comparison of the probability density functions of the right-truncated Erlangand exponential distributions. Each histogram corresponds to 50 cells. The best-fit parameters for the right-truncated Erlang distribution are: (a) RTE(17 , . − ,
20 h); (b) RTE(5 , . − ,
20 h); (c) RTE(7 , .
79 h − ,
30 h);(d) RTE(11 , .
89 h − ,
30 h); (e) RTE(5 , .
30 h − ,
40 h); (f) RTE(8 , .
62 h − ,
40 h). (g) Probability density functionsof the right-truncated Erlang and exponential distributions, for RTE(8 , . − ,
80 h) and Exp(0 .
05 h − ). Theexponential distribution has mean 20 h, and the right-truncated Erlang distribution has mean 20.0000 h (to fourdecimal places) Experimental data and models for cell cycle durations
Standard deterministic mathematical models of cell proliferation, such as exponential (S1) and logistic (S2)growth models, d P ( t )d t = rP ( t ) , (S1)d P ( t )d t = rP ( t ) (cid:18) − P ( t ) K (cid:19) , (S2)where P ( t ) is the population density at time t , r is the intrinsic growth rate, and K is the carryingcapacity density, are based on cell cycle durations with an exponential distribution Exp( θ ) for rate θ .Realising Equations (S1) and (S2) as stochastic pure birth processes yields continuous-time homogeneousMarkov chains with exponentially-distributed durations of cell cycle residence [9]. Experimental investi-gations, however, suggest that the duration of the cell cycle, and in particular each cell cycle phase, isnot exponentially distributed [10–14]. Rather, it is often found that the hypoexponential distribution,characterised as the sum of k independent exponential random variables with distinct rate parameters λ i , for i = 1, . . . , k , is a reasonable distribution for cell cycle duration [10–14]. The hypoexponentialdistribution generalises the Erlang distribution for which the exponential random variables have the samerate parameter λ . Since the exponential distribution allows for a relatively large probability of arbitrarilysmall cell cycle durations, ordinary differential equations such as (S1) and (S2) tend to overestimate thepopulation growth rate, particularly for slow-proliferating cells.We consider only bounded distributions for their greater biological realism. A particular example ofa bounded delay kernel that may be relevant when modelling the cell cycle is the probability density forthe right-truncated Erlang distribution. The probability density for the Erlang distribution is g ( z ) = λ k z k − e − λz ( k − , for z ∈ [0 , ∞ ) , (S3)where k is the shape parameter, λ is the rate parameter, and the mean is k/λ . Restricting the Erlangdistribution to the bounded interval [0 , U ], where U ∈ (0 , ∞ ), gives the right-truncated Erlang distributionwhich has probability density g ∗ ( z ) = g ( z ) Z U g ( w ) d w , for z ∈ [0 , U ] . (S4)We also show in Figure S1 the best fit to the data for the right-truncated Erlang distribution RTE( k, λ, U ),where the shape k and rate λ correspond to the Erlang distribution, and U is the right-truncation point4f the Erlang distribution. Each best fit is obtained with the MATLAB nonlinear least-squares solver lsqnonlin [15] with the trust-region-reflective algorithm [16]. These data demonstrate that the durationsof G1 and S/G2/M, and therefore of the complete cell cycle, are well approximated by a right-truncatedErlang distribution. In Figure S1(g) we compare the probability density functions of the exponentialand right-truncated Erlang distributions, which clearly demonstrates that these cell cycle durations arenot exponentially distributed. Using the lsqnonlin solver we find that for each data set the norm of theresidual from fitting the right-truncated Erlang distribution is no greater than the norm of the residualfrom fitting the Erlang distribution. This outcome is expected as the cumulative density of the Erlangdistribution tends to one exponentially after the maximum data value. We could also left-truncate theErlang distribution to ensure a positive lower bound for the distribution, depending on whether a betterfit would be achieved with particular data. 5 Complete proof of Theorem 2
It suffices to prove that S ( t ) ≥ F ( t ) ≥ t >
0, for then it follows from Lemma 3 that S ( t ) ≤ K and F ( t ) ≤ K for all t >
0. Define t and t by t = inf { t > | S ( t ) < } and t = inf { t > | F ( t ) < } .We consider the infima in the extended real numbers, so that t , t ∈ (0 , + ∞ ], where either infimum isequal to + ∞ if the corresponding set is empty. The proof consists of four separate cases.Case 1: Let β S − β F ≥ S ( t ) < t > t ∈ R . If t < t then choose t ∈ ( t , t ) such that S ( t ) <
0, d S ( t ) / d t | t = t <
0, and the delayssatisfy S ( t ) ≥ F ( t ) ≥
0. Then, since P ( t ) ≤ K by Lemma 3 and since F ( t ) ≥
0, Equation (2)gives d S ( t ) / d t | t = t ≥
0, a contradiction.If t ≥ t then, since F (0) > F ( t ) = 0, there exists t ∈ (0 , t ) such that d F ( t ) / d t | t = t <
0. Then,since P ( t ) ≤ K by Lemma 3, since the delays satisfy S ( t ) ≥ F ( t ) ≥
0, and since S ( t ), F ( t ) ≥ F ( t ) / d t | t = t ≥
0, a contradiction. We conclude that S ( t ) ≥ t > β S − β F < S ( t ) < t > t ∈ R . If t ≤ t then, since S (0) > S ( t ) = 0, there exists t ∈ (0 , t ) such thatd S ( t ) / d t | t = t <
0. Then, since P ( t ) ≤ K by Lemma 3, since the delays satisfy S ( t ) ≥ F ( t ) ≥ S ( t ), F ( t ) ≥
0, Equation (2) gives d S ( t ) / d t | t = t ≥
0, a contradiction.Suppose now that t > t , and choose t ∈ ( t , t ) such that F ( t ) <
0, d F ( t ) / d t | t = t <
0, and the delayssatisfy S ( t ) ≥ F ( t ) ≥
0. Then, since P ( t ) ≤ K by Lemma 3, and since S ( t ) ≥
0, Equation (2)gives d F ( t ) / d t | t = t ≥
0, a contradiction.Suppose now that t = t . Since S (0) > S ( t ) = 0, there exists t ∈ (0 , t ) such that d S ( t ) / d t | t = t <
0. Then, since P ( t ) ≤ K by Lemma 3, since the delays satisfy S ( t ) ≥ F ( t ) ≥
0, and since S ( t ), F ( t ) ≥
0, Equation (3) gives d S ( t ) / d t | t = t ≥
0, a contradiction. We conclude that S ( t ) ≥ t > β S − β F ≥ F ( t ) < t > t ∈ R . If t < t then, since F (0) > F ( t ) = 0, there exists t ∈ (0 , t ) such thatd F ( t ) / d t | t = t <
0. Then, since P ( t ) ≤ K by Lemma 3, since the delays satisfy S ( t ) ≥ F ( t ) ≥ S ( t ), F ( t ) ≥
0, Equation (3) gives d F ( t ) / d t | t = t ≥
0, a contradiction.Suppose now that t > t , and choose t ∈ ( t , t ) such that S ( t ) <
0, d S ( t ) / d t | t = t <
0, and the delayssatisfy S ( t ) ≥ F ( t ) ≥
0. Then, since P ( t ) ≤ K by Lemma 3, and since F ( t ) ≥
0, Equation (2)gives d S ( t ) / d t | t = t ≥
0, a contradiction.Suppose now that t = t . Since F (0) > F ( t ) = 0, there exists t ∈ (0 , t ) such that d F ( t ) / d t | t = t <
0. Then, since P ( t ) ≤ K by Lemma 3, since the delays satisfy S ( t ) ≥ F ( t ) ≥
0, and since S ( t ),6 ( t ) ≥
0, Equation (3) gives d F ( t ) / d t | t = t ≥
0, a contradiction. We conclude that F ( t ) ≥ t > β S − β F < F ( t ) < t > t ∈ R . If t < t then choose t ∈ ( t , t ) such that F ( t ) <
0, d F ( t ) / d t | t = t <
0, and the delayssatisfy S ( t ) ≥ F ( t ) ≥
0. Then, since P ( t ) ≤ K by Lemma 3 and since S ( t ) ≥
0, Equation (3)gives d F ( t ) / d t | t = t ≥
0, a contradiction.Suppose now that t ≥ t . Since S (0) > S ( t ) = 0 there exists t ∈ (0 , t ) such that d S ( t ) / d t | t = t <
0. Then, since P ( t ) ≤ K by Lemma 3, since the delays satisfy S ( t ) ≥ F ( t ) ≥
0, and since S ( t ), F ( t ) ≥
0, Equation (2) gives d S ( t ) / d t | t = t ≥
0, a contradiction. We conclude that F ( t ) ≥ t >
0. 7
Complete proof of Theorem 4
We first show that f defined in Equation (14) satisfies the following Lipschitz condition on every boundedsubset of C : for all M >
L > ρ , ψ ∈ C ([ − b U , , R ) with k ρ k , k ψ k ≤ M wehave k f ( ρ ) − f ( ψ ) k ≤ L k ρ − ψ k .To further simplify the notation in Equation (14) we define κ = (2 α S − r S , κ = 2(1 − α F ) r F , κ = 2(1 − α S ) r S , κ = (2 α F − r F , and κ = ( β S − β F ) /K . Now, f ( ρ ) − f ( ψ ) = (cid:16) κ ρ S + κ ρ F (cid:17)(cid:18) − ( ρ S + ρ F )(0) K (cid:19) − κ ρ S (0) ρ F (0) − (cid:16) κ ψ S + κ ψ F (cid:17)(cid:18) − ( ψ S + ψ F )(0) K (cid:19) + κ ψ S (0) ψ F (0) (cid:16) κ ρ S + κ ρ F (cid:17)(cid:18) − ( ρ S + ρ F )(0) K (cid:19) + κ ρ S (0) ρ F (0) − (cid:16) κ ψ S + κ ψ F (cid:17)(cid:18) − ( ψ S + ψ F )(0) K (cid:19) − κ ψ S (0) ψ F (0) = (cid:0) κ ( ρ S − ψ S ) + κ ( ρ F − ψ F ) (cid:1)(cid:18) − ( ρ S + ρ F )(0) K (cid:19) + (cid:16) κ ψ S + κ ψ F (cid:17)(cid:18) ( ψ S − ρ S )(0) + ( ψ F − ρ F )(0) K (cid:19) + κ ( ψ S − ρ S )(0) ψ F (0) + κ ( ψ F − ρ F )(0) ρ S (0) (cid:0) κ ( ρ S − ψ S ) + κ ( ρ F − ψ F ) (cid:1)(cid:18) − ( ρ S + ρ F )(0) K (cid:19) + (cid:16) κ ψ S + κ ψ F (cid:17)(cid:18) ( ψ S − ρ S )(0) + ( ψ F − ρ F )(0) K (cid:19) + κ ( ρ S − ψ S )(0) ψ F (0) + κ ( ρ F − ψ F )(0) ρ S (0) so, using the triangle inequality, we obtain k f ( ρ ) − f ( ψ ) k ≤ (cid:13)(cid:13)(cid:13)(cid:13)(cid:13)(cid:13)(cid:13)(cid:13)(cid:13)(cid:13) κ ( ρ S − ψ S ) (cid:18) − ( ρ S + ρ F )(0) K (cid:19) κ ( ρ S − ψ S ) (cid:18) − ( ρ S + ρ F )(0) K (cid:19) (cid:13)(cid:13)(cid:13)(cid:13)(cid:13)(cid:13)(cid:13)(cid:13)(cid:13)(cid:13) (cid:13)(cid:13)(cid:13)(cid:13)(cid:13)(cid:13)(cid:13)(cid:13)(cid:13)(cid:13) κ ( ρ F − ψ F ) (cid:18) − ( ρ S + ρ F )(0) K (cid:19) κ ( ρ F − ψ F ) (cid:18) − ( ρ S + ρ F )(0) K (cid:19) (cid:13)(cid:13)(cid:13)(cid:13)(cid:13)(cid:13)(cid:13)(cid:13)(cid:13)(cid:13) + (cid:13)(cid:13)(cid:13)(cid:13)(cid:13)(cid:13)(cid:13)(cid:13)(cid:13)(cid:13) κ ψ S (cid:18) ( ψ S − ρ S )(0) + ( ψ F − ρ F )(0) K (cid:19) κ ψ S (cid:18) ( ψ S − ρ S )(0) + ( ψ F − ρ F )(0) K (cid:19) (cid:13)(cid:13)(cid:13)(cid:13)(cid:13)(cid:13)(cid:13)(cid:13)(cid:13)(cid:13) + (cid:13)(cid:13)(cid:13)(cid:13)(cid:13)(cid:13)(cid:13)(cid:13)(cid:13)(cid:13) κ ψ F (cid:18) ( ψ S − ρ S )(0) + ( ψ F − ρ F )(0) K (cid:19) κ ψ F (cid:18) ( ψ S − ρ S )(0) + ( ψ F − ρ F )(0) K (cid:19) (cid:13)(cid:13)(cid:13)(cid:13)(cid:13)(cid:13)(cid:13)(cid:13)(cid:13)(cid:13) + (cid:13)(cid:13)(cid:13)(cid:13)(cid:13)(cid:13)(cid:13)(cid:13)(cid:13) κ ( ψ S − ρ S )(0) ψ F (0) κ ( ρ S − ψ S )(0) ψ F (0) (cid:13)(cid:13)(cid:13)(cid:13)(cid:13)(cid:13)(cid:13)(cid:13)(cid:13) + (cid:13)(cid:13)(cid:13)(cid:13)(cid:13)(cid:13)(cid:13)(cid:13)(cid:13) κ ( ψ F − ρ F )(0) ρ S (0) κ ( ρ F − ψ F )(0) ρ S (0) (cid:13)(cid:13)(cid:13)(cid:13)(cid:13)(cid:13)(cid:13)(cid:13)(cid:13) ≤ q κ + κ (cid:16) K | ρ S (0) | + 1 K | ρ F (0) | (cid:17) | ( ρ S − ψ S ) | + q κ + κ (cid:16) K | ρ S (0) | + 1 K | ρ F (0) | (cid:17) | ( ρ F − ψ F ) | + 1 K q κ + κ | ψ S || ( ψ S − ρ S )(0) + ( ψ F − ρ F )(0) | + 1 K q κ + κ | ψ F || ( ψ S − ρ S )(0) + ( ψ F − ρ F )(0) | + √ | κ | | ( ρ S − ψ S )(0) | | ψ F (0) | + √ | κ | | ( ρ F − ψ F )(0) | | ρ S (0) |≤ q κ + κ (cid:16) K k ρ S k + 1 K k ρ F k (cid:17) k ρ S − ψ S k + q κ + κ (cid:16) K k ρ S k + 1 K k ρ F k (cid:17) k ρ F − ψ F k + 1 K q κ + κ k ψ S k (cid:0) k ψ S − ρ S k + k ψ F − ρ F k (cid:1) + 1 K q κ + κ k ψ F k (cid:0) k ψ S − ρ S k + k ψ F − ρ F k (cid:1) √ | κ |k ρ S − ψ S kk ψ F k + √ | κ |k ρ F − ψ F kk ρ S k≤ (cid:18)q κ + κ (cid:16) MK (cid:17) + q κ + κ (cid:16) MK (cid:17) + 2 MK q κ + κ + 2 MK q κ + κ + 2 √ | κ | M (cid:19) k ρ − ψ k = (cid:18)q κ + κ + q κ + κ (cid:19)(cid:18) MK (cid:19) + 2 √ | κ | M ! k ρ − ψ k , so we can set L to be L = (cid:18)q κ + κ + q κ + κ (cid:19)(cid:18) MK (cid:19) + 2 √ | κ | M and then f satisfies the Lipschitz condition. Then [17, Page 32, Theorem 3.7] provides local existenceand uniqueness of solutions for the system (2) and (3). Since our solutions of interest are bounded byTheorem 2, it follows from [17, Page 37, Proposition 3.10] that the solutions are continuable to all positivetime. 10 Alternative proof of Proposition 6
The functions S and F are non-negative by Theorem 2, so the following inequalities show that theequilibrium point ( S * , F * ) = (0 , ∈ C ([ − b U , ∞ ) , R ≥ ) for the system (2) and (3) is stable if and only ifthe equilibrium point P * = 0 ∈ C ([ − b U , ∞ ) , R ≥ ) for Equation (12) is stable: k ( S, F ) k ≤ k P k ≤ √ k ( S, F ) k . (S5)The history function φ for the system (2) and (3) is strictly positive by Equation (9), and φ S + φ F < K byEquation (10). It follows by continuity that S ( t ) > F ( t ) >
0, and S ( t ) + F ( t ) < K in a neighbourhoodof (0 , P ( t ) / d t > P * = 0, which is therefore anunstable equilibrium point. We conclude that ( S * , F * ) = (0 ,
0) is also an unstable equilibrium point byEquation (S5). 11
Graphical illustration of Cauchy’s argument principle
In Figure S2 we graphically illustrate our application of Cauchy’s argument principle. We let Γ be theclosed rectangular contour in Figure S2(a). Due to the integrals in the transcendental characteristicequation G ( λ ) in Equation (25), a very large number of numerical integrations are required to calculate G ( λ ) along a contour. So, we instead use the discrete delay version of G ( λ ), denoted by G δ ( λ ), whichgives similar qualitative behaviour: G δ ( λ ) = λ − λ (cid:16) (2 α S − re − λr F τ S + (2 α F − e − λr F τ F (cid:17) + (2 α S + 2 α F − re − λr F ( τ S + τ F ) , (S6)obtained from G ( λ ) with the Dirac kernels g S ( z ) = δ ( z − τ S ) and g F ( z ) = δ ( z − τ F ) for discrete delays τ S and τ F . Figures S2(b), (d), and (f) show the images G δ (Γ) for three different sets of parameters for G δ ( λ ), and Figures S2(c), (e), and (g) show the respective close-up views around the origin. Note thateach coloured segment of G δ (Γ) in (b)–(f) is the image of the same-coloured segment of the rectangularcontour Γ in (a) under G δ ( λ ). To calculate the winding number of G δ (Γ) with respect to the origin wecount the net number of times that G δ (Γ) winds counter-clockwise around the origin, assigning +1 foreach time G δ (Γ) winds around the origin in a counter-clockwise direction, and − G δ (Γ)winds around the origin in a clockwise direction. Figure S2 illustrates that the behaviour of G δ (Γ) canbe complicated near the origin, so when calculating the winding number of G δ (Γ) with respect to theorigin we need to ensure that we also account for the possibility of G δ (Γ) winding around the origin in aclockwise direction. 12 igure S2: Graphical illustration of our application of Cauchy’s argument principle for G δ ( λ ) in Equation (S6),when 2 α S + 2 α F − >
0. (a) A rectangular contour Γ in the right half-plane does not intersect the imaginary axis.(b) The image contour G δ (Γ) when G δ has the parameters α S = 0 . α F = 1, r F = 1, and r = 0 .
01, with a close-upview of the origin in (c). (d) The image contour G δ ( γ ) when G δ has the parameters α S = 0 . α F = 1, r F = 0 . r = 1, with a close-up view of the origin in (e). (f) The image contour G δ (Γ) when G δ has the parameters α S = 1, α F = 1, r F = 1, and r = 0 .
01, with a close-up view of the origin in (g). Each coloured segment of G δ (Γ)in (b)–(f) is the image of the same-coloured segment of the rectangular contour Γ in (a) under G δ ( λ ) eferences [1] Vittadello ST, McCue SW, Gunasingh G, Haass NK, Simpson MJ. Examining go-or-grow usingfluorescent cell-cycle indicators and cell-cycle-inhibiting drugs. Biophysical Journal. 2020;118:1243–1247.[2] Haass NK, Beaumont KA, Hill DS, Anfosso A, Mrass P, Munoz MA, et al. Real-time cell cycleimaging during melanoma growth, invasion, and drug response. Pigment Cell & Melanoma Research.2014;27:764–776.[3] Davies MA, Stemke-Hale K, Lin E, Tellez C, Deng W, Gopal YN, et al. Integrated molecular andclinical analysis of AKT activation in metastatic melanoma. Clinical Cancer Research. 2009;15:7538–7546.[4] Hoek KS, Schlegel NC, Brafford P, Sucker A, Ugurel S, Kumar R, et al. Metastatic potentialof melanomas defined by specific gene expression profiles with no BRAF signature. Pigment CellResearch. 2006;19:290–302.[5] Smalley KSM, Contractor R, Haass NK, Kulp AN, Atilla-Gokcumen GE, Williams DS, et al. Anorganometallic protein kinase inhibitor pharmacologically activates p53 and induces apoptosis inhuman melanoma cells. Cancer Research. 2007;67:209–217.[6] Smalley KSM, Contractor R, Haass NK, Lee JT, Nathanson KL, Medina CA, et al. Ki67 expressionlevels are a better marker of reduced melanoma growth following MEK inhibitor treatment thanphospho-ERK levels. British Journal of Cancer. 2007;96:445–449.[7] Spoerri L, Beaumont KA, Anfosso A, Haass NK. Real-time cell cycle imaging in a 3D cell culturemodel of melanoma. Methods in Molecular Biology. 2017;1612:401–416.[8] Sakaue-Sawano A, Kurokawa H, Morimura T, Hanyu A, Hama H, Osawa H, et al. Visualizingspatiotemporal dynamics of multicellular cell-cycle progression. Cell. 2008;132:487–498.[9] Allen LJS. An Introduction to Stochastic Processes with Applications to Biology. 2nd ed. Taylor &Francis Ltd.; 2010.[10] Weber TS, Jaehnert I, Schichor C, Or-Guil M, Carneiro J. Quantifying the length and variance ofthe eukaryotic cell cycle phases by a stochastic model and dual nucleoside pulse labelling. PLoSComputational Biology. 2014;10:e1003616. 1411] Yates CA, Ford MJ, Mort RL. A multi-stage representation of cell proliferation as a Markov process.Bulletin of Mathematical Biology. 2017;79:2905–2928.[12] Vittadello ST, McCue SW, Gunasingh G, Haass NK, Simpson MJ. Mathematical models incorporat-ing a multi-stage cell cycle replicate normally-hidden inherent synchronization in cell proliferation.Journal of the Royal Society Interface. 2019;16:20190382.[13] Chao HX, Fakhreddin RI, Shimerov HK, Kedziora KM, Kumar RJ, Perez J, et al. Evidence thatthe human cell cycle is a series of uncoupled, memoryless phases. Molecular Systems Biology.2019;15:e8604.[14] Gavagnin E, Ford MJ, Mort RL, Rogers T, Yates CA. The invasion speed of cell migration modelswith realistic cell cycle time distributions. Journal of Theoretical Biology. 2019;481:91–99.[15] MATLAB lsqnonlin. Solve nonlinear least-squares (nonlinear data-fitting) problems (R2019b). Ac-cessed February 2020.; 2019. Available from: https://mathworks.com/help/optim/ug/lsqnonlin.htmlhttps://mathworks.com/help/optim/ug/lsqnonlin.html