On the characteristic equation λ = α1 + (α2 + α3λ)e(-λ) and its use in the context of a cell population model.

In this paper we characterize the stability boundary in the (α1, α2)-plane, for fixed α3 with −1 < α3 < +1, for the characteristic equation from the title. Subsequently we describe a nonlinear cell population model involving quiescence and show that this characteristic equation governs the (in)stability of the nontrivial steady state. By relating the parameters of the cell model to the αi we are able to derive some biological conclusions.

Introduction

The characteristic equationfor , with parameters , corresponds to the delay differential equationThe fact that (1.1) can have, for , solutions with therefore leads to the dictum that delayed negative feedback can lead to oscillatory instability.

Equation (1.1) can be analyzed in great detail (see e.g. Section III.3 in Èl’sgol’ts and Norkin 1973; Hayes 1950; Chapter 13.7 in Bellman and Cooke 1963; Chapter XI in Diekmann et al. 1991). The outcome can be conveniently summarized in the diagram depicted in Fig. 1 (corresponding to what Èl’sgol’ts and Norkin 1973 call the method of D-partitions; also see Breda 2012).

Fig. 1

The numbers specify the number of roots of (1.1) with for in the corresponding region of the parameter plane

In the context of specific models, (1.2) usually arises by linearization of a nonlinear equation around a steady state. In such a situation are (sometimes rather complicated) functions of the original model parameters. As emphasized in Chapter XI in Diekmann et al. (1991) and the recent didactical note (Diekmann and Korvasova 2013), one can combine the detailed knowledge embodied in Fig. 1 with an analysis of the parameter map, that relates the original parameters to , in order to obtain stability and bifurcation results for the steady state of the nonlinear equation. Even if the aim is to perform a one-parameter Hopf bifurcation study, it is much more efficient to derive first the two parameter picture of Fig. 1, see Diekmann and Korvasova (2013). We also refer to Insperger and Stépán (2011), Kuang (1993), Michiels and Niculescu (2014), Stépán (1989) for a systematic approach to deriving stability conditions for steady states of delay differential equations via an analysis of characteristic equations and to Cheng and Lin (2009) for potentially useful general theory.

In Sect. 3 we shall formulate a model for a population of cells that can either be on the way to division or quiescent. The model is very close in spirit to the one formulated and studied by Gyllenberg and Webb in their pioneering paper (Gyllenberg and Webb 1990). But our formulation is in terms of delay equations, more precisely a system consisting of one Renewal Equation (RE) and one Delay Differential Equation (DDE) (Diekmann and Gyllenberg 2012; Diekmann et al. 2010; 2007/2008), rather than in terms of PDE as in Gyllenberg and Webb (1990). The delay equation formulation enables a relatively painless derivation of the characteristic equation whose roots govern the (in)stability of the (unique) nontrivial steady state. In the relatively simple case considered here (see Alarcón et al. 2014; Borges et al. 2014 for a more general model formulation), the RE degenerates into a difference equation in continuous time! The characteristic equation corresponding to the linearization about the nontrivial steady state takes the formThe goal of the present paper is to investigate how the picture of Fig. 1 deforms if we let grow away from zero (but restrict to , for reasons explained around Lemma 2.2 below). See the supplementary material for a movie showing how the curves depicted in Fig. 1 move in the plane when ranges from to . Our main conclusion is that, from a qualitative point of view, nothing changes at all.

The stability region in the -plane

Our aim is to characterizeby a precise description of its boundary. We shall show in Theorem 2.1 that this boundary consists of the half-linecorresponding to being a root of (1.3), and the curvewhere corresponding to , being a root of (1.3). Note that and intersect in , corresponding to being a double root, see Lemma 2.4 below. The stability region is depicted in Fig. 2.

Fig. 2

Stability region for the characteristic equation (1.3). The boundary of consists of the half-line and the curve

Theorem 2.1

Let . is the connected open subset of that contains and has as its boundary.

The proof has quite a few components, so we build it up gradually by deriving auxiliary results. Our first step is to exclude that roots can enter the right half plane at when and are varied. This leads to the constraint .

Lemma 2.2

Let and let be a root of (1.3) with . Then

Proof

We write (1.3) in the formand take absolute values at both sides. Usingwe obtain the estimate (2.5).

To further illustrate the importance of the restriction , we formulate a corollary of Theorem 1.1 in Chapter 3 of Kuang (1993):

Lemma 2.3

When , Eq. (1.3) has infinitely many roots in the right half plane .

We refer to Lemma 2.12 below for a description of the limiting behaviour for and . In order to facilitate various formulations, we adopt the convention that in the rest of this section the inequalitieshold, unless stated otherwise.

Lemma 2.2 implies (by way of Rouché’s Theorem, see Lemma 2.8 in Chapter XI of Diekmann et al. 1991) that at the boundary of the Eq. (1.3) has either a root or a pair of roots .

Lemma 2.4

is a root of (1.3) iff . It is a double root iff in addition (and a triple root iff in addition ).

is a complex conjugate pair of roots of (1.3) iff with defined by (2.4).

.

Substituting into (1.3) we obtain . Taylor expanding we can write (1.3) in the form and conclude that for the root is a double root iff and a triple root iff in addition .

In terms of real variables and such that the complex equation (1.3) amounts to the two real equations , where corresponds to the real part and is given by while corresponds to the imaginary part and is given by For the equations reduce to Solving the second for we obtain . Next we can solve the first for . This yields .

Since and , it is clear that and for .

Our next step is to identify for any at least one point in the -plane that belongs to . Clearly Eq. (1.3) is very simple if , but this point lies on . The idea is to focus on a neighborhood, i.e., first consider and next perturb.

Lemma 2.5

For the roots of (1.3) are given by and, if , while for , is the one and only root.

Consider and small. The only root of (1.3) that can possibly lie in the right half plane is given by As a consequence, the points with , small, belong to .

For Eq. (1.3) reduces to . So is a root and all other roots satisfy .

So for we have one root on the imaginary axis and, if , countably many other roots that are at a uniformly positive distance away from the imaginary axis in the left half plane. According to the Implicit Function Theorem, there exists for small a root if we keep , and this is the only root in a small ball B around . According to a variant of Rouché’s Theorem, see Lemma 2.8 of Chapter XI of Diekmann et al. (1991), there are for small no roots in the open set (Note that here we use Lemma 2.2 to compensate for the non-compactness of the closure of this set.)

As we saw in Part 2 of the lemma above, it is helpful to know how the root moves in if we move away from the line . The following result shows that if we cross transversally, away from its endpoint, from below to above, a real root of (1.3) crosses the imaginary axis from left to right.

Lemma 2.6

Suppose that , for some , is with for some . Then (1.3), for small , has a real root , which is of the formIf additionally and , then for small negative and for small positive .

We omit the elementary proof.

In a similar spirit, we want to know how the roots move in if we move away from the curve . In Section XI.2 of Diekmann et al. (1991) it is explained that the crucial quantity is the sign of , wherewith defined by (2.8) and by (2.9). So in the present case we haveand for . In order to define what we mean by “to the left” and “to the right”, we provide with the orientation of increasing . According to Proposition 2.13 in Chapter XI of Diekmann et al. (1991), we can draw the following conclusion.

Lemma 2.7

When we cross from right to left at a point corresponding to , a pair of complex conjugate roots of (1.3) crosses the imaginary axis from left to right.

Motivated by Lemma 2.4(2) we define for the intervals and the curveswith defined by (2.4). In the next lemma we collect some observations that help to understand the shape of the curves .

Lemma 2.8

For defined by (2.4) we have

for

is periodic with period and decreases from to on and increases from to on and on,

Now note that and that Since is a linear function of , it follows that for the inequality holds for .

At with m odd, switches from positive to negative values and . For even m, switches from negative to positive values, but .

Since , the sign of is the opposite of the sign of .

is obviously -periodic. so is an increasing function of on intervals on which and a decreasing function of on intervals on which . Moreover

Note that Lemma 2.8(1) implies that the curves , can also be parameterized by . If we do so for , we can in Lemma 2.7 replace “left” by “below” and “right” by “above”, just like in Lemma 2.6.

Lemma 2.9

There are no intersections of the curves .

The curves are situated in the quarter plane while the curves and are situated in the quarter plane except for the starting point of which lies on the boundary.

The curves are ordered as shown in Fig. 1.

Assume that for some and it holds that then necessarily i.e., Multiplying this identity by we find that necessarily . By substituting in the identity we deduce that in addition the identity has to hold. But also implies that , so we arrive at the conclusion that . Since we parameterize with , the only possibility is .

From Lemma 2.8(4) we know that for If this amounts to , i.e., and . If it amounts to , i.e., and . In case of the inequality holds for but not for , where the inequality becomes equality.

Intersections of the curves with the -axis are characterized by cf (2.4a). So the corresponding value of equals with for and and for . Accordingly the ordering is determined by so by k.

At last we are ready to present the

Proof of Theorem 2.1

From Lemma 2.5(2) we know that at least a part of the negative -axis belongs to . As noted just before Lemma 2.4, the a priori estimate of Lemma 2.2 provides the compactness needed to apply Rouché’s theorem in order to deduce that is to a large extent characterized by (1.3) having a root on the imaginary axis. In other words, is contained in the union of the line and the curves , , cf. Lemmas 2.4(1) and 2.4(2) and the definitions (2.3) and (2.11). On account of Lemma 2.9 we can now conclude that the connected component of containing the negative -axis has as its boundary.

In principle might have other components. Indeed, roots that enter the right half plane, when crossing or , could return to the left half plane as one of the curves is crossed. There are at least two ways to see that, actually, this does not happen. The first is by extending Lemma 2.7 to the curves ; more precisely, by noting that the sign of is opposite to the sign of and that accordingly more and more roots move into the right half plane if we keep decreasing after crossing or if we keep increasing after crossing L. So also the numbers shown in Fig. 1 extend to .

The second way is by showing that the roots that enter the right half plane upon crossing remain “caught” in the strip , so cannot possibly return to the left half plane when one of the curves is crossed. Indeed, let be a root of (1.3) and assume that and . Then necessarily and hence the equation (recall (2.9)) readsshowing that either or . But since we have and since we have , so is not possible. (What can, and does happen though, is that they become real and that subsequently one of the two real roots returns to the left half plane when the parameter point crosses the line from below to above at a point with , cf. Lemma 2.6.)

We conclude that (1.3) has a root with if is in the connected open subset of that has as its boundary and that does not contain the negative -axis.

It follows from Theorem 2.1 and Lemma 2.9(2) that the setbelongs to for (the index u is meant to express “uniformly” in ). By analyzing the limiting behavior of for and , we shall show that this is sharp, i.e., is the largest set that belongs to for . Again we need some auxiliary results.

Lemma 2.10

For the inequality holds for .

For there exists such that

for .

with Now note that and that when .

If we find from the above expression for that h increases for small positive and hence takes positive values for small positive . Since , h changes sign at least once on and the exact number of sign changes of h must be an odd number. If h changes sign three or more times, the derivative must also change sign at least three times. Now requires that . We claim that whenever then necessarily . It follows that can change sign only once on , that the same holds for and, therefore, for h. So h changes sign exactly once and since so does . It remains to prove the claim. First observe that since the left hand side equals zero for and . Next note that So if then

and both factors are positive on .

For the point , where and meet, moves to (0, 0) and any point on moves to on the line . The fact that for one cannot interchange the limits and is irrelevant, since the lines and intersect in (0, 0). Note that decreases strictly from 0 to as increases from 0 to , see Lemma 2.8(1).

For , the point moves to and any point , on moves to on the line . In this case, one cannot interchange the limits and . For slightly above , there is a small interval such that moves from close to (0, 0) for all the way to for , see Lemma 2.10. So the limit set of is the union of and . Or, in other words, converges to , but if we use to parameterize the part of the boundary, the convergence is rather non-uniform. We summarize our conclusions in

Theorem 2.11

If for all in an interval of either the form or of the form for some , then .

For completeness, we formulate a result about the limiting behaviour of the curves . As a first indication of what to expect, recall from the proof of Lemma 2.9(3) that the intersection with the -axis hasand conclude that the intersection point converges to (0, 0). As in the proof of Lemma (2.10)(2) one can show that changes sign exactly once on defined in (2.10). So geometrically it seems clear that the bounds from Lemma 2.9(2) become sharp in the limit as converges to 1. We now prove that this is indeed the case.

Lemma 2.12

When either or , the curves all converge towhile the curves all converge to

We provide the proof for and , in order to illustrate some details more clearly. We trust that the reader believes that all the other cases are covered by essentially the same arguments.

The curve is parameterized by . So and increases from to . From (2.4) we conclude that for any fixed both and converge to as . So converges to a point on the half-line , . Since decreases from for to 0 for , any point on this line is in fact the limit of points on as .

Let r be an arbitrary negative real number. By Lemma 2.8(2) we know that exists such that . Note that we must have for . Hence as . So let us put where, of course, . Using Taylor expansion of and , we find from the equation that . Next the expression for yields . We conclude that converges to as .

The limiting stability diagram is depicted in Fig. 3. Lemma 2.12 and this diagram are in complete accordance with the results that one obtains by studying (1.3) directly for . The two real equationsyield that for while putting yieldsSo requires k to be even and then while requires k to be odd and then .

Fig. 3

Partitioning of the -parameter plane according to the number of roots of (1.3) in the right half plane for the limiting cases and

A cell population model involving quiescence

Motivated by Alarcón et al. (2014), we assume that the cell cycle incorporates a checkpoint for the prevailing environmental condition in the sense that, upon passage of this point, a cell commits itself to division with probability , while going quiescent with probability . If we allow that one can interpret as the probability that the cell undergoes checkpoint triggered apoptosis. Here, however, we shall assume that and that .

Quiescent cells have a probability per unit of time to reactivate and then progress towards division. Once reactivated, they differ in no way from the cells that never went quiescent. So quiescence amounts to a variable delay determined by the course of the environmental condition and an element of chance.

We assume that the proliferation step takes a fixed amount of time , in the following sense:A dividing cell produces two daughter cells. But we allow for a uniform death rate and accordingly the expected number of progeny arriving at the checkpoint after time equals . See Alarcón et al. (2014) for a variant that incorporates a more general probability distribution for cell cycle duration as well as a more general survival probability. Also see Adimy and Chekroun (Adimy and Chekroun) for a detailed analysis of the situation that does not depend on E and for references to various papers on hematopoietic cell models, many of them originating from the pioneering work of Mackey (1978).

if a cell commits itself to proliferation at the checkpoint, its (surviving) offspring arrives at the checkpoint after time

if a quiescent cell is reactivated, it likewise takes exactly time for its offspring to arrive at the checkpoint.

Let Q(t) denote the quantity of quiescent cells at time t. Let p(t) denote the quantity of cells that, per unit of time, set out on division at time t. The mathematical formulation of the assumptions described above takes the form

To complete the model formulation, we need to specify the dynamics of E and, in particular, the feedback law that describes the impact of the cell population on the environmental condition.

For concreteness, think of E as oxygen concentration. The equationexpresses that E is determined by the balance of inflow, outflow and consumption. If the constants and are all big relative to and the range of and G, we can make a quasi-steady-state-assumption and replace (3.2) by the explicit expression for E, in terms of Q and the history of p, obtained by putting the right hand side of (3.2) equal to zero and solving for E. Before writing this expression, let us discuss the issue of scaling and the concomitant reduction in the number of parameters.

By scaling of time, we can achieve that . Since the system (3.1) is linear in p and Q, it does not change when we scale both of these variables with the same factor. A particular choice of this factor amounts to replacing by and by with . Finally, in (3.1) the variable E only occurs as an argument of and G and we have not yet specified these functions. So scaling E by the factor does not lead to any loss of generality. The upshot is the scaled system that is based on the quasi-steady-state-assumption. The remaining parameters are and the two functions G and (note carefully that throughout the rest of the paper ).

Our main general results concerning (3.3) areWe use the systematic methodology of Diekmann et al. (2003) to derive Theorem 3.1. Linearization in the trivial steady state amounts to putting in the equations for p and Q. So the linearized system reads (when, abusing notation, we use the same symbols to denote the variables) In order to reveal the link between the (in)stability of the trivial steady state and the (non)existence of a nontrivial steady state, consider, for fixed positive E, the linear system

Theorem 3.1 below about existence and uniqueness of a nontrivial steady state

(in)stability of the nontrivial steady state is determined by the position in of the roots of (1.3) with given by (3.20) below

At steady state the value E should be such that (3.5) has a nontrivial (and positive) steady state. If it has, it has a one-parameter family (since (3.5) is a linear system). The parameter should then be tuned so that E, when expressed by the steady state version of (3.3c), does have the required value.

So we study the linear system (3.5) with parameter E. The simplest approach is based on the biological interpretation and runs as follows. In a constant environment a cell that goes quiescent has probabilityto become reactivated (rather than dying while quiescent). Hence the overall probability that a cell arriving at the checkpoint sets out on division, equalsAccordingly the expected number of progeny arriving at the checkpoint is given by the formulaOn the basis of (3.4) we now conclude that the trivial steady state is stable if and unstable if .

If eitherorwe find (by calculating the derivative of the right hand side of (3.6)) that is a strictly increasing function of E.

In that case, there is at most one root for the equationThe Eq. (3.3c) shows that only roots in [0, 1] yield candidates for steady states. The assumptionguarantees that the unique root, if it exists, belongs to [0, 1). The existence is indeed guaranteed if, in addition to (3.8), we assume thatIncidentally, note that cannot possibly hold if , so by imposing (3.8) we implicitly require that . Also note that from (3.6) it follows that and hence we automatically haveif .

We denote by the unique solution of (3.7). Let us look for constant solutions of (3.5). These should satisfyThe determinant of the matrix is equal toso we see right away that the condition for is both necessary and sufficient for a solution to exist. Then the vectoris, for any real number K, in the null space of the matrix, so is a constant solution of (3.5). Finally, we close the feedback loop by requiring thatSince iffholds, one haswhere we defineWe then find thatUsing (3.11)–(3.14) one obtains and asWe summarize the discussion above in

Theorem 3.1

Let and G be continuously differentiable functions such that, for ,Assume that, with defined by (3.6),Then Eq. (3.7) has a unique solution . Moreover, (3.3) has a unique steady state with p and Q given by (3.15).

Remark 3.2

It is biologically obvious that if the population goes extinct, while if it grows beyond any bound, despite the negative feedback via and/or G. A proof is provided in the Appendix.

Our next aim is to study the stability of the nontrivial steady state. The first step in this direction consists of linearization of system (3.3). We putand focus on smallx, y. The Eq. (3.3c) implies thatFor we therefore haveThus we deduce that the linearized system is given byDefineThen from (3.15)Recall that (and hence ).

The characteristic equation has the formwhere is the matrix with entries

A straightforward calculation now establishes that (3.18) is exactly of the form (1.3) with

Note that (recall (3.10)). A straightforward computation, using (3.12) and (3.13), establishes thatAs observed in between (3.9) and (3.10) the first factor is negative. The expression (3.17) shows that is positive if . We conclude that for and that corresponds to the transcritical bifurcation at which (and hence ) loses its stability (in principle could yield saddle-node bifurcations as well, but since the internal equilibrium is unique, these do not actually happen).

The impact of (a one parameter study)

Proliferating and quiescent cells consume, in the model, oxygen in the proportion . Here we allow in principle , but a more reasonable range is (since if quiescent cells would consume more oxygen than proliferating cells, there does not seem to be any benefit in going quiescent). The aim of this section is to investigate the impact of on the stability of the nontrivial steady state.

As does not depend on , neither does . The formulas (3.20c) and (3.21) above next show that both and remain constant when is varied. Since (recall (3.13)), decreases if is increased, so increasing amounts to moving North-West in the -plane along the line defined by (3.21).

Proposition 4.1

For the nontrivial steady state is asymptotically stable.

By inserting in (3.20b) we find . Since , the corresponding point is located in the 2nd quadrant of the -plane, which belongs to . If we increase we move along the line defined by (3.21) away from its intersection with defined by (2.3), so we remain in .

The intersection of the line (3.21) and the curve is found by solving the equationfor . Note that decreases from 0 for to for . Hence, since , a unique solution exists.

Once the root of (4.1) is found, we can compute the corresponding value of by inserting the root into (2.4a). Let us call the result . The value of corresponding to is, according to (3.20a), given by . This provides us with a simple test:In the next section we show that it is far more efficient (as well as far more illuminating) to implement a two parameter version of this test. For now we just conclude that increasing promotes stability of the nontrivial steady state.

if the nontrivial steady state is asymptotically stable for .

if there exists such that the nontrivial steady state is asymptotically stable for but unstable for .

Two-parameter case studies

The model described in Sect. 3 involves two regulatory mechanisms: the partitioning between quiescent and proliferating cells at the checkpoint, as described by the dependence of on E, and the duration of the quiescent period, as described by the dependence of G on E. We will study these mechanisms separately and in turn, meaning that we first consider the case that does not depend on E and next the case that G does not depend on E. We use respectively and as a parameter, which measures how strongly the system reacts locally near the steady state to changes in the environmental condition.

The feedback to the environmental condition is described by Eq. (3.3c). The parameter captures the relative impact on the environmental condition of proliferating and quiescent cells. We will take as the second parameter.

As shown below, in both cases , as defined by (3.20c), does not depend on . Our strategy will be to use (3.20a) and (3.20b), and the definition of A and , in order to express (or ) and in terms of and . These expressions then allow us to picture the stability boundary in the -plane by mapping and as defined in, respectively, (2.2) and (2.3), to the -plane.

The following observations are inspired by the analysis of Sect. 4. The shape of shows that instability is promoted by letting increase and decrease i.e., by moving South-East in the -plane. If then defined by (3.20b) decreases if increases. Such an increase of has, for , the effect that decreases as well. But for only slightly bigger than the effect on is much larger than the effect on , so we might expect that making bigger leads to instability. Next (3.17) tells us that having either or large in is expected to promote instability. In a nutshell: we expect that steep response promotes instability.

Regulation via the length of the quiescent period

We assume that is independent of E. It follows right away from (3.20c) that is a constant depending on and , but not on . We fix and such that (3.10) holds and assume that the strictly increasing function G is such thatholds, which is equivalent to assuming (3.16). Using (3.6) and (3.20c) we write (3.7) in the formand accordingly defineIf we eliminate A from the pair of equations (3.20a) and (3.20b), we obtain the identitySolving for we findwhich upon substitution of (5.1) becomes an explicit expression for in terms of , , and .

When does not depend on E, (3.17) simplifies toand if we substitute this into (3.20a) and solve for we obtainwhich, by (5.1) and (5.2) is an explicit expression in , , , , and the function G.

For given and G the expressions (5.3) and (5.4) allow us to transform the stability boundary from in the plane to the plane. The result, depicted in Fig. 4a for and , clearly demonstrates that a steep response promotes instability. Of course is not really a free parameter and as a consequence Fig. 4 should be regarded as an illustration of a general phenomenon. In order to obtain further biological insights, one would need to consider the dependence on mechanistic parameters that shape the function G.

Fig. 4

a Stable and unstable parameter regions in the -plane for the case of regulated duration of quiescence. The parametric curve in the -plane is transformed to the outermost curve in the -plane (stability boundary). The curve inside the instability region corresponds to . The equilibrium becomes unstable for small and large . b Graph of the imaginary part along the stability boundary. On the dashed curve while on the continuous curve

Regulation via the fraction of cells that become quiescent

In our second case study we assume that G is independent of E. Recalling (3.6) we can write (3.7) in the formSince the right hand side does not depend on or , we conclude from (3.20c) that once again is a constant (now determined by and G). Then we getEliminating from (3.20a) and (3.20b) we find exactly as before the expression (5.3) for . Solving (3.20a) for we obtainwhich only differs from the right hand side of (5.4) by a factor fully determined by .

We conclude that, apart from a scaling of the axis, the stability domain is exactly the same as in the situation of Sect. 5.1. In other words, these stability considerations do NOT yield any information that, as we had hoped when embarking on our investigation, helps to decide on the basis of observable fluctuations whether the regulation works via initiation or via termination of quiescence. So it remains a wide open question whether model considerations are at all useful when trying to decide about this issue?

Discussion

Physiologically structured population models lead to delay equations (Diekmann et al. 2010, this paper). A first step in the subsequent analysis is, as a rule, an investigation of the dependence of steady states, and in particular their stability, on parameters. Thus characteristic equations enter the scene.

Characteristic equations come, in a sense, with their own natural parameters. As emphasized in Diekmann et al. (1991), and more recently echoed in Diekmann and Korvasova (2013), it is attractive to single out two parameters and determine curves in the two parameter plane corresponding to roots lying on the imaginary axis, so to roots being critical, i.e., neither contributing to stability nor to instability. If the characteristic equation has more than two parameters, this leads to two dimensional slices of a higher dimensional parameter space. With a bit of luck one can sometimes understand the full picture in terms of parameterized families of two dimensional sections. Here we have been lucky indeed.

The natural parameters of the characteristic equation are themselves functions of the model parameters. So after one has obtained a picture in terms of the natural parameters, it still remains to analyse how natural parameters change when model parameters change. In Sect. 4 we did exactly this: we studied how and change when varies between 0 and 1. But an efficient and attractive alternative is to single out TWO model parameters and to map the stability boundary to the corresponding plane of model parameters, as indeed we did in Sect. 5.

Our motivation for studying the characteristic equation (1.3) came from the cell model described in Alarcón et al. (2014) and in Sect. 3. So it was tempting to organise our results differently, in particular to begin with the model, derive the characteristic equation and then study it. Our decision to put, instead, the characteristic equation itself in the spotlight is rooted in the belief that this characteristic equation arises in other contexts and that, once the analysis in terms of natural parameters has been done, other applications only require the study of the map sending model parameters to natural parameters and its inverse. In other words, we hope (and expect) that Sect. 2 is the most useful part of the paper.

Both the initiation of quiescence and the termination of quiescence involve environmental signals. Here we have assumed that one and the same signal is involved, viz. the concentration of an essential resource like oxygen. Consumption of the resource creates a nonlinear feedback loop. Does it matter whether regulation occurs via initiation of quiescence or via termination of quiescence? Is it possible to infer from observed dynamics which of the two mechanisms is the dominant one? Here we have shown that this is not easy, if possible at all. Indeed, we found that destabilization of the steady state hinges onbut is rather independent of the precise way in which the feedback acts.

little difference in oxygen consumption of proliferating and quiescent cells

steep response to differences in oxygen concentration near the steady state value

In the present paper we have neither discussed the well-posedness of system (3.3) nor the justification of the Principle of Linearized Stability. Both Alarcón et al. (2014) and Borges et al. (2014) deal with these issues for distributed delay variants of the model. In Alarcón et al. (2014) it is also shown that one can take the limit in the characteristic equation for the delay kernel tending to a Dirac mass and arrive at the characteristic equation considered here. But the limit is not yet considered for the delay equations themselves. In work in progress, S. M. Verduyn Lunel and O. Diekmann are considering duality for neutral equations and this will probably make it unnecessary to consider limits. Concerning the Principle of Linearized Stability, the recent Diekmann and Korvasova (2016) does provide inspiration, but details certainly require attention.

Supplementary material 1 (gif 454 KB)