Mathematical Biology

J. Math. Biol. DOI 10.1007/s00285-015-0874-3

Discrete limit and monotonicity properties of the Floquet eigenvalue in an age structured cell division cycle model Stéphane Gaubert1 · Thomas Lepoutre2

Received: 25 February 2014 / Revised: 12 March 2015 © Springer-Verlag Berlin Heidelberg 2015

Abstract We consider a cell population described by an age-structured partial differential equation with time periodic coefficients. We assume that division only occurs within certain time intervals at a rate κ for cells who have reached minimal positive age (maturation). We study the asymptotic behavior of the dominant Floquet eigenvalue, or Perron-Frobenius eigenvalue, representing the growth rate, as a function of the maturation age, when the division rate κ tends to infinity (divisions become instantaneous). We show that the dominant Floquet eigenvalue converges to a staircase function with an infinite number of steps, determined by a discrete dynamical system. This indicates that, in the limit, the growth rate is governed by synchronization phenomena between the maturation age and the length of the time intervals in which division may occur. As an intermediate result, we give a sufficient condition which guarantees that the dominant Floquet eigenvalue is a nondecreasing function of the division rate. We also give a counter example showing that the latter monotonicity property does not hold in general. Keywords Cell cycle · Circadian rhythms · Structured PDEs · Perron-Frobenius theory · Floquet eigenvalue · Delay differential equations Mathematics Subject Classification

B

Primary 58F15 · 58F17; Secondary 53C35

Thomas Lepoutre [email protected] Stéphane Gaubert [email protected]

1

INRIA and CMAP, École Polytechnique, CNRS UMR 7641, CMAP, École Polytechnique, 91128 Palaiseau cedex, France

2

INRIA and Université de Lyon, Université Claude Bernard Lyon 1, CNRS UMR 5208, Institut Camille Jordan, 43 blvd. du 11 novembre 1918, 69622 Villeurbanne cedex, France

123

S. Gaubert, T. Lepoutre

1 Introduction 1.1 Population models in periodic environment We first review the motivations of this work. In population dynamics, one often needs to take into account the effect of the time-heterogeneity of the environment. There has been for a long time a great interest in the effect of seasonality, arising in particular from the question of harvesting (see Park et al. 1998 for instance). For models with time periodic coefficients, there has also been a more recent regain of interest for the question of an adapted definition of the classical parameter R0 and the question of epidemic threshold in particular by Bacaër (2007) and Bacaër and Guernaoui (2006), Bacaër and Ouifki (2007), Bacaer and Abdurahman (2008). However, our primary motivation originates from the study of biological systems with a smaller scale and a smaller period: the tissue scale and the period of a day. It has been established for years that living organisms exhibit so-called circadian rhythms (rhythms of around one day) (De Mairan 1729). In particular, the cell cycle is governed by circadian rhythms, which turned out to be of major interest in cancerology. Filipski et al. (2002, 2005) showed in mice that disruption (through abnormal light and food scheduling) of the circadian rhythms had a huge influence on tumor growth (in a nutshell, the more severe the disruption is, the faster the tumor grows). Moreover, one can take advantage of circadian rhythms to deliver drugs following a periodic schedule. This is the principle of chronotherapy (Mormont and Levi 2003). Based on these facts, Clairambault, Perthame and the authors (Clairambault et al. 2003, 2006, 2007a, b, 2009, 2011) studied models of the cell cycle in which a periodic forcing represents both the circadian rhythm and the periodic treatment. Those models were introduced in Clairambault et al. (2003, 2007b) and may be written as follows. We consider I phases of the cell cycle, in which the cells are represented with an age structure. The density of cells in phase i at time t with age x, denoted by n i (t, x), satisfies ⎧ ∂t n i (t, x)+∂x n i (t, x)+[di (t, x)+ K i→i+1 (t, x)]n i (t, x) = 0, x > 0 i < I ⎪ ⎪ ⎪ ⎨ ∂t n I (t, x) + ∂x n I i(t, x) + [d I (t, x) + K I →1 (t, x)]n 1 (t, x) = 0, x > 0, ∞ ⎪ n 1 (t, 0) = 2 0 K I →1 (t, x)n I (t, x)d x, ⎪ ⎪ ∞ ⎩ n i (t, 0) = 0 K i−1→i (t, x)n i−1 (t, x)d x, i > 1.

(1)

The number of phases I is typically 4. The situation in which I ∈ {1, 2}, which arises when phases of the detailed model are aggregated, is also of interest. Here, the cells in phase i are aging with time and may die at rate di or go to the next phase (coefficient K i→i+1 ). For the last phase of the cell cycle, the coefficient K I →1 represents the division rate. The boundary terms are representing the reintroduction of such cells into the next phase of the cell cycle. All the time dependencies are taken to be T −periodic (the period T representing the day). As the model is linear, it exhibits an exponential growth, see (Michel et al. 2005; Perthame 2007) for more information. More precisely, it is established that solutions behave as

123

Discrete limit and monotonicity properties …

n(t, x) ∼ N (t, x)eλ F t , as t → ∞

(2)

where N is the nonnegative eigenvector associated to the principal eigenvalue of the system: ⎧ ⎪ λ N (t, x)+∂x Ni (t, x)+[di (t, x)+ K i→i+1 (t, x)]Ni (t, x) = 0, x > 0 i < I ⎪ ⎪ F i ⎪ ⎪ ⎪ ⎨λ F N I (t, x) +∂x N I i(t, x) + [d I (t, x) + K I →1 (t, x)]N1 (t, x) = 0, x > 0, ∞ N1 (t, 0) = 2 0 K I →1 (t, x)N I (t, x)d x, ⎪  ⎪ ∞ ⎪ ⎪ ⎪ Ni (t, 0) = 0 K i−1→i (t, x)Ni−1 (t, x)d x i > 1, ⎪ ⎩ N  0, N (t + T, x) = N (t, x). i i i

(3)

The scalar λ F , which is known as the Floquet eigenvalue, determines the asymptotic behavior of the system, see (2). We refer to (Michel et al. 2005; Clairambault et al. 2006) for the question of existence and uniqueness of eigenelements. This model was initially designed in order to understand the experiment of Filipski et al., showing the influence on the growth rate of the loss of periodicity of the coefficients. This led to a series of results on the comparisons of periodic models and models with averaged coefficients over a time period (Clairambault et al. 2006, 2007a, 2009, 2011), in the framework of renewal models. One of the conclusions of the results was that the classical Kingman log-convexity inequality for the spectral radius of nonnegative matrices could be extended to such models. In particular, consider a renewal equation with periodic coefficients ⎧ ∂n ∂n ⎪ + + d(t, x)n(t, x) = 0, ⎪ ⎪ ⎨ ∂t ∂x  ∞ n(t, 0) = B(t, x)n(t, x)d x, ⎪ ⎪ 0 ⎪ ⎩ n(0, x) = n 0 (x) given

(4)

(we restrict to a single population instead of a system for readability). Then the growth rate was shown to depend on death rates d in a convex way and on birth rate B in a geometrically convex way (convex dependence on log B). This is in favor of the interpretation of the success of chronotherapy more as a toxicity reduction than an increase of efficiency (the latter being then interpreted as an indirect effect). We need also to recall the use of similar linear models, approached by delay equations techniques, by Bernard and Herzel (2006), Bernard et al. (2010). The second of these references deals with therapy optimization. We also refer to Thieme (1984), Diekmann et al. (1986), Magal (2001) for work on renewal models with periodic coefficients. The question of optimization of the Floquet eigenvalue has been addressed on such models by Billy et al. (2013), and interestingly on more complex models (size structured models) by Coron et al. (2015), where the concerned ’division’ is the breakup of polymers in prion polymerization model. In this case the time periodicity is not natural but the effect of the experimental device (PMCA, which consists of alternating period of high and low fragmentation). In this case, even without time heterogeneity, the eigenprob-

123

S. Gaubert, T. Lepoutre

lem is already difficult (Doumic Jauffret and Gabriel 2010) and the Floquet problem necessitates a reduction of the model to a small ODE system (3-compartments). 1.2 The question of transition rates The main medical reason to focus on transition rates is the use of cytostatic drugs in cancer therapy. Such drugs are slowing the cell cycle, their effect may be represented by a reduction of a coefficient K i→i+1 in the model. One of the main difficulties here comes from the fact that such coefficients appear twice in the equations, both as a death rate and as a birth rate. It was already noticed by Clairambault et al. (2006) that the effect of averaging the transition coefficients did not seem obvious and it was established in Clairambault et al. (2009) that no general comparison involving the transition rates can be expected. Indeed when each transition rate K i→i+1 (t, x) is T replaced by its time average T −1 0 K i→i+1 (t, x)dt, one can obtain either a faster of slower growth. 1.3 A simplified division model As we want to obtain a better understanding of the contribution of the transition and division rates, we focus on the simplest version of (1) 

∂t n(t, x) + ∂x n(t, x) + K (t, x)n(t, x) = 0, ∞ n(t, 0) = 2 0 K (t, x)n(t, x)d x.

(5)

The coefficient K is then interpreted as a division rate. In this simplified model, we do not take into account death rates. Indeed, time periodic age independent death rates or time independent death rates could be incorporated without changing the main features of our results. However, death rates depending on both variables (age and time) could drive to specific synchronizations between transition and death rates which are out of the scope of this paper. We are interested in the properties of the principal eigenvalue λ F which determines the growth of the population. It is determined as the sole eigenvalue associated to nonnegative eigenvectors: the direct eigenvector ⎧ ∂N ∂N ⎪ ⎪ ⎪ ∂t + ∂ x + (λ F + K (t, x))N (t, x) = 0, ⎨  ∞ n(t, 0) = 2 K (t, x)N (t, x)d x, ⎪ ⎪ ⎪ 0 ⎩ N (t + T, x) = N (t, x), N  0, N ≡ 0.

(6)

and the adjoint eigenvector 

123

∂φ ∂φ − + (λ F + K (t, x))φ(t, x) = 2K (t, x)φ(t, 0), ∂t ∂x φ(t + T, x) = φ(t, x), φ > 0. −

(7)

Discrete limit and monotonicity properties …

For uniqueness issues, these equations are usually completed with the following normalization: 1 T

 0

T



∞ 0

N d xdt =

1 T



T 0



∞

N (t, x)φ(t, x)d xdt = 1.

(8)

0

Various properties of the Floquet eigenvalue have been studied in Bacaer and Abdurahman (2008), Clairambault et al. (2006, 2007a, 2009, 2011) with the comparison with time independent systems (in which K (t, x) is replaced by suitable time averages) in mind. Here, we are interested in properties of the eigenvalue (and more generally of the dynamics) that are intrinsically related to the fact that the coefficients are timedependent. We consider in Sect. 3 a division rate of the form K (t, x) = κψ(t)1[a,+∞[ (x).

(9)

Here, κ > 0 is an amplitude parameter and the function ψ(t)  0 expresses the dependence of the division rate as a function of time. We denote by 1 the indicator function of a set. The term 1[a,+∞[ (x) indicates that the cells only divide when x  a, i.e., when they are older than a “maturity age” a, which we prefer here to refer to as the maturation age, to avoid the confusion with maturity structured models. Existence results for the first eigenvalue of renewal equations with periodic coefficients are often established with very strong assumptions on the coefficients (in order to ensure λ > 0), see for instance Theorem 5.1 in Michel et al. (2005). Since the assumptions (especially Assumption (5.2)) are generally not satisfied in our setting, we establish in the appendix an existence result for the Floquet eigenvalue for division models with separated variables.

1.4 Main results In this section, we state our main results. Proof will be given in Sects. 2 and 3. Our goal is to explore the behavior of the growth rate when the amplitude parameter κ in (9) tends to ∞. In this case, the ’survival’ probability (to be understood as the probability x to reach a certain age without undergoing division) P(t + x, x) = e− 0 K (t+s,s)ds for an individual born at time t converges either to 0 or 1. This means that at time t and age x, either none or all cells have divided. The growth rate is represented by the Floquet eigenvalue, λ F , which is the unique scalar λ solution of (6), (7). Hence, we shall investigate the asymptotic behavior of λ F as a function of κ. To this end, we first study the monotonicity of the Floquet eigenvalue with respect to the division coefficient K (t, x). Since every division replaces one old cell by two new ones, it is natural to think that the growth rate is a nondecreasing function of K (t, x). Our first result is an explicit counter example, showing that this apparently intuitive property is not valid due to the time and age dependence of the division rate.

123

S. Gaubert, T. Lepoutre

Proposition 1 There exist division rates such that ∀(t, x), 0  K 1 (t, x)  K 2 (t, x), and the corresponding Floquet eigenvalues satisfy λ1F > λ2F . This means that for these division rates K 1 , K 2 , the lowest division rate leads to the fastest growth! However, the following theorem identifies a subclass of transition rates for which the Floquet eigenvalue is indeed a monotone function of the transition rate. Theorem 1.1 (Comparison principle) If K 1  0 satisfies  v →

v

t

K 1 (s, s − v)ds is nondecreasing for any t,

and λ1F , λ2F  0, then for all K 2  0, (∀(t, x) K 2 (t, x)  K 1 (t, x)) ⇒ λ2F  λ1F , (∀(t, x) K 2 (t, x)  K 1 (t, x)) ⇒ λ2F  λ1F . The proof relies on a disaggregation idea of independent interest (related ideas have been seen in Perthame and Touaoula (2008) for a time homogeneous size structured model). The condition on K 1 has a natural interpretation. Indeed the quantity e−

t v

K 1 (s,s−v)ds

represents the probability for an individual born at time v of having aged up to time t without dividing. The condition then ensures that this probability increases with v. This means that individuals that were born a longer time ago are more likely to have divided by time t. When K has the form (9) with inf ψ > 0, the intuition predicts that the growth rate goes to the limit log 2/a as κ tends to infinity, because every cell will divide shortly after reaching age a when the parameter κ is large. Our next result confirms that this is the case. Proposition 2 Suppose that ψ is a bounded T -periodic function with inf ψ > 0. Let λκF (a) be the principal Floquet eigenvalue associated to the division rate given by (9). Then, for all κ > 0, 0  λκF (a) 

123

log(2) , a

(10)

Discrete limit and monotonicity properties …

and lim λκF (a) =

κ→∞

log(2) . a

(11)

When ψ vanishes, meaning that the division of cells is blocked at certain times of the day, we shall see that complex synchronization phenomena appear. We shall assume that ψ has a square wave shape, i.e., that ψ is a T -periodic function such that, for some 0 < τ < T ,  1 for t ∈ [0, τ ) ψ(t) = . (12) 0 for t ∈ [τ, T ) The following is the main result of this paper. It shows that, as the scaling factor κ of the division rate tends to infinity, the Floquet eigenvalue converges to an explicit staircase shape function, which will be shown in Sect. 3.2.1 to have a simple interpretation in terms of synchronization of reproduction phenomena with the circadian clock. We denote by x the smallest integer which is not less than a real number x. Theorem 1.2 (Discrete limit) Suppose that the division rate is given by K (t, x) = κψ(t)1[a,∞) (x), where ψ is given by (12) and let λκ (a) be the associated principal Floquet eigenvalue. Denote Wτ = [τ, T ] + NT and let Na be the integer defined by Na a = min Wτ ∩ Na. Then, the limit of the Floquet eigenvalue λκ (a) as κ → +∞ is given by λ∞ (a) =

Na log 2. Na a/T T

(13)

Moreover, we have the estimate

λ∞ (a)  λκ (a)  λ∞ (a) 1 − e−κr (a) ,

(14)

where the function r (a) (the convergence rate) satisfies

ar − a , τ , ar = sup{a  , λ∞ (a  ) = λ∞ (a)}, r (a)  min 2

(15)

and r is positive outside a closed countable set. The function a → λ∞ (a) is a staircase with a countable number of steps (see Fig. 1 below and Figs. 8 and 9 in Sect. 3.5). The quantity ar − a which controls the rate represents the distance to the right end of the step containing a (on Fig. 1 for instance, we can observe that the speed of convergence depends on this distance).

123

S. Gaubert, T. Lepoutre

4log2

3log2

2log2

log2

0.0 0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9 1.0 1.1 1.2 minimal age for division a

Fig. 1 Convergence for τ = 2/3 and T = 1. Limit in dashed blue and κ = 10, 30, 50 (in black, green and red respectively) (color figure online)

There is a special case in which the limit of the Floquet eigenvalue can be rewritten in terms of integer parts. If τ/T < 1/2, and 0 < a < T , then, we have λ∞ (a) =

τ a

log 2.

We shall derive intuitively Formula (13) in Sect. 3.2.1, by introducing a discrete dynamical system, equipped with a multiplicative functional, the asymptotic geometric mean of which yields λ∞ (a). It allows one to interpret Na as the number of cell cycles that may occur consecutively before the periodic control blocks the division. The existence of such a simple asymptotic formula reflects the fact that the PerronFrobenius operator describing the evolution of the population degenerates as κ → ∞. The rate r (a) appears naturally in the proof of convergence. This proof relies on Theorem 1.1 (the comparison principle), on the generational model already used in the proof of this theorem, and on explicit bounds of the density, for a special initial condition leading to a trajectory “close” to that of the discrete dynamical system. The nature of this degeneration may be best understood by considering an elementary analogy with matrices. Let us call a n × n nonnegative matrix M unambiguous if, for every i, j ∈ {1, . . . , n} and every k  1, in the expansion (M k )i j =

Mii1 Mi2 i3 . . . Mik−1 j

1i 1 ,...i k−1 n

there is at most one non-zero term. Then, it is readily seen that the Perron root of M, which coincides with limk M k 1/k for any norm  · , is nothing but the maximal geometric mean (Mi1 i2 . . . Mik i1 )1/k where the maximum is taken over all k  1 and all sequences i 1 , . . . , i k of elements of {1, . . . , n}. Therefore, if a nonnegative matrix Mκ converges to an unambiguous matrix M∞ , the Perron-root of Mκ (which depends continuously of the entries of Mκ ) must converge to the maximal geometric mean of a

123

Discrete limit and monotonicity properties …

circuit of the digraph of M∞ . Theorem 1.2 is an infinite dimensional analogue of this phenomenon, in which a Perron-Frobenius operator degenerates to a “nonambiguous” infinite dimensional operator. The periodic orbits of the discrete dynamical system that we construct play the role of the “circuits” in the finite dimensional case. The convergence rate in Theorem 1.2 is reminiscent of large-deviation or Laplace asymptotics. However, this theorem differs in its essence of already known asymptotics results in Perron-Frobenius theory with such large deviations flavor, including zerotemperature asymptotics of Ising type models or Frenkel-Kantorova models (Akian et al. 1998; Thieullen and Garibaldi 2012). In these problems, the free energy per site, which is the limit of the Perron eigenvalue, in a log-scale, turns out to be the ergodic constant of a “one player deterministic game”, i.e., a controlled deterministic dynamical system. In the present case, the dynamical system which determines the limit involves no control, so that Theorem 1.2 may be thought of as a large deviation result with a 0-player limit. This is due to the previous mentioned unambiguity phenomenon. When ψ vanishes, we may consider ψ ε (t) := ψ(t) + ε, with ε > 0, and define the growth rate λ( , κ) with an obvious notation. A comparison of the previous theorems shows that the limits in κ and ε do not commute log(2) = lim lim λ( , κ) = lim lim λ( , κ) = λ∞ (a). κ→∞ ε→0 ε→0 κ→∞ a

1.5 Outline of the paper The paper is organized as follows. In Sect. 2 we discuss the issues of the monotonicity of the Floquet eigenvalue λ F with respect to the division rate K . We prove here Proposition 1 and Theorem 1.1. In Sect. 3 we prove the main result of the paper and focus on the case K (t, x) = κψ(t)1[a,∞) (x). We start by proving the simple case when ψ is lower bounded and then establish the proof of Theorem 1.2, case where ψ is periodic square wave. First we give some intuitions on the construction of λ∞ in Sect. 3.2.1. Then we establish that λ∞ is an greater than λκF (a). Equation (14) which is the core of Theorem 1.2 is established in Sect. 3.4 for the case r (a) > 0. The case r (a) = 0 is treated using a continuity argument established in Sect. 3.2.2. Numerical illustrations of the convergence are given in Sect. 3.5. We end up the paper with some concluding remarks on the possible generality of our results. The appendix is dedicated to the proof of an existence theorem for λ F which is needed for this article and was not available in the literature.

2 Monotonicity with respect to the division rate We first focus on the question of the monotonicity of the Floquet eigenvalue with respect to the division rate.

123

S. Gaubert, T. Lepoutre

2.1 A counter example in which the growth rate is not a monotone function of the division rate We next show that increasing the division rate may decrease the growth rate of the population. The counter example is based on the following exactly solvable model. We assume that the period is T = 1, and fix α ∈ (0, 1). It will be convenient to write I1 = [0, α) + Z,

I2 = [α, 1) + Z,

and χ j (t) := 1 I j (t),

j = 1, 2, so that χ1 + χ2 = 1.

We shall assume that the division rate is of the form K (t, x) = χ1 (t − x)K 1 (t) + χ2 (t − x)K 2 (t). Integrating along a characteristic, we get  n(t, x) = n(t − x, 0) exp −

x

K (t − x + y, y)dy

0

and so  n(t, x) = n(t − x, 0) exp −

x

K j (t − x + y)dy , if t − x ∈ I j .

0

In other words, when an individual is born at time t ∈ I j , then its next division will occur with a rate K j (t). This model has an intuitive interpretation. We may assume for instance that I1 represents the days, whereas I2 represents the nights. If being born during the night or during the day influences the fertility (division rate) of the individuals, we arrive at this model. We may then consider that there are actually two subpopulations n j (t, x), j = 1, 2, corresponding respectively to individuals born during the day or during the night. We have n j (t, x) = n(t, x)χ j (t − x),

j = 1, 2

They satisfy the equations: ⎧ ⎨ ∂t n j + ∂x n j + K j(t)n j = 0, j = 1, 2 ∞ (K 1 (t)n 1 (t, x) + K 2 (t)n 2 (t, x))d x ⎩ n j (t, 0) = 2χ j (t) 0

Denoting by 

∞

P j (t) = 0

123

n j (t, x)d x

j = 1, 2.

Discrete limit and monotonicity properties …

the total population of type j at time t, we get that the vector P(t) := (P1 (t), P2 (t)) satisfies the ODE: d P(t) = M(t)P(t) dt where M(t) :=

2χ1 (t)K 2 (t) (2χ1 (t) − 1)K 1 (t) (2χ2 (t) − 1)K 2 (t) 2χ2 (t)K 2 (t)



We now assume that the functions K 1 (t) and K 2 (t) are constant during the day and during the night, meaning that K 1 (t) = a1 χ1 (t) + b1 χ2 (t),

K 2 (t) = a2 χ1 (t) + b2 χ2 (t),

where the coefficients ai , bi are constant. With this form of the coefficients, we obtain M(t) = χ1 (t)Ma + χ2 (t)Mb where Ma :=

a1 2a2 0 −a2



,

Mb :=

−b1 0 2b1 b2

.

An important point here, still with the above interpretation is the following: a2 and b1 could be considered in some sense as transition coefficients between populations n 1 and n 2 . Since the linear dynamics M(t) switches between the constant dynamics Ma and Mb , which are successively exercised during the intervals [0, α) and [α, 1), we get that     P(1) = exp (1 − α)Mb exp α Ma P(0). Denoting by ρ the Perron root (spectral radius) of a nonnegative matrix, we arrive at the following expression for the Floquet eigenvalue

    λ F = log ρ exp (1 − α)Mb exp α Ma .

(16)

Let us now consider the situation in which a2 = b2 = 0, meaning that an individual born during the night is sterile (division will never occur). Then, the second columns of the matrices Ma and Mb vanish, and so     exp (1 − α)Mb exp α Ma =



exp(αa1 − (1 − α)b1 )

0 1



123

S. Gaubert, T. Lepoutre

Fig. 2 The counter example. The Floquet eigenvalue as a function of the division rates b1 and b2 , when a1 = 10 and a2 = 0.1. When b2 is small, increasing b1 (and so, increasing the division rate K (t, x)) decreases the growth rate

(the value of the off diagonal entry, denoted by , is irrelevant). Since the spectral radius of a nonnegative triangular matrix is the maximum of its of its diagonal entries, we get   λ F = max αa1 − (1 − α)b1 , 0 . Hence, λ F is an increasing function of the division rate a1 , but a decreasing function   of the division rate b1 . This establishes Proposition 1. It should be noted that the preceding counter-example subsists when a2 and b2 are sufficiently small, because the spectral radius is a continuous function of the parameters. This is confirmed by Fig. 2, in which the Floquet eigenvalue is plotted as a function of b1 and b2 for a1 = 10 and a2 = 0.1. Remark 1 This counter example can be understood intuitively in Fig. 3. The values of the division rate K (t, x) are represented (Left, the values are constant on each cell). Consider now a chain of descendants of the same individual (i.e., this individual, one of his children, one of the children of this children, etc.). This can be represented by a Markov process X t with state space in R+ (the age), which jumps from age x to age 0 with rate K (t, x). Each time the path hits the axis x = 0 corresponds to a division. If the path hits x = 0 during the night, since b2 = 0, the corresponding individual is sterile (no division will ever occur, which is reflected that the path is now an unending line of slope 1 in the (t, x) plane). The production of sterile individuals occurs when the process jumps from a cell labeled “b1 ” to the x axis, and this occurs precisely with

123

Discrete limit and monotonicity properties …

x

x

b1 b2

a1 b1 a1 α

b2

a2 1

t

t

Fig. 3 The counter example explained. The values of K (t, x) (left). The night is shown in grey. A typical path (right, broken line). The child born during the night is sterile (final segment of the path), and the division coefficient b1 determines the creation of this child

rate b1 , which explains intuitively why the growth rate, λ F , is a decreasing function of b1 . Remark 2 The analysis of this system does not rely on the result of the appendix for the existence of the Floquet eigenvalue, since here the exponential growth is determined explicitly. Indeed, the preceding study merely determines the growth rate of the “aggregated” population P(t), which is determined as the Perron eigenvalue of a matrix. However, the growth of n(t, x) is readily derived from the one of P(t). To see this, it is convenient to introduce the adjoint positive Floquet eigenvector of M(t), whose existence and uniqueness are guaranteed by Perron-Frobenius theory (as soon as the coefficients a1 , a2 , b1 , b2 are positive): −

     d  φ1 φ2 + λ φ1 φ2 + φ1 φ2 M(t) = 0. dt

We have d (φ(t), P(t)) = λ(φ(t), P(t)) dt which can be read as, denoting φ(t, x) = φ1 (t)χ1 (t − x) + φ2 (t)χ2 (t − x), d dt



∞ 0

 n(t, x)φ(t, x)d x = λ

∞

n(t, x)φ(t, x)d x, 0

It is straightforward to check that φ1 and φ2 are positive (still assuming that the coefficients a1 , a2 , b1 , b2 are positive). Hence, there exists M  m > 0 such that m  φ(t, x)  M, and therefore,    λt 0 λt me n d x  n(t, x)d x  Me n 0 d x.

123

S. Gaubert, T. Lepoutre

x t v

v

t u

K(s, s − u)ds

s

t

u

K(s, s − v)ds

Fig. 4 The sufficient monotonicity condition: the survival probability at time t (the exponential of the opposite of the integral) increases with the birth time v

This shows that the Floquet eigenvalue λ computed from the aggregated dynamical system does determine the behavior of the solutions. 2.2 A sufficient condition for monotonicity The preceding counter example is due to a “birth day penalty” (the date of birth of the individual influences critically its division rate). We shall avoid such a pathological behavior by making the following assumption:  (t, v) →

t v

K (s, s − v)ds is nonincreasing in v.

(17)

This condition is fulfilled when the division rate K (t, x)  0 depends only on one of the two variables t and x. It is also fulfilled if for any t, K (t, ·) is a nondecreasing function (and particularly, in the case of separated variables K (t, x) = ψ(t)B(x) when B is nondecreasing). This is illustrated in Fig. 4. The condition requires the integral of K over the characteristic shown on the figure to be nonincreasing in v. Note that this has a probabilistic interpretation in terms of the Markov process X t introduced in the preceding section. Indeed,  t

exp − K (s, s − v)ds v

represents the ‘survival’ probability at time t of an individual born at time v, i.e., the probability, conditional to X v = 0, that the process makes no jump to point 0 before time t. The condition requires this probability to be a nondecreasing function of v, in other words, that by delaying the birth, the survival probability at a given time cannot decrease. Theorem 2.1 Provided K 1 or K 2 satisfies (17) and K 1  K 2 , then, for any n 0 ∈ L 1 (R+ ), if we denote by n i the solution of (5) with K = K i , we have, for any t  0,  0

123

∞

 n (t, x)d x  1

0

∞

n 2 (t, x)d x.

Discrete limit and monotonicity properties …

To establish this result, we introduce the following system of PDEs: ⎧ ∂t n i + ∂x n i + K (t, x)n i = 0, i  0, ⎪ ⎪ ⎪ ⎪ ⎨ n 0 (t, 0) = 0 ∞

(18)

K (t, x)n i−1 (t, x)d x, i  0 n i (t, 0) = ⎪ ⎪ ⎪ 0 ⎪ ⎩ n 0 (t = 0, x) = n 0 (x), n i (t = 0, x) = 0, for i  1.

 We easily check that n = i 0 2i n i is the solution to (5). Intuitively, the term 2i n i represents the number of individuals of the “generation” i, and so, we shall refer to this model as the “generational model” in the sequel. We denote by (n i1 )i 0 , (resp. (n i2 )i 0 ) the solution to (18) with K = K 1 (resp. K = K 2 ). We also introduce: S j (t) :=

 ij

∞

n i (t, x)d x, for j  0.

0

A short computation leads to 

∞

n(t, x)d x = S0 (t) +

0

2 j−1 S j (t).

j 1

Lemma 2.2 We have:  ∞ S0 (t) = S0 (0) = n 0 (x)d x, 0   ∞ 0 S1 (t) = n (x) 1 − exp − 0

(19) t



K (s, x + s)ds d x,

(20)

0

and for j  1,  S j (t) =

t

 n j−1 (s, 0) 1 − exp −

s=0

t−s



K (s + y, y)dy

ds.

(21)

0

Proof It is straightforward to show that for every j  0, d S j (t) = n j (t, 0), dt and so 

t

S j (t) =

n j (s, 0)ds + S j (0).

(22)

0

This leads immediately to the first statement, since for all s > 0, n 0 (s, 0) = 0.

123

S. Gaubert, T. Lepoutre

To establish the next statements, we shall use the characteristics: n 0 (t, x) = 0, if x  t,

 n 0 (t, x) = n (x − t) exp −

t

0

K (s, x − t + s)ds , if x  t,

(23)

0

and for i  1, n i (t, x) = 0, if t  x,

 n i (t, x) = n i (t − x, 0) exp −



x

K (t − x + y, y)dy , if x  t.

(24)

0

The lack of symmetry between the cases i  1 and i = 0 is due to the boundary conditions: n 0 (0, x) = n 0 (x) is known, whereas for i  0, n i (t, 0) can be inductively assumed to be known, using n i−1 . Since n 1 (0, x) = 0, we have S1 (0) = 0. Using (22), we get  S1 (t) =



t

t



∞

n 1 (s, 0)ds = K (s, x)n 0 (s, x)d xds s=0 s=0 x=0  s

 t  ∞ 0 K (s, x)n (x − s) exp − K (u, x − s + u)du d xds = s=0 x=s 0  s

 t  ∞ 0 = K (s, y + s)n (y) exp − K (u, y + u)du dyds, s=0

y=0

0

and, interchanging the order of integration,  S1 (t) = =



∞

y=0  ∞

n (y) 0

t

 K (s, y + s) exp −

s=0

K (u, y + u)du dsdy

0

 n 0 (y) 1 − exp −

0

s

t



K (s, y + s)ds

dy,

0

which is the second statement. Assume now that j  2. Arguing as above, but using this time (24) instead of (23), we get  S j (t) =

t



∞

K (s, x)n j−1 (s, x)d xds

s=0 x=0  t  s

 x

K (s, x)n j−1 (s − x, 0) exp − K (s − x + y, y)dy d xds s=0 x=0 0  x

 t  t = K (s, x)n j−1 (s − x, 0) exp − K (s − x + y, y)dy dsd x x=0 s=x 0  x

 t  t−x K (u + x, x)n j−1 (u, 0) exp − K (u + y, y)dy dud x = =

x=0 u=0

123

0

Discrete limit and monotonicity properties …

 =

t

 n j−1 (u, 0)

u=0

t−u

 K (u + x, x) exp −

x=0

x

K (u + y, y)dy d xdu,

0

 

which leads to the third statement.

Theorem 2.1 will be obtained as an immediate consequence of the following lemma. Lemma 2.3 Suppose K 1 satisfies (17), then K 2  K 1 ⇒ S 2j (t)  S 1j (t), ∀t, j K 2  K 1 ⇒ S 2j (t)  S 1j (t), ∀t, j. Proof The proof is performed by induction on j. We shall only consider the case in which K 2  K 1 , the other case being similar. For j = 0, 1, the conclusion follows readily from the two first statements in Lemma 2.2. Assume now that j  2. We set  Q (t, v) := 1 − exp −

t−v

k

K (v + y, y)dy k

0

Using (21), we get  t  t S 1j (t) − S 2j (t) = n 1j−1 (v, 0)Q 1 (t, v)dv − n 2j−1 (v, 0)Q 2 (t, v)dv, 0 0  t  t (n 1j−1 (v, 0) − n 2j−1 (v, 0))Q 1 (t, v)dv + n 2j−1 (v, 0)(Q 1 (t, v) − Q 2 (t, v))dv = 0 0  t t  d (S 1j−1 (v) − S 2j−1 (v)) Q 1 (t, v) = (S 1j−1 (v) − S 2j−1 (v))Q 1 (t, v) 0 − dv 0  t 2 1 2 n j−1 (v, 0)(Q (t, v) − Q (t, v))dv, + 0  t  t d 1 1 2 (S j−1 (v) − S j−1 (v)) Q (t, v) + n 2j−1 (v, 0)(Q 1 (t, v) − Q 2 (t, v))dv. =− dv 0 0 d The condition (17) ensures the nonpositivity of dv Q 1 (t, v). Moreover, K 2  K 1 readily implies Q 2  Q 1 . Using the induction hypothesis, we deduce that S 2j  S 1j .

We can now give the proof of Theorem 1.1 as a consequence of Theorem 2.1. We assume that K i satisfies λiF  0. In this case, we can establish the boundedness of the adjoint eigenvector φi . Indeed using the characteristic formula we obtain, as equation (5.8) in Michel et al. (2005) φi (t + x, x)  ∞ s i =2 K i (t + s, s)e− x K i (t+u,u)du−λ F (s−x) φi (t + s, 0)ds  2φ(., 0)∞ . x

123

S. Gaubert, T. Lepoutre

Therefore, since for all solution of Eq. (5) with coefficient K i satisfy 

∞

n i (t, x)φi (t, x)d x = eλ F t i

0



∞

n i (t = 0, x)φ(0, x)d x,

0

we have in particular 

∞



∞

n i (t, x)d x  c

0

n i (t, x)φi (t, x)d x  ceλ F t . i

0

Now we assume by contradiction for instance that K 1  K 2 and λ1F < λ2F . Taking as initial condition n 1 (t = 0, x) = n 2 (t = 0, x) = N1 (t = 0, x), then we have by Theorem 3 for t > 0 

∞

n 1 (t, x)d x   ∞ 1 N1 (t, x)d x  = eλ F t

0

0

∞



∞

n 2 (t, x)d x  c

0

n 2 (t, x)φ 2 (t, x)d x = ceλ F t . 2

0

In particular, since N1 is T −periodic, this leads to an inequality of the form eλ F nT  Ceλ F nT 1

2

for all n > 0 contradicting the fact λ1F < λ2F . The alternative case K 2  K 1 is proved by the same argument. This ends the proof of Theorem 1.1.   Remark 3 Lemma 2.3 has an interpretation in terms of stochastic order or majorization (Marshall and Olkin 1979). Recall that for R-valued random variables Z 1 , Z 2 , the stochastic order st is such that Z 2 st Z 1 holds if and only if P(Z 2  t)  P(Z 1  t) for all t ∈ R. Now, let X (t) denote the Markov process representing the age of a single individual, defined in Remark 1, and let Y (t) denotes  its number  ∞ of jumps at time t, i.e., the “generation” of an individual. Then, S j (t) = i  j 0 n i (x, t)dt = P(Y (t)  j), and so, denoting by Y 1 (t), Y 2 (t) the two processes obtained in this way from K 1 , K 2 , we see that the properties of Lemma 2.3 are equivalent to K 2  K 1 ⇒ Y 2 (t) st Y 1 (t), ∀t K 2  K 1 ⇒ Y 2 (t) st Y 1 (t), ∀t.

3 Asymptotics of the growth rate We now assume that the division rate is of the form K (t, x) = κψ(t)1[a,+∞[ (x)

123

Discrete limit and monotonicity properties …

an determine the limit of the Floquet eigenvalue of System (5) as the division rate tends to infinity, i.e., as the scaling parameter κ tends to infinity. 3.1 When ψ is positive: Proof of Proposition 2 The Floquet eigenvalue λκF satisfies ⎧ ∂   N κ (t, x) + ∂∂x N κ (t, x) + λκF + κψ(t)1[a,+∞[ (x) N κ (t, x) = 0, ⎪ ⎪ ∂t ⎨ ∞ N κ (t, 0) = 2κψ(t) a N κ (t, x)d x, ⎪ ⎪ ⎩ κ N > 0, T -periodic.

(25)

The existence of the Floquet eigenvector N κ is derived in the appendix as a consequence of the Krein-Rutman theorem. We normalize the Floquet eigenvector by requiring that  T ∞ N κ (t, x)d xdt = 1. (26) 0

0

To establish Proposition 2, we shall think of the Floquet eigenvalue λκF = λκF (ψ) as a function of the time modulation ψ. Since the division rate is a nondecreasing function of ψ, using the comparison principle (Theorem 1.1), we get λκF (ψ)  λκF (ψ)  λκF (ψ),

(27)

where ψ and ψ denote the constant functions equal to the minimum or maximum of ψ, respectively. When the function ψ(t) is equal to a constant α, the Floquet eigenvector can be chosen to be independent of t, so that N κ (t, x) = N κ (x). Then, the Floquet eigenvalue λκF (α) is actually an ordinary eigenvalue. Indeed, it is the only scalar μ solution of the system ⎧ d   N κ (x) + μ + κα1[a,+∞[ (x) N κ (x) = 0, ⎪ ⎪ d x ⎨ ∞ N κ (0) = 2κα a N κ (x)d x, ⎪ ⎪ ⎩ κ N >0

(28)

Integrating the latter differential equation, we get   N κ (x) = N κ (0) exp − μx − κα(x − a)+ . Using the boundary condition, and eliminating N κ (0), we determine μ via the following implicit equation as in Clairambault et al. (2009): 2καe−μa = 1. κα + μ

123

S. Gaubert, T. Lepoutre

Fig. 5 Convergence of the Floquet eigenvalue to log 2/a. Limit in dashed blue and κ = 15, 20, 50 (in black, green and red respectively) (color figure online)

This gives immediately 0 < μ < log 2/a and thereby the bounds 0 < λκF (ψ)  λκF (ψ)  λκF (ψ) < log 2/a. So the first part of Proposition 2 is proved. To prove the limit, it suffices to notice that solution of 2καe−μa =1 κα + μ have a limit log 2/a when κ → +∞ as long as α > 0 is fixed. Then taking α = ψ, ψ we obtain λκF (ψ), λκF (ψ) → log 2/a as κ → +∞, Together with Eq. (27), this completes the proof of the proposition 2.   We give a numerical illustration of this theorem in Fig. 5. These numerical simulations were produced along the lines of Clairambault et al. (2009) using a monotone finite difference scheme. The Floquet eigenvalue was computed by applying the power algorithm to the one day discretized evolution operator. 3.2 When ψ can vanish: Proof of Theorem 1.2 When ψ can vanish, the conclusion of Proposition 2 fails; its proof uses the assumption that min(ψ) > 0. Then, we next show that more complex behaviors appear. In particular, the limit of λkF can be a discontinuous, staircase like, function of the maturation age a, and we shall see that the limits κ → +∞ and ψ → 0 (locally) may not commute.

123

Discrete limit and monotonicity properties …

3.2.1 Intuitive derivation of the formula via a deterministic jump Markov process We consider the original PDE (5), taking T = 1, 

∂ ∂ ∂t n(t, x) + ∂ x n(t, x) + κψ(t)1[a,∞[ (x)n(t, x) ∞ n(t, 0) = 2κψ(t) a n(t, x)d x.

= 0,

(29)

where ψ is a 1-periodic square wave, such that if t is the integer part of t, ∃ 0 < τ < 1, ψ = t → 1[0,τ [ (t − t).

(30)

We next derive the limit of the growth rate λ when κ → +∞ (within this section, λ denotes the Floquet eigenvalue λκF , thought of as a function of a and τ ). This derivation will be purely formal for the moment, leading to an ansatz the validity of which will be proved in the next sections. Intuitively, whenever the cell is old enough in the cycle (x > a) and the time is favorable (0  t − t  τ ), a division may occur. The idea is to consider that at time t = 0, we have only cells of age 0 (n(0, ·) = δ0 ), and to look formally for an asymptotic solution n(t, ·) ∼ 2m(t) δα(t) , κ → ∞ of the PDE (29), meaning that the population at time t consists mostly of 2m(t) individuals all of age α(t), as κ → ∞. If for some integer time K , we have α(K ) = 0, so that n(K , ·) ∼ 2m(K ) δ0 , using the time-periodicity and linearity of (29) we may expect that n( pK , ·) ∼ 2 pm(K ) δ0 , for all integers p  1. Commuting the limits as κ → ∞ and p → ∞, we arrive at the following expression for the limit of the growth rate as κ → ∞, λ∞ (a) :=

m(t) m(K ) = lim . t→∞ K log 2 t log 2

(31)

The idea is now that the trajectory (α(t))t 0 , representing the age of the population, is produced by a degenerate (deterministic) time-dependent càdlág jump process, in which the state space is R+ , and a particle moves to the right with unit speed, until it reaches an age and a time at which division is permitted, meaning that α(t)  a and t ∈ [0, τ [+N. At every time satisfying this condition, it jumps to point 0 (age is reset to zero and a division occurs). The number m(t) counts the number of jumps up to time t and represents the number of divisions of the cell which occurred up to that time. It suffices for our argument to consider the case in which α(0) = 0. Then, the division times can be readily computed: if we start from age 0 at time 0, then, a division should occur at every instant a, 2a, . . . , ka as long as it is permitted, meaning that a, 2a, . . . , ka ∈ [0, τ [+N. Let K a := inf{k ∈ N; ka ∈ [τ, 1[+N},

(32)

123

S. Gaubert, T. Lepoutre

8δ0

4δ

4δ0

2δ 0 δ0

4 δ 1−2a

2δa

δa 1

τ

0

Fig. 6 Example of the derivation of the formula. Here 2a < τ < 3a < 1 and we obtain n(1, ·) = 8n(0, ·). Hence we expect λ∞ (a) = 3 log(2)

so that K a a is the first time at which a division is not permitted, or K a = ∞. When, for the first time, we reach a possible division time K a a that is not permitted, then the corresponding division occurs later, namely at the next period: K a a , where · denotes the upper integer part. Then, at time K a a the population divides (for the K a th time) and we have n(K a a , ·) = 2 K a δ0 . Therefore, (31) yields ∞ (a)K

eλ

a a

= 2 K a , λ∞ (a) =

Ka log 2. K a a

(33)

When K a = ∞, we shall adopt the convention that λ∞ (a) = (log 2)/a. This is illustrated in Fig. 6. Working with a semi-open interval [τ, 1[ in (32) is consistent with the càdlág nature of the process α(t), however, it leads to the possibility that K a = ∞, which is somehow unpleasant. Hence, we will use in the sequel an equivalent formulation of λ∞ , obtained by replacing [τ, 1[ by the closed interval [τ, 1], Na := min{k ∈ N; ka ∈ Wτ }, Wτ := [τ, 1] + N. Lemma 3.1 We have Na < ∞. Moreover, λ∞ (a) =

123

Na log 2. Na a

(34)

Discrete limit and monotonicity properties …

Proof If a is irrational, the finiteness of Na follows from the density of the sequence ka modulo 1, in [0, 1[; whereas if a is rational, we have ka ∈ 1 + N for some integer k, so that Na  k is finite. It remains to show that the expressions (33) and (34) coincide. Assume first that Na a ∈ [τ, 1[+N. Then, by definition of Na , ka ∈]0, τ [+N for 1  k < Na , and so, ka ∈ [0, τ [ for 0  k < Na , which implies that K a = Na . Assume now that Na a ∈ 1+N. Then, we deduce from ka ∈]0, τ [+N, for 1  k < Na that (Na + k)a ∈]0, τ [+N, for all 1  k < Na , and so K a = ∞, which by our convention regarding (33) leads to (34). In the remaining subsections, we show that (a) limκ→∞ λκ (a) does coincide with the function λ∞ (a) (Theorem 1.2). Throughout the proof, the parameter τ will be fixed. The proof of the Theorem 1.2 relies on the following strategy: (i) prove that the ansatz λ∞ (a) (staircase formula) is an upper bound for λκ (a), for every κ > 0, (Sect. 3.3); (ii) derive technical properties of the limit and of λ∞ (a); in particular monotonicity and right continuity, (Sect. 3.2.2); (iii) assume that a is not the right end of a step of λ∞ , which is established to be true outside a discrete set by Lemma 3.6, (iv) consider the generational formulation of the division equation with initial data N κ (0, x) and show that at time Na a , almost every individual has reached generation Na and deduce the convergence for such an a. (v) conclude that the limit  of λκ as κ → ∞ is given by the staircase formula for λ∞ (a), for all a > 0, using the fact that both λ∞ and limκ→∞ λκ are right continuous functions of the parameter a which coincide outside a discrete set. 3.2.2 Preliminaries: technical properties of λ∞ and of  = limκ→∞ λκ We shall need a number of technical observations regarding the function λ∞ defined in (33). Lemma 3.2 The function a → λ∞ (a) is nonincreasing. Proof Let (α(t))t 0 denote the trajectory of the jump process defined above, starting from α(0) = 0, and let m(t) denote the number of jumps (divisions) up to time t (included). We have m(t) = sup {k ∈ N; ∃0  a1  · · ·  ak  t, ai+1 − ai  a, ai ∈ [0, τ [+N} (number of divisions since time 0) Since the constraint ai+1 − ai  a gets stronger as a increases, it is then obvious that m is a nonincreasing function of a. We have, by definition of λ∞ , Na λ∞ m(t) = = . t→+∞ t Na a log 2 lim

(35)

Since for every t, m(t) is a nonincreasing function of a, the same is true for limt m(t)/t = λ∞ / log 2.  

123

S. Gaubert, T. Lepoutre

Lemma 3.3 The limit (a) = limκ→∞ λκ (a) satisfies (a) = sup λκ (a), κ>0

and the map a → (a), is nonincreasing and right continuous. Proof Since λκ (a)  log 2/a, and since, by the comparison principle (Theorem 1.1), λκ (a) is a nondecreasing function of κ, (a) = sup λκ (a)  κ>0

log 2 . a

Furthermore, since for any fixed κ, λκ (a) is nonincreasing with a, so is (a). Finally, since for any fixed κ > 0, λκ (a) is continuous with respect to a, (a), which is a supremum of a family of lower semi-continuous functions, is lower semi-continuous. Combined with monotonicity, this leads to the convenient property: ∀a > 0,

(a) = lim inf (x) = lim (a). x→a

x→a+0

That is, (a) = limκ→∞ λκ (a) is right continuous.

 

Lemma 3.4 The function a → λ∞ (a) is right continuous. Proof By Lemma 3.1, Na a = min Wτ ∩ Na =: f (a), λ∞ (a) =

f (a) Na log 2 = log 2, Na a a f (a)

where we just denote f (a) = Na a. To show that the function λ∞ is right continuous, it suffices to prove that g(a) =

f (a) ,  f (a)

is right continuous. If g(a) = 1, this means that f (a) =  f (a) , and therefore, for ε > 0 small enough, we can have n(a + ε) ∈ [na, na + τ [ for all n < Na and Na = Na+ε , f (a + ε) = f (a) + Na ε and  f (a) =  f (a + ε) . Thereby, for ε > 0 a small enough, we have g(a + ε) = g(a) +  fN(a) ε, which leads to the right continuity of g at a. If g(a) = 1, then, we can check that lim f (a + ε) = +∞.

ε→0+

and so lim g(a + ε) = 1.

ε→0+

 

123

Discrete limit and monotonicity properties …

These observations leads to the following useful property (stating basically that λ∞ is piecewise constant). Lemma 3.5 For all a > 0, Ia = {λ∞ (x) = λ∞ (a)} is an interval with nonempty interior. Moreover, inf Ia ∈ Ia . Finally, if we denote E τ = {a > 0, a = sup Ia }, we have • if a ∈ E τ then λ∞ is continuous at a, • if a ∈ E τ then Na a is an integer. We shall see that the set E τ contains exactly the right ends of steps of λ∞ where no jump occur. An infinite number of steps accumulate at the right of every point of E τ . On Fig. 1, we can see for instance that 1/2 belongs to E 2/3 . Proof Step 1: Ia is an interval Firstly, this set is necessary an interval by monotonicity: if a  ∈ Ia , then for any convex combination of a, a  , we have λ∞ (a) = λ∞ (max(a, a  ))  λ∞ (θa + (1 − θ )a  )  λ∞ (max(a, a  )) = λ∞ (a). To prove that it is nontrivial, we use the definition of Na . By definition a, 2a, . . . (Na − 1)a ∈ Wτ . As Wτ is a closed set, its complementary is open and thus for ε > 0 small enough, ]a − ε, a + ε[∩Wτ =]2(a − ε), 2(a + ε)[∩Wτ = ](Na − 1)(a − ε), (Na − 1)(a − ε)[∩Wτ = ∅. Finally, for ε > 0 small enough, we have necessarily ]Na (a − ε), Na a] ⊂ Wτ or [Na a, Na (a + ε)[⊂ Wτ (possibly both). Thereby, for any a, there exists ε > 0 such that [a, a + ε[ or ]a − ε, a] is contained in Ia which has then a nonempty interior (contains an interval of size ε). The fact that inf Ia belongs to Ia is just a consequence of right continuity. Step 2: continuity on E τ The continuity follows from right continuity of λ∞ and the fact that in this case, λ∞ (sup Ia ) = lim x→sup Ia −0 λ∞ (x), that is λ∞ is also left continuous. Step 3: Na a = Na a , for all a ∈ E τ Assume that a satisfies Na a − 1 + τ  Na a < Na a . Then, for some ε > 0, we have Na a − 1 + τ  Na a  < Na a for any a  ∈]a, a + ε[. Thereby, we can claim by definition of Na that Na   Na ,

Na  a   Na a .

Suppose now that in addition that a ∈ E τ . By definition, as a = sup Ia , we have λ∞ (a  ) = λ∞ (a) for a  ∈]a, a + ε[. Furthermore for a  ∈]a, a + ε[, λ∞ takes its

123

S. Gaubert, T. Lepoutre

values in 

N log 2, p

N  Na ,

 p  Na a ,

therefore, for a  ∈]a, a + ε[, we have |λ∞ (a) − λ∞ (a  )| 

1 log 2 (Na a )2

and this contradicts the continuity of λ∞ at a. This ends the proof of lemma 3.5 The set E τ has a last interesting property. Lemma 3.6 The set E τ is discrete. Proof Let a ∈ E τ . Thanks to Lemma 3.5, we know that λ∞ is constant on Ia and therefore a  ∈ E τ for a  ∈] inf Ia , a[. It remains to prove the existence of ε > 0 such that ]a, a + ε[∩E τ = ∅. To prove this, let m ∈ N\{0} satisfy τ  1 − m1 . Choose ε = 2N1a m . Let a  = a + ε ∈]a, a + ε[. We have 0 < ε < ε. For every integer k, we have k Na a  = k Na a + k Na ε  . As ε  then,

1 2Na ,

there exists an integer k such that k Na ε  1 < (k + 1)Na ε . We have

τ 1−

1  1 − 2Na ε < (k − 1)Na ε < k Na ε  1. m

This leads to, adding (k − 1)Na a, (k − 1)Na a + τ < (k − 1)Na a  < (k − 1)Na a + 1. Since Na a ∈ N, this means that (k − 1)Na a  ∈ int(Wτ ). If we assume that a  ∈ E τ , we have Na  a  ∈ N and then by construction 

   0, a  , . . . , Na  a  ∩ Wτ = 0, Na  a 

It follows then that a  N ∩ Wτ = (Na  a  )N, which in particular implies a  N ∩ int (Wτ ) = ∅. This contradicts the existence of p. We have proved that ] inf Ia , a + ε[∩E τ = {a}, which ends the proof of Lemma 3.6.

123

 

Discrete limit and monotonicity properties …

3.3 Upper bound on λκ We next show that λ∞ is an upper bound of λκ , for each κ > 0. This will turn out to be a relatively simple consequence of the heuristic argument, involving the earliest possible times for division, which led to the formula for λ∞ . The proof relies the generational reformulation of the division equation. Let a > 0, τ ∈]0, 1[, κ > 0 be fixed. We take the eigenvector as an initial data. ⎧ ⎪ (x)n i (t, x) = 0, ⎨∂t n i + ∂x n i + κψ(t)1  ∞[a,∞[ n i+1 (t, 0) = κψ(t) a n i (t, x)d x, ⎪ ⎩ n 0 (x) = N κ (0, x).

.

We note by Imax (t) the maximal integer i (generation) such that n i (t, ·) ≡ 0. At time t = 0, we have Imax (0) = 0. Since cells need to have an age bigger than a to change generation, we have necessarily Imax (t)  1 for t < a and more generally, Imax (t)  n, if t < na. Imax (t)  Na if t < (Na )a. We introduce the notation pa = Na a − 1. The integers Na , pa are then the minimal integers such that pa + τ  Na a  pa + 1. We use now the properties of the equation: since no change of generation may occur between times pa + τ and pa + 1, we have Imax (t)  Na if pa + τ  t < pa + 1. Now, from the definition of N κ and λ := λκ , we have ∞

2i n i (t, x) = eλt N κ (t, x).

i=0

In particular, for t < pa + 1, we have Imax

(t)

2i n i (t, x) =

Na

2i n i (t, x) = eλt N κ (t, x).

0

i=0

We integrate with respect to x and obtain eλt



∞ 0

N κ (t, x)d x =

Na

0



∞

2i 0

n i (t, x)d x  2 Na

Na 

0

∞

n i (t, x)d x.

(36)

0

123

S. Gaubert, T. Lepoutre

But we also have from the properties of the generational formulation and the definition of Imax ,  ∞  ∞ Na  ∞ ∞  ∞

n i (t, x)d x = n i (t, x)d x = n 0 (0, x)d x = N κ (0, x)d x. 0

0

0

0

0

0

(37) Combining (36–37), we get e

λt



κ



κ

N (t, x)d x  2

Na

0

∞

N κ (0, x)d x, t < pa + 1.

0

As κ is fixed, we can use the continuity of   exp λκ ( pa + 1)



∞

∞ 0

N κ (t, x)d x. We have therefore 

κ

N ( pa + 1, x)d x  2

Na

0

∞

N κ (0, x)d x.

0

As N κ ( pa + 1, x) = N κ (0, x), we have   exp λκ ( pa + 1)  2 Na . And therefore, as expected λκ (a) 

Na log 2 = λ∞ (a). pa + 1  

3.4 Bounding λ from below Step 1: Deriving inequality (14). (with ε instead of r .) We derive the bound through the following technical lemma, in which the sequence θ1 , . . . , θk represents possible successive division times of an individual of age 1 − τ at time 0. Lemma 3.7 Suppose there exist θ1 , . . . θ Na and ε  0, satisfying the following properties: ⎧ ⎨ θ1  max(a − 1 + τ, 0), ∀i  Na − 1, θi+1 − (θi + ε)  a, (38) ⎩ ∀i ∈ {1, . . . , Na } θi + ε < θi  + τ, then the following inequality holds: κ (a)( p +1) a

eλ

 (2 − e−κε ) Na .

(39)

Intuitively, this lemma shows that if there is a “tube” of positive width ε, limited from below by such a sequence θ1 , . . . , θ Na and included in the set of times at which division is permitted, then, “enough” individuals can follow this tube, so the Floquet eigenvalue λκ (a) can be bounded below by (Na log 2)/( pa + 1), up to a correction of order O(exp(−κ )) as κ → ∞.

123

Discrete limit and monotonicity properties …

age α(t) ar = 0.3 a = 0.22 −1 + τ

θ1

θ2

θ3

t

Fig. 7 Illustration of Lemma 3.7 and of Lemma 3.12: the trajectory of the jump process α(t) starting from α(−1 + τ ) = 0 is shown in bold for a = 0.22 and τ = 0.6 (division is not permitted in the light ray regions). Here, Na = 3 and pa = 0, so λ∞ (a) = 3 log 2. The minimal division age a can be increased up to ar = 0.3, leaving λ∞ (a) unchanged. The conditions of Lemma 3.7 can be satisfied with a positive ε, because, for all a  ∈ [a, ar [, the combinatorial type of the trajectory of the jump process α(t) is unchanged (trajectory shown by a thin line, for a  close to ar ); choosing θ1 := 0, the times θ2 and θ3 can be taken to be the abscissæ of the left of the dark grey regions of width ε < (ar − a)/2 = 0.04

The existence of the sequence θ1 , . . . , θ Na will be established in Step 2 below when ε = 0. We shall see in Step 3 that such a sequence does exist, for some ε > 0, provided that a does not coincide with the right end ar of a step of the staircase function λκ . This is illustrated in Fig. 7. Proof As for the upper bound, we consider the generational reformulation 18, in which K (t, x) = κψ(t)1[a,∞[ (x) and the eigenvector N κ (t = 0, x) is the initial data. From the definition, we have

  2i n i ( pa + 1, x) = N κ (0, x) exp λκ (a)( pa + 1) . i

Especially,



∞

2i

  n i ( pa + 1, x)d x = exp λκ (a)( pa + 1)



0

i

∞

N κ (0, x)d x.

0

We recall that N κ (0, x) = 0 for x ∈ [0, 1 − τ [. As we look for a bound from below, thanks to the comparison principle, it is sufficient to consider the individuals that are the youngest at the beginning, that is the individuals starting at age (1 − τ ). As K (t, x) is nondecreasing with respect to x, we know that if we choose K  (t, x)  K (t, x), and n i with obvious notations, then we have thanks to the comparison principle

 i

2

i

0

∞

n i ( pa

  + 1, x)d x  exp λκ (a)( pa + 1)



∞

N κ (0, x)d x.

(40)

0

We choose ⎧ ⎨ κ if t ∈]θ1 , θ1 + ε], x  a K  (t, x) = κ if t ∈]θi + (i − 1)ε, θi + iε], x  a, ⎩ 0 otherwise.

123

S. Gaubert, T. Lepoutre

From the definition of (θ1 , . . . θ Na ) and ε, we have K  (t, x)  K (t, x). The most   important property of the transition rate K  is the following. Lemma 3.8 Suppose that n i are defined by the generational dynamics with transition rate K  , then for any 1  i  Na , supp n i (θi , ·) ⊂ [a, +∞[. Proof We have for x < a, θi −x ∈ 0. Therefore for any j  0 n j+1 (θi − x, 0) =





∞

0

i ]θi +(i −1)ε, θi +iε] and thereby K

 (θ −x, ·) i

=

K  (θi − x, y)n j (θi − x, y)dy = 0,

and using the characteristics n j+1 (θi , x) = n j+1 (θi − x, 0)e−

x 0

K  (θi −x+s,s)ds

= 0.

Finally, as supp n 00 ⊂ [1 − τ, +∞[ and as supp n 0 (t, ·) = supp n 00 + t, we have since θi  a − 1 + τ for any i supp n 0 (θi , ·) ⊂ [a, +∞[.   The latter inclusion means that for any j > 0, at time θ j , every individual is mature enough to divide. Lemma 3.9 The population n i satisfies the following equality: defining Ii (t) = ∞ n (t, x)d x, we have i 0  ∀ j > 0, ∀i,

I0 (θ j + ε) = e−κε I0 (θ j ), Ii+1 (θ j + ε) = e−κε Ii+1 (θ j ) + (1 − e−κε )Ii (θ j ).

Proof It is essentially based on the previous lemma. By construction, we have ε < a. Since we already proved that supp n i (θ j , ·) ⊂ [a, +∞[, we can claim supp n i (θ j +ε, ·) ⊂ [0, ε] ∪ [a + ε, +∞[ for i  1, supp n 0 (θ j +ε) ⊂ [a + ε, +∞[ From the definition of the n j , we have, ∀i, ∀x  a + ε, n i (θ j + ε, x) = n i (θ j , x − ε)e−

ε 0

K  (θ j +s,x−ε+s)ds

= n i (θ j , x − ε)e−κε . Therefore,  x a+ε

123

n i (θ j + ε, x)d x = e−κε

 x a

n i (θ j , x)d x = e−κε Ii (θ j ).

Discrete limit and monotonicity properties …

We have also, for all i  0, 

 x ε

n i+1 (θ j + ε, x)d x = 0

 = =

ε ε

0 ε

n i+1 (θ j + ε − x, 0)e−

x 0

K  (θ j +ε−x+s,s)ds

dx

n i+1 (θ j + ε − x, 0)d x n i+1 (θ j + s, 0)d x.

0

Since supp n i (θ j , ·) ⊂ [a, ∞[ and ε < a, we also have  n i+1 (θ j + s, 0) = 0

 =

∞ ∞

K  (θ j + s, x)n i (θ j + s, x)d x κn i (θ j , x)e−κs = Ii (θ j )κeκs .

a

We integrate and obtain finally  x ε

n i+1 (θ j + ε, x)d x = (1 − e−κε )Ii (θ j ).  

This completes the proof of the lemma.

To conclude the proof of Lemma 3.7, we shall need the following binomial type representation. From Lemma 3.9, using the fact that Ii (θ j + ) = Ii (θ j+1 ), for all j, we can deduce ∀i, Ii (θ j + ε) = I0 (0)



j (1 − eκε )i (e−κε ) j−i , i

  where np denotes the binomial coefficients (with value 0 if p > n). We finally estimate



∞

2i

i

0

n i ( pa + 1, x)d x =

2i Ii ( pa + 1).

i

It follows from the construction of K  that Ii ( pa + 1) = Ii (θ Na + ε). Therefore,

i



∞

i

2

0

n i ( pa

+ 1, x)d x = I0 (0)

i

i

2

Na (1 − eκε )i (e−κε ) Na −i i

= (2(1 − e−κε ) + e−κε ) Na = (2 − e−κε ) Na .

123

S. Gaubert, T. Lepoutre

Inequality (40) can then be read κ (a)( p +1) a

I0 (0)(2 − e−κε ) Na  eλ

.  

This ends the Proof of Lemma 3.7.

Step 2: Construction of the sequence θi , with ε = 0 We need now to prove the existence of the sequence θi used in Lemma 3.7. First we prove the following Lemma 3.10 There exists a sequence θi as in Lemma 3.7 with ε = 0. The following tk ’s define the dynamics of successive division times of a cell initially of age x changing generation/dividing whenever it reaches age a at a time which is allowed. These will be used afterwards through the construction of θi . ⎧ if x  a, ⎨0 if a − x ∈ [τ, 1[ +N, t0 (x) = 0, t1 (x) = a − x (41) ⎩ a − x otherwise  if ti (x) ∈ [τ, 1[ +N, t (x) + a ti+1 (x) = i ti (x) + a otherwise We now describe the behavior of the tk ’s. First of all, from the definition of Na , pa , we have the following properties Lemma 3.11 The sequence ti satisfies the following properties: • for i < Na , ti (0) = ia, t Na (0) = pa + 1, • if ti0 (x) ∈ N , then for 0  n < Na , ti0 +n (x) = ti0 (x) + na, • for x > a, for any i  0, ti+1 (x) = ti (0), • x → ti (x) is nonincreasing for any i, • ∀i, ti (x)  ti (0)  ti+1 (x). Furthermore if x > 1 − τ or (x = 1 − τ and Na a = pa + 1), then we have t Na (x) < pa + 1. Proof The first and second point is an immediate consequence of the construction of Na , pa . The third point follows from the fact that if x  a, then t1 (x) = 0, and of the definition of Na . The last two points of the list are straightforward. The main difficulty states in the last statement. Firstly, this is obvious if x  a since in this case t Na (x) = t Na −1 (0). Now, suppose x < a, let the ti (x) be defined by (41) and i 0 := inf i. ti (x)∈N

Then we have by definition ∀n < Na , ti0 +n (x) = ti0 (x) + na, ∀i < i 0 , ti (x) = ia − x, ti0 (x) = i 0 a − x

123

Discrete limit and monotonicity properties …

As a consequence, we know that if i 0 < Na (taking n = Na − i 0 ) or if i 0 > Na , we have t Na (x) ∈ N and t Na (x)  t Na (0) ∈ N, therefore, in this case t Na (x) < t Na (0) = pa + 1. If i 0 = Na , then, since x > 1 − τ , Na a − x  pa + 1 − x < pa + τ. Therefore, we have necessarily t Na (x) < pa + τ < pa + 1 which ends the proof of Lemma 3.11 and thereby of Lemma 3.10.   Step 3: Construction of the sequence θi , with ε > 0. Here, we use the notation ar = sup{a  , λ∞ (a  ) = λ∞ (a)} and al = inf{a  , λ∞ (a  ) = λ∞ (a)}. ar −a 2

Lemma 3.12 Suppose there exists a sequence θ1 , . . . , θ Na and that ⎧ ⎨ θ1  max(a − 1 + τ, 0), ∀i, [θi , θi + ε] ⊂ [θi , θi  + τ ], ⎩ ∀i < Na , θi+1  θi + ε + a.

 ε  0, such

Then there exists a sequence θ1 , . . . , θ N a ⎧  θ  max(a − 1 + τ, 0), ⎪ ⎪ ⎨ 1 ∀i, [θi , θi + ε ] ⊂ [θi , θi  + τ ], ⎪ ⎪ ⎩ ∀i < N , θ   θ  + ε + a  . a

with ε = min(ε +

i+1

i

a−a  2 , τ ).

Proof We start the proof with an obvious but necessary remark: if λ∞ is constant on ]al , ar [, then • ετ • ar − al  τ We build the sequence θi following this procedure:

a − a . θi = max θi , θi − 2 With this construction, we have θ1  a  − 1 + τ and since ε −

a−a  2

 ε,

θi + ε  max(θi  + ε , θi + ε)  θi  + τ.

123

S. Gaubert, T. Lepoutre 

Note that we also have θi +ε  θi +ε + a−a 2 . This helps us to check the last property:  θi+1 − θi − a  − ε  θi+1 −

a − a a − a − θi − a  − ε − = θi+1 − θi − ε  0. 2 2  

This ends the proof of the lemma.

To complete the proof of the Theorem 1.2, we combine the previous arguments. If λ∞ is constant on [al , ar [, then for a < ar we can build a sequence θi as in Lemma 3.7 for any ε < ar 2−a (we first build a sequence for ar − 0 and ε = 0 and then use Lemma 3.12 to build a sequence with ε = ar − a − 0). Using Lemma 3.7, we get the inequality (39) for any ε < (ar − a/2). By a standard continuity argument, the same is true for ε = ar 2−a = r (a) as defined in (14). Finally, taking the logarithm of (39) with ε = r (a), we get λκ (a)  λ∞ (a)

log(2 − e−κr (a,τ ) . log 2

This leads to (14) thanks to the following inequality, which is a consequence of the concavity of the function log: ∀x ∈ [0, 1],

log(2 − x)  (log 2)(1 − x).

This ends the proof of the Theorem 1.2.

 

3.5 Numerical illustrations We next illustrate Theorem 1.2 by numerical experiments, that we produced using again monotone finite difference schemes (see the discussion concerning Fig. 8 above). It should be noted that the bound of the rate of convergence given in this theorem vanishes when Na a is a multiple of the period T . This is confirmed by the numerical experiments: the corresponding values of a are discontinuity points of the staircase function λ∞ (a), and the convergence is seen to be slow at the left of each of these points (Figs. 1, 8 and 9).

4 Concluding remarks We presented in this manuscript some results about a simplified time dependent cell division model. We showed that the principal eigenvalue which leads the growth of the system exhibits unexpected properties. We showed in particular that the growth rate is determined, in the limit of a large division rate κ, by synchronization phenomena involving the maturation age and the length of the interval in which division is allowed. We believe this could be generalized to more detailed cell cycle models, corresponding to refinements of the present discrete dynamics. This might lead to some interesting

123

Discrete limit and monotonicity properties …

3log2

2log2

log2 0.2

0.3

0.4

0.5

0.6

0.7

0.8 0.9 1.0 minimal age for division a

Fig. 8 Convergence for τ = 1/2 and T = 1. Limit in dashed blue and κ = 10, 30, 50 (in black, green and red respectively) (color figure online) 2log2

log2

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9 1.0 1.1 1.2 minimal age for division a

Fig. 9 Convergence for τ = 1/2 and T = 1. Limit in dashed blue and κ = 10, 30, 50 (in black, green and red respectively) (color figure online)

new questions concerning the synchronization between phases of the cycle. More precisely, we believe our tools could be applied to some cycle models of the form ⎧ ⎪ ∂t n i + ∂x n i + K i ψi (t)1[ai ,+∞) (x)n i (t, x)d x, ⎪ ⎪ ⎨ ∞ n i+1 (t, 0) = K i ψi (t) ai n i (t, x)d x, ⎪ ⎪ ⎪ ⎩n (t, 0) = K ψ (t)  ∞ n (t, x)d x. 1 I I I aI The functions ψi could be then chosen of the form ψi (t) = 1(τ l ,τ r ) (t). i

i

123

S. Gaubert, T. Lepoutre

In this case, establishing the discrete dynamics would be harder (only one the τi could be translated to 0) because the behavior should take into account many parameters: • minimal duration of each phases (ai ) • size of the possible transition interval (τir − τil ), l • synchronization of those intervals (position of τi+1 with respect to τil,r and ai ). The discrete dynamics could help us to analyze the effect of the slowing of one of the phases (for instance by modifying one of the ai ). Whereas the general analytical computation of the limiting growth rate appears to be a challenging combinatorial problem in the case of several phases, specific computations for given values of the time parameters, like the one of Fig. 6, can be easily performed, allowing one to detect the synchronizations that rule the growth rate. Also, we believe that the convergence to a limit governed by a discrete dynamical system, which we established here only for one phase models, is a general result, which holds for more complex models. This is left for future work. Acknowledgments The authors thank Jean Clairambault and Benoît Perthame for numerous insightful remarks and discussions. Most of this work has been achieved during the PhD of the second author at Laboratoire Jacques Louis Lions (Paris) and INRIA Rocquencourt.

Appendix A: Existence theory for the Floquet eigenvalue We next extend the proofs of Michel et al. (2005, section 5) to the present situation in which division rates may vanish. We shall use the following notation: • τh B = x → B(x + h), extended to take the value zero on [0, h] if h < 0, • L∞ per (0, T, X ) is the space of bounded T -periodic functions taking values in X , similarly, Cper (0, T, X ) is the space of T -periodic continuous (with respect to the time variable) functions taking values in X . Theorem A.1 Assume that a > 0, κ > 0, that ψ is nonnegative, bounded, not identically zero and that B positive, bounded, satisfying  ∀t, 0

∞



∞

ψ(t − x)B(x)d x =

ψ(t + x)B(x)d x = +∞.

(42)

0

1 + Then, there exists a unique λ F > 0 such that there exists (N , φ) in L ∞ per (0, T, L (R )) ∞ + × Cper (0, T, L (R )) satisfying (6–7–8). Furthermore, φ is unique.

Before starting the proof, we give a few remarks on the hypotheses: • we do not need regularity assumptions on the function ψ, ∞ • when min(ψ) > 0, the condition (42) is equivalent to 0 B = +∞, in the case min(ψ) = 0 it is for instance satisfied if min(B) > 0 or at least liminf+∞ B > 0, it is not optimal but other conditions would need assumptions on a and κ. It could

123

Discrete limit and monotonicity properties …

be replaced for instance by  ∀t,

∞



∞

κψ(t − x)B(x)d x,

a

κψ(t + x)B(x)d x > log 2,

a

but, as we are studying asymptotic properties, we rather restrict the study to a case where existence does not depend on a nor on κ, The proof of Theorem A.1 is based on the method of characteristics and on the Krein-Rutman Theorem, along the lines of Michel et al. (2005), but we need more precisions in order to relax the regularity assumptions on ψ. For the sake of simplicity  ∞ we take κ = 1. We first recall that N (t, x), φ(t, x) are give We set P(t) := a B(x)N (t, x)d x. Using the methods of characteristics, we have the following integral equations: 

x

, N (t, x) = N (t − x, 0)e−λx− 0 K (t−x+s,s)ds y ∞ − x (λ+K (t+y  −x,y  )dy  ) φ(t, x) = 2 x K (t + y − x, y)e φ(t + y − x, 0)dy.

(43) in particular, we notice that (provided λ  0) φ is bounded by 2φ(t, 0) L ∞ per (0,T ) since we can compute  φ(t, x)  2φ(t, 0) L ∞ per (0,T )



 2φ(t, 0) L ∞ per (0,T )  2φ(t, 0) L ∞ . per (0,T )

∞ x

∞

K (t + y − x, y)e− K (t + y, y + x)e−

y x

x 0

(λ+K (t+y  −x,y  )dy  ) K (t+y  ,y  +x)dy 

dy

dy

0

Therefore we are led to consider the eigenvalue problems in terms of integral equations on N (t, 0), φ(t, 0).  N (t, 0) = 2 =:

x

ψ(t − x + s)B(s)ds d x

a

 ψ(t − x)B(x)P(t − x) exp −λx −

a λ L2 (P)(t),  ∞

φ(t, 0) = 2 =:

 ψ(t)B(x)N (t − x, 0) exp −λx −

a λ L1 (N (·, 0))(t),  ∞

P(t) = 2 =:

∞

ψ(t − x + s)B(s)ds d x

x

a

 ψ(t + x)B(x)φ(t + x, 0) exp −λx −

a λ L3 (φ(·, 0))(t).

x

ψ(t + s)B(s)ds d x

a

We have defined three linear operators on L ∞ per (0, T ). These are well defined as soon λ λ as λ > 0. Moreover, we can see L2 and L3 as operators on the space C per (0, T ) of T-periodic continuous functions. We have

123

S. Gaubert, T. Lepoutre

Lemma A.2 Under the assumptions of Theorem A.1, for any λ > 0, Lλ2 and Lλ3 are nonnegative, compact linear operators on C per (0, T ). Proof The non negativity and linearity are obvious. We next show the compactness. If we fix λ, a > 0 then for f continuous and T -periodic with  f   1, if we write g = Lλ3 ( f ), we have 

∞

g(t + h) = 2 =2

a+h  ∞

 ψ(t + x)B(x − h) f (t + x) exp λ(h − x) −

a





 ψ(t + x)B(x − h) f (t + x) exp λ(h − x) −

a+h  x



a+h

 ψ(t + s)B(s)ds d x

x

 ψ(t + s)B(s)ds d x

a+h



a

ψ(t + s)B(s)ds d x.

x

ψ(t + x)B(x − h) f (t + x) exp λ(h − x) −

−2

a+h

This leads to g(t + h) − g(t)

 x  a+h ψ(t + x)B(x − h) f (t + x) exp λ(h − x) − ψ(t + s)B(s − h)ds d x = −2 a a+h

 ∞  x ψ(t + x)B(x − h) f (t + x) exp −λx − ψ(t + s)B(s − h)ds d x + 2(eλh − 1) a a+h

 x  ∞ ψ(t + x)(B(x − h) − B(x)) f (t + x) exp −λx − ψ(t + s)B(s − h)ds d x +2 a a+h 

 ∞ x x − a+h ψ(t+s)B(s−h)ds −λx e ψ(t + x)B(x) f (t + x)e − e− a ψ(t+s)B(s)ds d x, +2 a

= I1 + I2 + I3 + I4 .

We deal with each Ii separately, |I1 |  2ψBeλh h, |I2 |  2ψB(eλh

(44)

e−λa , − 1) λ

(45)

to make I3 small, we make the following remark: for any R > a,  |I3 |  2ψτh B − B L 1 ([a,R]) + 4ψB

∞

e−λx d x.

(46)

R

Observe that for all a and R, lim h→0+ τh B − B L 1 ([a,R]) = 0 (indeed, this is obvious if B is continuous, and, by density of continuous functions in L 1 ([a, R]), the same is true under the present assumptions). Since λ is positive, the second term in (46) can be made arbitrarily small by choosing R large enough, and then, the first term can be made arbitrarily small by choosing h small enough. It follows that for all > 0, we have |I3 |  ε for h small enough.

123

Discrete limit and monotonicity properties …

The same method can be applied to I4 . Finally we have the equicontinuity of Lλ3 (B) where B is the unit ball for the supremum norm. Thanks to Arzela-Ascoli theorem, Lλ3 is compact. We are now in position to apply Krein-Rutman theorem, there exist U2 , U3 nonnegative eigenvectors of Lλ2 , Lλ3 associated to their respective spectral radii ρ2 (λ), ρ3 (λ). We have now Lemma A.3 ρi (λ), Ui > 0,

(47)

Lλ1 (ψU2 )

(48) (49)

= ρ2 (λ)ψU2 ,

ρ3 = ρ2 .

Proof The second statement is a straightforward computation, we prove the first by contradiction: if U2 vanishes, then, for some t, 

∞

2

 ψ(t − x)B(x)U2 (t − x) exp −λx −

a

x

ψ(t − x + s)B(s)ds d x = 0,

a

as B > 0, it would mean thanks to T -periodicity ψ B ≡ 0 which would lead to U2 ≡ 0. Finally, as U2 does not vanish then Lλ2 (U2 ) > 0 and therefore, ρ2 (λ) > 0. The equality comes from the duality of operators L1 and L3 , we have 

T

ρ2



T

ψ(t)U2 (t)U3 (t)dt =

0

0

 L1 (ψU2 ).U3 (t)dt =

T

ψU2 L3 (U3 )dt

0



T

= ρ3

ψ(t)U2 (t)U3 (t)dt,

0

therefore ρ2 = ρ3 = ρ.

 

To end the proof, we need to find λ such that ρ(λ) = 1. Obviously, ρ is a decreasing function that vanishes at infinity. Since  x

ψ(t + x)B(x) exp − ψ(t + s)B(s)ds d x 2 a a

∞   x = 2, = 2 exp − ψ(t + s)B(s)ds 

∞

a

a

ρ → 2 as λ → 0. Therefore there exists a unique λ satisfying ρ(λ) = 1. Up to a renormalization, φ and P are unique, and therefore so is N . This ends the proof of Theorem A.1.

123

S. Gaubert, T. Lepoutre

References Akian M, Bapat R, Gaubert S (1998) Asymptotics of the Perron eigenvalue and eigenvector using max algebra. C R Acad Sci Paris 327(Série I):927–932 Bacaër N (2007) Approximation of the basic reproduction number R0 for vector-borne diseases with a periodic vector population. Bull Math Biol 69(3):1067–1091 Bacaer N, Abdurahman X (2008) Resonance of the epidemic threshold in a periodic environment. J Math Biol 57(5):649–673. doi:10.1007/s00285-008-0183-1 Bacaër N, Guernaoui S (2006) The epidemic threshold of vector-borne diseases with seasonality. The case of cutaneous leishmaniasis in Chichaoua, Morocco. J Math Biol 53(3):421–436 Bacaër N, Ouifki R (2007) Growth rate and basic reproduction number for population models with a simple periodic factor. Math Biosci 210(2):647–658 Bernard S, Herzel H (2006) Why do cells cycle with a 24 hour period? Genome Inform 17(1):72–79 Bernard S, Cajavec-Bernard B, Lepoutre F, Tvi Herzel H (2010) Tumor growth rate determines the timing of optimal chronomodulated treatment schedules. PLoS Comput Biol 6(3):e1000,712. doi:10.1371/ journal.pcbi.1000712 Billy F, Clairambault J, Fercoq O (2013) Optimisation of cancer drug treatments using cell population dynamics. In: Ledzewicz U, Schttler H, Friedman A, Kashdan E (eds) Mathematical methods and models in biomedicine lecture notes on mathematical modelling in the life sciences. Springer, New York, pp 265–309. doi:10.1007/978-1-4614-4178-6_10 Clairambault J, Laroche B, Mischler S, Perthame B (2003) A mathematical model of the cell cycle and its control. Tech. Rep. 4892, Inria Clairambault J, Michel P, Perthame B (2006) Circadian rhythm and tumour growth. C R Acad Sci 342(1):17– 22 Clairambault J, Gaubert S, Perthame B (2007a) An inequality for the Perron and Floquet eigenvalues of monotone differential systems and age structured equations. C R Math Acad Sci Paris 345(10):549–554 Clairambault J, Michel P, Perthame B (2007b) A mathematical model of the cell cycle and its circadian control. In: Mathematical modeling of biological systems. Birkhäuser, Boston, pp 247–259 Clairambault J, Gaubert S, Lepoutre T (2009) Comparison of Perron and Floquet eigenvalues in age structured cell division cycle models. Math Model Nat Phenom 4(3):183–209, also eprint arXiv:0812.0803 Clairambault J, Gaubert S, Lepoutre T (2011) Circadian rhythm and cell population growth. Math Comput Modell 53(7–8):1558–1567. doi:10.1016/j.mcm.2010.05.034. http://www.sciencedirect.com/ science/article/B6V0V-508X3CV-4/2/a7e4728e98fecac9b8b1d2c13d53fd90 (Mathematical Methods and Modelling of Biophysical Phenomena) Coron JM, Gabriel P, Shang P (2015) Optimization of an amplification protocol for misfolded proteins by using relaxed control. J Math Biol 70(1–2):289–327. doi:10.1007/s00285-014-0768-9 De Mairan J (1729) Observation Botanique. In: Histoire de l’Académie royale des sciences, pp 35–36 Diekmann O, Heijmans HJAM, Thieme HR (1986) On the stability of the cell-size distribution. II. Time-periodic developmental rates. Comput Math Appl Part A 12(4–5):491–512 (hyperbolic partial differential equations, III) Doumic Jauffret M, Gabriel P (2010) Eigenelements of a general aggregation-fragmentation model. Math Models Methods Appl Sci 20(05):757–783. doi:10.1142/S021820251000443X Filipski E, King VM, Li X, Granda TG, Mormont M, Liu X, Claustrat B, Hastings MH, Levi F (2002) Host circadian clock as a control point in tumor progression. J Natl Cancer Inst 94(9):690–697 Filipski E, Innominato P, Wu M, Iacobelli XLS, Xian L, Levi F (2005) Effects of light and food schedules on liver and tumor molecular clocks in mice. J Natl Cancer Inst 97(7):507–517. doi:10.1093/jnci/dji083 Magal P (2001) Compact attractors for time-periodic age-structured population models. Electron J Differential Equations pp No. 65, p. 35 (electronic) Marshall A, Olkin I (1979) Inequalities: theory of majorization and its applications. Academic Press, New York Michel P, Mischler S, Perthame B (2005) General relative entropy inequality: an illustration on growth models. J Math Pures et Appl 84(9):1235–1260 Mormont MC, Levi F (2003) Cancer chronotherapy: principles, applications, and perspectives. Cancer 97(1):155–169. doi:10.1002/cncr.11040 Park E, Iannelli M, Kim M, Anita S (1998) Optimal harvesting for periodic age-dependent population dynamics. SIAM J Appl Math 58(5):1648–1666. doi:10.1137/S0036139996301180 Perthame B (2007) Transport equations in biology. Birkhäuser, Boston

123

Discrete limit and monotonicity properties … Perthame B, Touaoula T (2008) Analysis of a cell system with finite divisions. Boletín SEMA 44:55–79 Thieme HR (1984) Renewal theorems for linear periodic Volterra integral equations. J Integral Equ 7(3):253– 277 Thieullen P, Garibaldi E (2012) Description of some ground states by Puiseux techniques. J Stat Phys 146:125–180

123