This is the accepted version of the journal article: Busse, Jan-Erik; Cuadrado Gavilán, Sílvia; Marciniak-Czochra, Anna. «Local asymptotic stability of a system of integro-differential equations describing clonal evolution of a self-renewing cell population under mutation». Journal of mathematical biology, Vol. 84, Issue 1-2 (January 2022), art. 10. DOI 10.1007/s00285-021-01708-w This version is available at https://ddd.uab.cat/record/268838 under the terms of the license Journal of Mathematical Biology manuscript No. (will be inserted by the editor) Local asymptotic stability of a system of integro-differential equations describing clonal evolution of a self-renewing cell population under mutation Jan-Erik Busse · Sı́lvia Cuadrado · Anna Marciniak-Czochra? Received: date / Accepted: date Abstract In this paper we consider a system of non-linear integro-differential equations (IDEs) describing evolution of a clonally heterogeneous population of malignant white blood cells (leukemic cells) undergoing mutation and clonal selection. We prove existence and uniqueness of non-trivial steady states and study their asymptotic stability. The results are compared to those of the system without mutation. Existence of equilibria is proved by formulating the steady state problem as an eigenvalue problem and applying a version of the Krein-Rutmann theorem for Banach lattices. The stability at equilibrium is analysed using linearisation and the Weinstein-Aronszajn determinant which allows to conclude local asymptotic stability. Keywords Selection mutation process · integro-differential equations · cell differentiation model · stationary solutions · asymptotic stability Mathematics Subject Classification (2010) MSC 92D25 1 Introduction The paper is devoted to the analysis of a system of integro-differential equa- tions (IDEs) describing clonal evolution of a self-renewing cell population. The ? Corresponding Author S.C. has been partially supported by the grants MTM2017-84214-C2-2-P and RED2018- 102650-T from MICINN. Research of A.M.-C. and J.-E.B. has been part of SFB 873 (Project B08) supported by German Research Foundation (DFG). Jan-Erik Busse · Anna Marciniak-Czochra Institute of Applied Mathematics Interdisciplinary Center for Scientific Computing (IWR) and BIOQUANT Center, Heidelberg E-mail: jan-erik.busse@web.de, anna.marciniak@iwr.uni-heidelberg.de Sı́lvia Cuadrado Departament de Matemàtiques, Universitat Autònoma de Barcelona E-mail: silvia@mat.uab.cat 2 Jan-Erik Busse, Sı́lvia Cuadrado, Anna Marciniak-Czochra?? population is heterogeneous with respect to the self-renewal ability of dividing cells, a property that may change due to cancerous mutations. The model is applied to describe dynamics of acute leukemias, an important type of malig- nant proliferative disorders of the blood forming system. Acute leukemias show a considerable inter-individual genetic heterogeneity and a complex genetic relationship among different clones, i.e. subpopulations consisting of genetically identical cells [32]. Similarly to the healthy hematopoi- etic system, the leukemic cell bulk is maintained by cells with stem-like proper- ties that can divide and give rise to progeny cells which either adopt the same cell route as the parent cell (undergo self-renewal) or differentiate to a more specialised cell type [32]. There exists theoretical [39] and experimental [34, 47] evidence suggesting that the self-renewal ability of leukemic stem cells has a significant impact on disease dynamics and patient prognosis [43]. Increased self-renewal confers a competitive advantage on cancer cell clones by leading to aggressive expansion of both stem and non-stem cancer cells and can be responsible for the clonal selection observed in experimental and clinical data [46]. The latter has been investigated using mathematical models of evolution of an arbitrary number of leukaemic clones coupled to a healthy cell lineage [8, 40]. A mathematical proof of clonal selection has been shown in Ref. [8] exploiting the analytical tractability of the model with a continuum of het- erogeneous clones differing with respect to the stem cell self-renewal ability. A similar result has been recently obtained in Ref. [30] for an extended model with two-parameter heterogeneity with respect to cancer stem cell self-renewal fraction and proliferation rate. It was shown that while increased proliferation rates may lead to a rapid growth of respective clones, the long-term selection process is governed by increased self-renewal of the most primitive subpopu- lation of leukemic cells [40, 30]. Mathematical analysis of the model provided an understanding of the link between the observed selection phenomenon and the nonlocal mode of growth control in the model, resulting from description of different plausible feedback mechanisms. Moreover, comparison of patient data and numerical simulations of the model allowing emergence of new clones suggested that self-renewal of leukemic clones increases with the emergence of clonal heterogeneity [42]. An open question is whether a mutation process with phenotypic heterogeneity in the course of disease may change the observed se- lection effect. To address this question, we propose an extension of the basic clonal selection model from Ref. [8] to account for the process of mutations. The model from Ref. [8] takes the form    ∂ a(x)    ∂t v(t, x) = 2 1 − 1+kρ(t) pu(t, x) − dv(t, x),    2a(x)  ∂ ∂t u(t, x) = 1+kρ(t) − 1 pu(t, x), (1) v(0, x) = v0 (x),     u(0, x) = u0 (x),  R where x ∈ Ω ⊂ R, a ∈ C(Ω), p, k, d ∈ R+ and ρ(t) = Ω v(t, x) dx. Stability of a system of equations describing clonal evolution of a cell population 3 The model describes evolution of one healthy cell lineage and an arbi- trary number of leukemic clones, where the structural variable x ∈ Ω repre- sents a continuum of possible cell clones (e.g. characterised by different gene expression levels) differing with respect to the self-renewal ability of divid- ing cells. We follow the convention that x = x0 ∈ Ω corresponds to healthy cells whereas different leukemic clones are characterised by different values of x ∈ Ω \{x0 }. The description of cell differentiation within each cell line is given by a two-compartment version of the multi-compartment system established in [33], mathematically studied in [26] and applied to patient data in [39]. The model focuses on self-renewal of primitive cells u(t, x) and their differentia- tion to mature healthy cells v(t, x0 ) or leukemic blasts v(t, x), x ∈ Ω \ {x0 }, which do not divide and for which the death rate has been denoted by d. The two-compartment architecture is based on a simplified description of the multi-stages differentiation process. A dividing cell gives rise to two progeny cells. A progeny cell is either more specialised than the mother cell, i.e., it is differentiated, or it is a copy of the mother cell (the case of self-renewal). The proliferation rate is denoted by the constant p. The function a(x) describes the fraction (probability) of self-renewal of cells of clone x, where dependence on x reflects the clonal heterogeneity. The feedback signal that promotes the self- 1 renewal of dividing cells is modelled using a Hill function 1+kρ(t) , where the parameter k > 0 is related to the degradation rate of the feedback signal [33]. This formula has been derived from a simple model of cytokine dynamics using a quasi-stationary approximation [21], motivated by biological findings presented in [27, 28]. Implementing other plausible regulation mechanisms led to a similar model dynamics that can reproduce the clinical observation [41, 44]. In [8] it has been shown that the solution (u, v) of system (1) converges in the space of positive Radon measures M+ (Ω) to a measure with support contained in the set of maximal values of the self-renewal fraction function a. In particular, if the set of maximal values of a consists solely of discrete points, the solution (u, v) converges to a sum of Dirac measures. In [30], it has been shown that a similar result holds for a model with multiple compartments under an additional assumption preventing Hopf bifurcation that may occur in the model with at least three-stages maturation structure [26]. The purpose of this paper is to extend the clonal selection model (1) to include a process causing phenotypic heterogeneity with respect to x. The process is called “mutation”, although mechanistically, it amounts to a con- tinuously distributed change of the phenotype in each round of division. Bio- logically, there exist at least two ways such process can be accomplished: (1) A continuous variable (“quantitative trait” in population genetics vocabulary [22]) can be a result of a superposition of discrete effects of mutations in a number of genes. Otherwise, it may be an effect of a discrete mutation in one gene, modulated by controlled or random gene expression varying from cell to cell, and hence “random” or “diffuse” [3]. Mathematically, this distinction 4 Jan-Erik Busse, Sı́lvia Cuadrado, Anna Marciniak-Czochra?? does not seem to matter much. As the aim of the model is to provide insights into multi-clonal cell dynamics under evolution of the stem cell self-renewal ability, we do not account for mechanistic details of the mutation process but focus on its mapping into the self-renewal parameter. There are essentially two ways to model mutations in a continuous setting. First, it can be done by adding a diffusive term, which may be introduced using a Laplacian, as in Ref. [35, 37, 29, 2]. Such models are based on the as- sumption that mutations can occur at all times within a cell’s life cycle and are not limited to division, as it is the case in epigenetic modifications. Con- sequently, such mutations do not change the overall number of individuals and only affect their distribution with respect to the structure variable. From a mathematical point of view the diffusive ansatz provides good properties of the obtained solution and allows using the library of methods for semi- linear parabolic equations. Nevertheless, as the references above indicate, there exists a difficulty with characterisation of the long-term behaviour. Alterna- tively, mutations can be modelled with an integral operator, see for instance [7, 9, 10, 23, 31]. This approach includes mutations that occur during prolifer- ation, which seems biologically realistic in case of genetic mutations [1]. In contrast to the diffusive ansatz, the mutation kernel approach can be used to model discontinuous changes in genotype or phenotype. For example, if for a given x, κ(x, y) = a1 1{(x−δ1 ,x+δ2 } + a2 1{(1−δ2 ,1) (y), then this is equivalent to modelling expected (mean) values of the effects of a change which avoids the interval (x + δ1 , 1 − δ2 ). This can be called a jump, while strictly speaking, the definition might require δ1 or δ2 being infinitesimally small. Mathematical disadvantage of the integral operators is the non-local structure which makes the analysis more complicated. Selection processes under mutation have been studied using both classes of mutation models, however using different methods to show convergence of the solutions, in the limit of small or rare mutations, to weighted Dirac measures. In case of a reaction-diffusion equation (RDE), the most common ansatz is to transform the RDE into a Hamilton-Jacobi equation, for which the viscosity solution provides the desired convergence result [35, 37]. For scalar equations, this transformation can be performed by setting ϕε (t, x) = ε ln(uε (t, x)) where uε (t, x) is the unknown of the RDE [18]. In case of a system of equations, the latter ansatz strongly depends on the structure of the model, since it is neces- sary to transform a system into a scalar Hamilton-Jacobi equation. In contrast, if mutation is modelled by an integral term, appropriate mathematical tools are provided by theory of positive semigroups and the infinite dimensional version of the Perron-Frobenius theorem. In this paper, we focus on the effect of rare mutations taking place during proliferation [1], and propose a model based on integro-differential equations where we assume that, during proliferation, a mutation occurs with probabil- ity ε. We prove, under suitable hypotheses, existence of locally asymptotically Stability of a system of equations describing clonal evolution of a cell population 5 stable steady states of the model, which converge, for ε → 0, to a weighted Dirac measure located at the point of maximum fitness of the corresponding pure selection model (1). The mathematical tools applied to analysis of the selection-mutation model are based on the ones used in [14, 15] and extended here to a system of two phenotype-structured populations in a non-compact domain. Similar results can be shown for small mutations that occur indepen- dently of the proliferation process [10]. The paper is structured as follows. In Section 2 we introduce the model with mutations and justify its well-posedness. Section 3 is devoted to existence and uniqueness of stationary solutions. Convergence of the steady states for the zero limit of mutation rate is studied in Section 4. In Section 5, we show local asymptotic stability of the stationary solution of the model. The paper is completed with two appendices providing technical proofs needed for the results in Section 5 and Section 3 respectively. 2 Selection-mutation model and its assumptions We extend system (1) to account for mutations described by an integral kernel operator and consider the following system of integro-differential equations,    ∂ a(x)   ∂t vε (t, x) = 2 1 − 1+kρ ε (t) puε (t, x) − dvε (t, x),      ∂ 2a(x) R u ε (t, x) = − (1 + ε) puε (t, x) + εp κ(x, y)uε (t, y) dy, ∂t 1+kρε (t) Ω (2) 0     v ε (0, x) = v (x), uε (0, x) = u0 (x),  R where ρε (t) = vε (t, x) dx. Ω As previously, uε (t, x) denotes the density of dividing cells structured with respect to the trait x that represents the expression level of genes influencing self-renewal ability of the cells, while vε (t, x) denotes the resulting mature cells for x = x0 ∈ Ω or leukemic blasts of clone x for x ∈ Ω \ {x0 }. The growth terms describing self-renewal and differentiation of dividing cells, regulated by a nonlocal nonlinear feedback from all non-dividing cells are taken from model (1) in Ref. [8]. Additionally, the model accounts for mutations that take place during proliferation at a rate ε ∈ (0, 1]. If a mutation occurs, then the proba- bility density that an individual with trait y mutates into one with trait x is denoted by κ(x, y). After scaling (multiplying ρε , vε and uε by k) one can assume that k = 1 and after scaling time, one can assume that the proliferation rate is p = 1 with d being replaced by θ = dp . Without loss of generality we consider then the 6 Jan-Erik Busse, Sı́lvia Cuadrado, Anna Marciniak-Czochra?? following system of integro-differential equations:    ∂ a(x)   ∂t vε (t, x) = 2 1 − 1+ρ ε (t) uε (t, x) − θvε (t, x),      ∂ 2a(x) R uε (t, x) = − (1 + ε) uε (t, x) + ε κ(x, y)uε (t, y) dy, ∂t 1+ρε (t) Ω (3)  v (0, x) = v 0 (x),    ε uε (0, x) = u0 (x),  R where ρε (t) = vε (t, x) dx. Ω In the remainder of this paper, we make the following assumptions: Assumption 1 1. Ω ⊂ R is open and bounded. 2. a ∈ C 1 (Ω) with 0 < a(x) < 1 for all x ∈ Ω and there exists x∗ ∈ Ω such that a(x∗ ) > 12 (otherwise, as shown in [8], the population dies out). Moreover, there exists a single point x̄ ∈ Ω where the maximal value of the self- renewal function is attained, i.e. x̄ = argmaxx∈Ω a(x), ā = a(x̄). 3. u0 , v 0 ∈ L1 (Ω) with u0 , v 0 > 0. R 4. κ ∈ C(Ω × Ω) is nonnegative and such that Ω κ(x, y)dx = 1 . Moreover, we assume that there exists δ > 0 such that the support of κ(·, y) contains the interval (y − δ, y + δ) for each y ∈ Ω. In particular, for small δ > 0 the model describes small mutations (small values of the distance between the traits) that occur with positive probability. System (3) can be written as     vε vε =A + f (u(t), v(t)), uε t uε where   −θ 0 A= R  0 −ε + ε Ω κ(x, y) · dy and f : L1 (Ω) × L1 (Ω) −→ L1 (Ω) 1  × L (Ω)    a(x) 2a(x) (vε (t), uε (t)) 7−→ 2 1 − 1+ρ ε (t) uε (t), 1+ρε (t) − 1 uε (t) . Since the operator A is the infinitesimal generator of a C 0 positive semigroup and f is locally Lipschitz, existence and uniqueness of local mild solutions of the initial value problem for (3) follow from application of the theory of semilinear evolution equations [36]. Notice that, since the right hand-side of (3) is a bounded and (locally) Lipschitz function in L1 (Ω) × L1 (Ω), the same result follows from an application of the Picard-Lindelöf theorem in Banach spaces [5]. Boundedness of total mass (see Appendix A) implies that solutions are global. Stability of a system of equations describing clonal evolution of a cell population 7 Numerical observation. Numerical simulations of the model suggest a selection effect, similar to that in the pure selection model (1). The difference is that, depending on the size of mutation frequency ε, we observe a distribution of the different cell clones around the one with the highest self-renewal fraction, see Figure 1. Convergence of the system to a solution concentrated around the most aggressive phenotype is a fast process and does not depend on initial data. The remainder of this paper is devoted to a rigorous proof of this observation. 3 Existence and uniqueness of non-trivial steady states In general, the stationary problem for selection-mutation equations can be reduced (see [9]) to a fixed point problem for a real function whose defini- tion depends on the existence and uniqueness of a dominant eigenvalue and a corresponding positive eigenvector of a certain linear operator (obtained by fixing the nonlinearity in the model). To solve the problem for system (3), we follow an approach proposed in Ref. [14] for the stationary problem of a predator-prey model consisting of an IDE coupled with an ordinary differen- tial equation (ODE). The structure of the ODE considered in [14] is a logistic type equation for which the steady state is given by a constant. Consequently, the steady state of the IDE depends on this constant that can be interpreted as a parameter. All together, the steady state problem can be reformulated as an eigenvalue problem, associated to the eigenvalue 0, for which a positive eigenfunction is sought. The latter still depends on the parameter given by the steady state of the ODE (which in [14] was the prey population at equilib- rium). To solve the coupled problem, it is necessary to choose the parameter in such a way that the eigenvalue problem and the steady state problem for the ODE are solved simultaneously. It results in solving a fixed point problem. In the remainder of this section, we adapt this approach to the cell population model (3) which consists of a system of two IDE’s. 3.1 Eigenvalue problem To find steady states of model (3), we consider the model obtained by inte- grating over the structure variable x the first equation in (3)  d R  a(x)     dt ρ ε (t) = 2 1 − 1+ρε (t) uε (t, x) dx − θρε (t),   Ω    ∂ 2a(x) R u (t, x) = − (1 + ε) uε (t, x) + ε κ(x, y)uε (t, y) dy,  ∂t ε 1+ρε (t) Ω (4)  uε (0, x) = u0 (x),    = v 0 (x) dx. R  ρε (0)   Ω Since the nonlinearity depends only on the total population of mature cells (the integral of the second variable), the integrated equation becomes an or- dinary differential equation for ρε (t). Consequently, the first component of a 8 Jan-Erik Busse, Sı́lvia Cuadrado, Anna Marciniak-Czochra?? Fig. 1: The simulations depict the solution uε (t, x) of model (2), for different values of ε. Going from top to bottom, the values of ε are 43 , 13 and 100 1 . The presented simulations were done for parameters k = 0.01, p = 0.9 and d = 0.2 and a mutation operator with a scaled Gaussian kernel. Stability of a system of equations describing clonal evolution of a cell population 9 steady state (ρε , uε (x)) of system (4) is a constant. Furthermore, the first com- ponent of the (corresponding) steady state of system (3) can be computed by inserting the steady state of system (4) (ρε , uε (x)) into the first equilibrium equation of (3) and solving it for vε . The (nontrivial) equilibria of system (4) are given by the solutions of  R  a(x)   0 = 2 1 − 1+ρε uε (x) dx − θρε ,  Ω  2a(x)  R (5)  0 = 1+ρε − (1 + ε) uε (x) + ε κ(x, y)uε (y) dy.  Ω Let us define for ρε > 0,   1 1 2a(x) Bε,ρε : L (Ω) → L (Ω), Bε,ρε uε (x) := − (1 + ε) uε (x), 1 + ρε Z Kε : L1 (Ω) → L1 (Ω), Kε uε (x) := ε κ(x, y)uε (y) dy, Ω 1 1 Cε,ρε : L (Ω) → L (Ω), Cε,ρε uε := Bε,ρε uε + Kε uε . (6) If a non-trivial steady state of system (4) exists, the first equation of system (5) provides a constant solution ρ ∈ (0, ∞). The second equation of system (5) can be then interpreted as an eigenvalue problem for Cε,ρ , which depends on the parameter ρ. Thus, we are looking for a function ϕε,ρ and a constant λε (ρ) such that Cε,ρ ϕε,ρ = λε (ρ)ϕε,ρ . (7) The first component of the steady state of system (4) is then given by the solution ρε of the equation λε (ρ) = 0. Denoting by ϕε,ρε the corresponding normalised eigenfunction, the second component of the steady state of (4) has the form uε = cε ϕε,ρε where cε ∈ (0, ∞) satisfies θρε cε = R   . a(x) 2 1− 1+ρε ϕε,ρε (x) dx Ω Let us observe that 0 is an eigenvalue with corresponding eigenfunction ϕε,ρ of Cε,ρ if and only if Bε,ρ ϕε,ρ + Kε ϕε,ρ = 0, ⇔ Kε ϕε,ρ = −Bε,ρ ϕε,ρ , −1  ⇔ Kε −Bε,ρ ψε,ρ = ψε,ρ 10 Jan-Erik Busse, Sı́lvia Cuadrado, Anna Marciniak-Czochra?? with ψε,ρ := −Bε,ρ ϕε,ρ . This means that ψε,ρ is an eigenfunction corresponding to eigenvalue 1 of the operator Tε,ρ : L1 (Ω) → L1 (Ω) given by −1 Tε,ρ u := Kε ◦ (−Bε,ρ u) Z 1+ρ (8) = ε κ(x, y) u(y) dy. (1 + ρ)(1 + ε) − 2a(y) Ω Remark 1 The equivalence of the eigenvalues 1 of the operator Tε,ρ (= Kε ◦ −1 (−Bε,ρ )) and 0 of the operator Cε,ρ (= Bε,ρ + Kε ) is exploited in [14], but the idea goes back to [6, Proposition 2.1]. We choose to study the eigenvalue problem for the operator Tε,ρ , because it allows a direct application of the version for Banach lattices of the well-known Krein-Rutmann theorem [16]. Alternatively, one can directly study the eigenvalue problem for the operator Cε,ρ that would require more work. Existence of a strictly dominant eigenvalue of Cε,ρ can be obtained from an application of a result of Greiner (Corollary 1.8 in Ref. [24]) that provides existence of an algebraically simple, strictly dominant eigenvalue of a perturbation B + K of the generator of a positive semigroup by a positive bounded irreducible operator K satisfying that there exists and integer n such that (KR(λ, B))n is compact for all λ with Reλ > s(B) and that s(B + K) > s(B). Remark 2 While analysis of the spectral bound of the operator Cε,ρ corre- sponds to study the Malthusian parameter (exponential growth rate) of the population, the approach chosen in this paper to analyse of the spectral ra- −1 dius of the so-called next-generation operator Tε,ρ (= Kε ◦ (−Bε,ρ )) leads to the basic reproduction number of the population [17], [45]. We proceed then to investigate the eigenvalue problem for the operator Tε,ρ . We define 2ā ρ̂ε := −1 (9) 1+ε where ā = maxx a(x) = a(x̄). Notice that Tε,ρ is a positive operator for ρ ∈ (ρ̂ε , +∞). Proposition 1 Let Tε,ρ be the linear operator defined in (8) for ρ ∈ (ρ̂ε , +∞). Its spectral radius r(Tε,ρ ) is an algebraically simple eigenvalue of Tε,ρ with a corresponding strictly positive eigenfunction. Moreover, r(Tε,ρ ) is the only eigenvalue of Tε,ρ having a positive eigenfunction. Proof By the Krein-Rutman theorem for Banach lattices, see [16, Theorem 12.3], the problem reduces to proving that Tε,ρ is a compact positive irre- ducible operator. Tε,ρ is a positive operator by definition and the choice of ρ. It is also evident that Tε,ρ is a bounded operator. Stability of a system of equations describing clonal evolution of a cell population 11 Let us start by proving that Tε,ρ is irreducible. Indeed, the assumption on the support of κ(·, y) implies that Z Z κ(x, y)dxdy > 0 Sc S for each measurable set S such that both S and S c have positive measure (see [10] for more details), whichRis a characterisation for the irreducibility ∞ of the kernel operator T f (x) := 0 κ(x, y)f (y)dy (see [38], V.6). Since, for 1+ρ ρ ∈ ρ̂ε , +∞) the function (1+ρ)(1+ε)−2a(y) is strictly positive, the operator Tε,ρ is irreducible. Applying [20, Corollary 5.1] with 1+ρ κ̃(x, y) := κ(x, y) , (1 + ρ)(1 + ε) − 2a(y) and  κ̃(x, y), y ∈ Ω κ̄(x, y) := , 0, y∈/Ω we obtain that Tε,ρ : L1 (Ω) → L1 (Ω) is compact if and only if for all ι > 0, there exist δ > 0, R > 0 such that for almost all x ∈ Ω and for every h ∈ R with |h| < δ Z Z |κ̄(x, y)| dy < ι, |κ̄(x, y + h) − κ̄(x, y)| dy < ι. R\BR (0) R Let ι > 0 be arbitrary but fixed. As Ω is bounded, let us choose R > 0 such that ι |Ω \ BR (0)| < . max 2 κ̃(x, y) (x,y)∈Ω Then, Z Z |κ̄(x, y)| dy = |κ̃(x, y)| dy ≤ max 2 κ̃(x, y) |Ω \ BR (0)| < ι. (x,y)∈Ω R\BR (0) Ω\BR (0) Due to the dominated convergence theorem and the continuity of κ̃, it holds Z Z |κ̄(x, y + h) − κ̄(x, y)| dy = |κ̃(x, y + h) − κ̃(x, y)| dy < ι, R Ω for |h| small enough, which completes the proof. 12 Jan-Erik Busse, Sı́lvia Cuadrado, Anna Marciniak-Czochra?? From the previous proposition we have that the operator Tε,ρ admits a strictly positive eigenfunction corresponding to the eigenvalue r(Tε,ρ ). What is left to show, in order to obtain equilibria, is that it is possible to choose ρ such that r(Tε,ρ ) = 1 and that this choice of ρ is unique (remember that showing r(Tε,ρ ) = 1 is equivalent to showing λε (ρ) = 0). The idea is to prove that r(Tε,ρ ) is continuous with respect to ρ and strictly monotone. The idea of using continuous dependence goes back to the late 1980s (see [6], [12]). Lemma 1 Let Tε,ρ be the linear operator defined in (8) for ρ ∈ (ρ̂ε , +∞). Its spectral radius, r(Tε,ρ ) is a continuous function of ρ. Moreover, there exists ε0 such that for ε < ε0 , r(Tε,ρ ) is strictly decreasing. Proof (Proof of Lemma 1) Continuity of r(Tε,ρ ) with respect to ρ follows from the continuity of a finite system of eigenvalues of a closed operator [25, Chapter IV, §3.5]. For the monotonicity, we use Gelfand’s formula for the spectral radius of a bounded linear operator A on a Banach space 1 r(A) = lim kAn k∞ n , n→∞ where k.k∞ is the operator norm. We have to show that 1 1 n n n n Tε,ρ1 ∞ > Tε,ρ 2 ∞ for ρ1 < ρ2 . A straightforward proof by induction provides the formula n−1 n R Q Tε,ρ u(x) = κ(x, yn ) κ(yi , yi+1 ) Ωn i=1 (1+ρ)n n Q u(yn ) dy1 · · · dyn . ((1+ρ)(1+ε)−2a(yi )) i=1 n In order to obtain monotonicity, we compute the derivative of Tε,ρ with respect to ρ. Another straightforward proof by induction (see Appendix C) shows that d n n n 1 dρ Tε,ρ < 0 for ε small enough. Thus Tε,ρ1 u > Tε,ρ2 u for all u ∈ L (Ω), u ≥ 0 and ρ1 < ρ2 . Then, taking the operator norm on both sides and using that the 1 function x 7→ x n is strictly monotone, we obtain 1 1 n n Tε,ρ 1 n ∞ ≥ Tε,ρ2 n ∞ for all n ∈ N ⇒ r(Tε,ρ1 ) ≥ r(Tε,ρ2 ) for ρ1 < ρ2 . The argument for strict monotonicity is the same as in Ref. [14]. Up to this point we have showed that the spectral radius r(Tε,ρ ) is a continuous and strictly decreasing function of ρ. Hence it is necessary to prove that there exists some ρ ∈ (ρ̂ε , ∞) such that r(Tε,ρ ) = 1. This is provided by Lemma 2 For all ε < ε0 , there exists a unique ρ ∈ (ρ̂ε , ∞) such that r(Tε,ρ ) = 1. Stability of a system of equations describing clonal evolution of a cell population 13 Proof We observe that Z ε lim Tε,ρ = κ(x, y)uε (y) dy := Tε . ρ→∞ 1+ε Ω Tε is a bounded operator, hence the spectrum is bounded. Since r(Tε ) ≤ kTε k∞ we obtain r(Tε ) < 1 . As r(Tε,ρ ) is continuous with respect to ρ, we can find ρ1 ∈ (ρ̂ε , ∞) such that r(Tε,ρ1 ) < 1. On the other hand, since Z 1+ρ kTε,ρ k∞ = sup ε κ(x, y) u(y) dy dx u∈L1 (Ω) ((1 + ρ)(1 + ε) − 2a(y)) kukL1 (Ω) =1 Ω 1 n n and r(Tε,ρ ) = lim Tε,ρ we have that limρ→ρ̂ε r(Tε,ρ ) = +∞. Lemma 1 n→∞ and the intermediate value theorem imply the statement. We can now formulate the theorem giving existence and uniqueness of steady states of system (3). Theorem 1 There exists some ε0 > 0 such that for all ε < ε0 there exists a unique, non-trivial steady state (vε , uε ) of system (3). Proof Combination of Proposition 1 and Lemma 2 yields the existence of a unique, non-trivial steady state (ρε , uε ) of system (4). Substituting it in the first equilibrium equation in system (3), we obtain the first component of its unique, non-trivial steady state (vε , uε ). 4 Convergence of the steady states Once we have proved, for ε small enough, existence of a stationary solution (vε , uε ) of system (3), we are interested in the behavior of this steady state when the mutation rate goes to zero. In order to study it, we recall the equivalence between the eigenvalue prob- lem for the operators Cε,ρ = Bε,ρ + Kε and Tε,ρ = Kε (Bε,ρ )−1 which is Bε,ρ ϕε,ρ + Kε ϕε,ρ = λε (ρ)ϕε,ρ ⇔ Kε (λε (ρ) − Bε,ρ )−1 ψε,ρ = ψε,ρ , for ψε,ρ = (λε (ρ) − Bε,ρ )ϕε,ρ . That is, ϕε,ρ is an eigenfunction of Cε,ρ corre- sponding to eigenvalue λε (ρ) if and only if 1 is an eigenvalue of Kε (λε (ρ) − Bε,ρ )−1 with eigenfunction ψε,ρ . The proof of convergence of the steady state consists of the following steps: We begin by showing that λε (ρ) is a strictly dominant eigenvalue of Cε,ρ with a corresponding strictly positive eigenfunction. Then we show that for 14 Jan-Erik Busse, Sı́lvia Cuadrado, Anna Marciniak-Czochra??   ε → 0, λε (ρ) → maxx∈Ω 2a(x)1+ρ − 1 pointwise, from which we conclude con- vergence of the corresponding eigenfunctions ϕε,ρ . Additionally, convergence of the eigenvalues λε (ρ) allows deducing convergence of ρε (first component of the steady state of system (4)), and ultimately convergence of the steady state. 4.1 Existence of eigenvalues Following [6, Theorem 2.2], in order to show that λε (ρ) is a strictly dominant eigenvalue of Cε,ρ with a corresponding strictly positive eigenfunction, it is sufficient to find some   2a(x) λ1 > s(Bε,ρ ) = max − (1 + ε) x∈Ω 1+ρ such that r(Kε (λ1 Id−Bε,ρ )−1 ) > 1, (where s(A) and r(A) denote respectively the spectral bound and the spectral radius of a linear operator A). We use the following characterisation of the spectral bound of the generator A of a strongly continuous positive semigroup s(A) ≥ sup{µ ∈ R : Af ≥ µf for some 0 < f ∈ D(A)}. (10) This property of the spectral bound is stated in Ref. [4] in C(K), the space of all real-valued continuous functions on a compact space K, but the proof also holds for any generator of a positive semigroup in a Banach lattice such that the spectral bound and the growth bound coincide, which is the case in Lp -space, 1 ≤ p < ∞, [49]. Proposition 2 There exists λ1 > s(Bε,ρ ) such that r(Kε (λ1 Id−Bε,ρ )−1 ) > 1. Proof Let us define function q : R × Ω → R by   2a(y) q(λ, y) := λ − − (1 + ε) . 1+ρ Then, we need to prove that Z ∃λ1 > s(Bε,ρ ) ∃g ∈ L1 (Ω), g > 0 ∀x ∈ Ω : ε κ(x, y)q(λ1 , y)−1 g(y) dy > g(x). Ω   2a(x) Observe that argminx∈Ω q(s(Bε,ρ ), x) = argmaxx∈Ω 1+ρ − (1 + ε) = x̄. Using the definition of q, we obtain  q(s(Bε,ρ ), x̄) = 0, ∂ (11) ∂y q(s(B ε,ρ ), y) |y=x̄ = 0. Stability of a system of equations describing clonal evolution of a cell population 15 Choosing g = χBδ (x̄) , for a small δ > 0, we estimate Z Z ε κ(x, y)q(λ1 , y)−1 g(y) dy = ε κ(x, y)q(λ1 , y)−1 dy Ω Bδ (x̄) Z ≥ε min 2 κ(x, y) q(λ1 , y)−1 dy.(12) x,y∈[x̄−δ,x̄+δ] Bδ (x̄) Expanding the function q using the Taylor formula up to the 0th order around x̄ yields q(λ1 , y) = q(λ1 , x̄) + o (ky − x̄k) . (13) Inserting equation (13) into inequality (12) leads to ε κ(x, y)q(λ1 , y)−1 g(y) dy ≥ R Ω 1 R ε min κ(x, y) q(λ1 ,x̄)+o(ky−x̄k) dy > 1, x,y∈[x̄−δ,x̄+δ]2 Bδ (x̄) by using the first equation of (11), choosing δ small enough and λ1 close enough to s(Bε,ρ ). 4.2 Convergence of the eigenvalues and the eigenfunctions Now that we have proved existence of a strictly dominant eigenvalue λε (ρ) of the operator Cε,ρ we can formulate a result about its limiting behavior. Lemma 3 Let λε (ρ) be the strictly dominant eigenvalue of the operator Cε,ρ defined in (6), then   ε→0 2a(x) λε (ρ) −→ max −1 . x∈Ω 1+ρ   Proof For notational simplicity, we denote µ(x) := 2a(x) 1+ρ − 1 ,   µε (x) := 2a(x) 1+ρ − (1 + ε) . We want to show that for all δ > 0 there exists ε0 such that for all ε < ε0 it holds λε (ρ) ∈ Bδ (µ(x̄)), where recall that x̄ denotes the point where the maximal value of the self- renewal function is attained. We begin by proving that λε (ρ) ≤ µ(x̄). The assumptions on κ imply that Z 1 0 0 ∀u ∈ L (Ω), u ≥ 0 ∃Ω ⊂ Ω ∀x ∈ Ω : κ(x, y)u(y) dy ≤ u(x) (14) Ω 16 Jan-Erik Busse, Sı́lvia Cuadrado, Anna Marciniak-Czochra?? where Ω 0 is not a Lebesgue null set. Assuming otherwise, Z κ(x, y)u(y) dy − u(x) > 0, Ω for R almost all x, and integrating R Rthe left hand side with respect to x, using Ω κ(x, y)dx = 1, we obtain Ω Ω κ(x, y)u(y)dy − u(x)dx = 0, which is a contradiction. Taking f1 (x) := χΩ 0 ϕε,ρ (x) where ϕε,ρ is the positive eigenfunction corre- sponding to the eigenvalue λε (ρ) and Ω 0 being a set of positive measure such that (14) holds for ϕε,ρ ,   2a(x) R λε (ρ)f1 (x) = 1+ρ − 1 f1 (x) − εf1 (x) + ε Ω 0 κ(x, y)f1 (y)dy   2a(x) ≤ 1+ρ − 1 f1 (x). Hence, by (10) we conclude that λε (ρ) ≤ µ(x) ≤ µ(x̄). We show now that for all δ > 0, there exists ε0 such that it holds λε (ρ) ≥ µ(x̄) − δ for any ε < ε0 . Let δ be such that µε (x) ≥ µε (x̄)−δ for x ∈ Ω 0 , where Ω 0 ⊂ Ω is a suitably chosen set. Let us consider a smooth function u with supp u ⊂ Ω 0 . Then, due to the positivity of Kε , it holds Cε,ρ u = Bε,ρ u + Kε u ≥ µε (x)u ≥ (µε (x̄) − δ)u. By inequality (10), we obtain λε (ρ) ≥ µε (x̄) − δ. Furthermore, we know |µ(x) − µε (x)| ≤ ε ≤ δ. Thus, we conclude ∀δ > 0 ∃ ε0 > 0 ∀ε < ε0 : λε (ρ) ∈ B2δ (µ(x̄)). Proposition 3 Let ϕε,ρ be the unique positive eigenfunction corresponding to the eigenvalue λε (ρ) of the operator Cε,ρ defined in (6). Then ϕε,ρ * δx̄ in M+ (Ω), as ε → 0, where * denotes narrow convergence and x̄ is the unique value where the maximum of the self-renewal function a(x) is attained. Stability of a system of equations describing clonal evolution of a cell population 17 Remark 3 We use the notion of the narrow convergence. A sequence (µn )n∈N ⊂ M(S), where S is a metric space, converges narrowly to a measure µ ∈ M(S) iff Z Z lim ψ(x) dµn (x) = ψ(x) dµ(x) ∀ψ ∈ Cb (S). n→∞ S S In probability theory, the latter convergence is also known as the weak con- vergence of measures. Narrow convergence is metrizable with the flat norm (bounded Lipschitz distance), i.e., the convergence result can be understood in a suitable metric in the space of positive Radon measures. We refer to [19] for different notions of convergence in measures. Proof To prove convergence of the steady state, we make an ansatz that the eigenfunctions form a Dirac sequence. By Proposition 2 and [6, Theorem 2.2], the eigenfunction ϕε,ρ of Cε,ρ is strictly positive, ϕε,ρ > 0 for all ε > 0. As an eigenfunction in L1 (Ω), it can be normalised to kϕε,ρ kL1 (Ω) = 1. It remains to show that Z ε→0 ϕε,ρ (x) dx −→ 0 Ωc for Ω c ⊂ Ω with x̄ ∈ / Ω c and dist(x̄, Ω c ) > 0.   According to Lemma 3, it holds λε (ρ) → maxx∈Ω 2a(x) 1+ρ − 1 for ε → 0. Hence, for any Ω c as defined above, it is possible to choose ε < ε0 such that   2a(x) ∀x ∈ Ω c : λε (ρ) > −1 . (15) 1+ρ Integrating the eigenvalue problem for Cε,ρ , we obtain R  2a(x)  0= 1+ρ − (1 + ε) ϕε,ρ (x) − λε (ρ)ϕε,ρ (x) dx Ωc R R +ε κ(x, y)ϕε,ρ (y) dy dx  Ω c Ω   ≤ maxc 2a(x) R 1+ρ − 1 − λ ε (ρ) ϕε,ρ (x) dx+ x∈Ω Ωc R R ε κ(x, y)ϕε,ρ (y) dy dx. Ωc Ω By inequality (15), the first term is negative and bounded, hence by rearrang- ing the inequality we obtain, for a constant C > 0, Z Z Z C ϕε,ρ (x) dx ≤ ε κ(x, y)ϕε,ρ (y) dy dx ≤ ε. Ωc Ωc Ω Consequently, ϕε,ρ is a Dirac sequence and converges subsequently to a Dirac measure concentrated in x̄. 18 Jan-Erik Busse, Sı́lvia Cuadrado, Anna Marciniak-Czochra?? 4.3 Convergence of the steady states Now we prove convergence of the steady states (ρε , uε ) for ε → 0. As a first step, we show Proposition 4 Let ρε be the unique zero of λε (ρ), hence the first component of the steady state of system (4) and let ρ¯0 := 2ā − 1. Then ε→0 ρε −→ ρ̄0 . Proof Notice that  ρ¯0 is the unique zero of the decreasing function λ0 (ρ) := 2a(x) max 1+ρ − 1 . Lemma 3 and the fact that λ0 (ρ) changes sign (recall that, by x∈Ω Assumption 1 there exists x∗ ∈ Ω such that a(x∗ ) > 21 ) prove the statement. The next result provides convergence of the family of steady states of system (4). Theorem 2 Let (ρε , uε ) be the family of stationary solutions of System (4). Then uε * ρ̄1 δx̄ in M+ (Ω), ρε → ρ̄0 in R, where * denotes narrow convergence in measures (see Remark 3), ρ̄0 = 2ā−1, ρ̄1 = θ(2ā−1) and x̄ is the isolated point where the maximum of the self-renewal function a(x) is attained x̄ = arg max a(x). x∈Ω Proof Convergence of ρε has already been proven in Proposition 4. Convergence of uε is done in two steps. By construction uε = cε ϕε,ρε , where ϕε,ρε is the unique (normalised) eigenfunction corresponding to the zero eigen- θρε value of the operator Cε,ρε defined in (6) and cε = R a(x) . 2(1− 1+ρε )ϕε,ρε (x) dx Ω Thus, we want to prove that both, the constants cε and the eigenfunctions ϕε,ρε , converge. The same argument as in Proposition 3 yields convergence of ϕε,ρε to the Dirac delta located at x̄. Because of this property of ϕε,ρε and convergence of ρε , it follows θρε ε→0 θρ¯0 cε = R   −→   = ρ̄1 . a(x) ā 2 1− 1+ρε ϕε,ρε dx 2 1− 1+ρ¯0 Ω This concludes the proof. We can finally formulate the result giving the behavior for small mutation rate of the equilibria of System (3). Stability of a system of equations describing clonal evolution of a cell population 19 Theorem 3 Let (vε , uε ) be the family of stationary solutions of System (3). Then, for ε → 0, (vε , uε ) * (ρ̄0 δx̄ , ρ̄1 δx̄ ) in M+ (Ω), where * denotes narrow convergence in measures (see Remark 3), ρ̄0 = 2ā−1, ρ̄1 = θ(2ā−1) and x̄ is the unique value where the maximum of the self-renewal function a(x) is attained. Proof Theorem 2 provides convergence of the second component of the steady state. The first component can be written as   2 a(x) vε (x) = 1− uε (x) =: gε (x)uε (x). θ 1 + ρε Proposition 4 implies that for any f ∈ Cc   ε→0 2 a(x) gε (x)f (x) −→ 1− f (x) =: g0 (x)f (x) strongly, θ 1 + ρ¯0 which, by Proposition 3.13 in [5], implies that ε→0 hgε (x)f (x), uε (x)i −→ hg0 (x)f (x), ρ̄1 δx̄ i, that is, vε * g0 (x̄)ρ̄1 δx̄ = (2ā − 1)δx̄ which concludes the proof. 5 Stability of the steady states Selection-mutation equations can be written, in a general way, in the form ∂ z(t, x) = Aε (F (z))z, (16) ∂t with F being a linear function from the state space to an m-dimensional space and such that Aε (E) is a linear operator for a fixed E = F (z). Assuming that equation (16) has a semilinear structure and the spectral mapping property holds (i.e., the growth bound of a semigroup is equal to the spectral bound of its generator, which is the case in L1 ), the principle of linearised stability [48] yields local asymptotic stability of a steady state if the spectrum of the corresponding linearisation is located entirely in the open left half plane. A stability result for equation (16) is provided in Ref. [11]. It is shown using the principle of linearised stability and the fact that, in case of finite dimensional nonlinearity, the linearised operator at the steady state is a degenerated perturbation of a known operator with spectral bound equal to 0. This reduces the computation of the spectrum of the linearisation to the 20 Jan-Erik Busse, Sı́lvia Cuadrado, Anna Marciniak-Czochra?? computation of zeroes of the so-called Weinstein-Aronszajn determinant [25]. System (4) can be written in form (16) with F : R × L1 (Ω) → R, (ρ, u) 7→ ρ, (17) and  R       −θ 2 1 − a(x) 1+ρ · dx ρ Ω Aε F = Aε (ρ) =  .     u 0 2a(x) − (1 + ε) + ε R κ(x, y) · dy 1+ρ Ω 0 Note that operator Aε (ρ) generates a C positive semigroup. Indeed, Aε (ρ) can be written in the following way  R    a(x) 0 2 1 − · dx ! −θ  0 1+ρ Aε (ρ) = + Ω R ,    2a(x) 0 1+ρ − (1 + ε) 0 ε κ(x, y) · dy Ω where the first term is the generator of a C 0 positive semigroup and the sec- ond term is a linear positive operator, thus the sum generates a C 0 positive semigroup. Linearising system (4) at the steady state zε := (ρε , uε ), i.e., taking a ∂ perturbation z = zε +z̄, applying the stationary equation ∂t zε = 0 and Taylor’s formula, we obtain  0   ρ̄ ρ̄ ∂ ū = (Ã ε + Sε ) , ∂t ū where  R  a(x)   −θ 2 1 − 1+ρε · dx Ω Ãε = Aε (ρε ) =  ,     2a(x) R 0 1+ρε − (1 + ε) + ε κ(x, y) · dy R 2a(x)u  Ω (18) ε (1+ρε ) 2 0 Sε = Ω . 2a(x)uε − (1+ρ ε) 2 0 with Ãε , Sε defined in R × L1 (Ω). We also define the following “limit” operators in R × M+ (Ω)  R    a(x) −θ 2 1 − 1+ ρ¯0 · dx Ã0 = A0 (ρ¯0 ) =  Ω  ,    2a(x) 0 1+ρ¯0 − 1 (19) 2āρ¯1 ! (1ρ¯0 )2 0 S0 = 2āρ¯1 . − (1+ ρ¯0 )2 δx̄ 0 Stability of a system of equations describing clonal evolution of a cell population 21 where ρ¯0 and ρ¯1 are given by Theorem 2. As mentioned before, our aim is to apply the stability result given in Ref. [11] to system (4). For the sake of completeness, we summarize the two relevant theorems [11, Theorem 1 and 2] into the following theorem (for m = 1 which is the case for our model): Theorem 4 Let zε be a non-trivial positive steady state of equation(16), where F is a linear function from the state space to an m-dimensional space and, for a fixed E = F (z), Aε (E) is a generator of a C 0 positive semigroup on the state space. Let Ãε + Sε be the linearisation of Aε at the equilibrium zε . Let ωε (λ), ω0 (λ) be the Weinstein-Aronszajn determinants for Ãε + Sε and Ã0 + S0 , respectively and D := {λ ∈ C |<(λ) ≥ 0, λ 6= 0 }. Let ωε (λ), ω0 (λ) be holomorphic functions in D such that ω0 (λ) does not vanish in D and ε→0 ωε (λ) −−−→ ω0 (λ) (20) uniformly in λ on compact sets in D. Additionally, assume that 1 ∃L > 0 ∀ |λ| > L : Sε R(λ, Ãε ) < . (21) ∞ 2 If 0 is a strictly dominant eigenvalue of Ãε with algebraic multiplicity 1, Pε is the projection onto the eigenspace of the eigenvalue 0 and   F (Pε Sε zε ) 6= 0 and lim inf+ λF (Ãε − λ)−1 Sε zε 6= 0, (22) (ε,λ)→(0 ,0) then for ε small enough the steady state zε is locally asymptotically stable. Theorem 5 Let Assumption 1 hold and additionally, let κ be separable in its variables, i.e., ∃κ1 , κ2 ∈ C(Ω) : κ(x, y) = κ1 (x)κ2 (y). Then, for ε small enough, the steady state (ρε , uε ) of system (4) is locally asymptotically stable. The proof of this theorem is a direct application of Theorem 4. Since it is technical, it is deferred to Appendix B. Theorem 6 Let Assumption 1 hold and additionally, let κ be separable in its variables, i.e. ∃κ1 , κ2 ∈ C(Ω) : κ(x, y) = κ1 (x)κ2 (y). Then, for ε small enough, the steady state (vε , uε ) of system (3) is locally asymptotically stable. 22 Jan-Erik Busse, Sı́lvia Cuadrado, Anna Marciniak-Czochra?? Proof By Theorem 5, it only remains to prove the result for the first component of the steady state vε (x). It holds ∂ vε (t, x) = F (uε (t, x), ρε (t)) − θvε (t, x) ∂t   a(x) with F (uε (t, x), ρε (t)) := 2 1 − 1+ρ ε (t) uε (t, x) and Z t vε (t, x) = v0 (x)e−θt + F (uε (τ, x), ρε (τ ))e−θ(t−τ ) dτ. 0 Then, since (vε (x), uε (x)) is an equilibrium of system (3), it follows Rt kvε (t, x) − vε (x)kL1 (Ω) = Ω v0 (x)e−θt + 0 F (uε (τ, x), ρε (τ ))e−θ(t−τ ) dτ R F (uε (x), ρε ) − dx θ Rt |v0 (x)e−θt |dx + e−θ(t−τ ) R R ≤ Ω 0 Ω F (uε (τ, x), ρε (τ )) R −F (uε (x), ρε ) dxdτ + Ω F (uε (x), ρε ) dx Rt e−θ(t−τ ) dτ − 1  0 θ which tends to zero as t → ∞ due to Theorem 5 (for more details on the second term see the proof of Lemma 7 in [8]). 6 Discussion Acute myeloid leukemia is an aggressive cancer of the blood forming system. The malignant cell population is composed of multiple clones that evolve over time. Clonal data reflect the mechanisms linked to clonal properties such as self-renewal ability and proliferation of leukemic stem cells, which govern treat- ment response and relapse. In particular, recent insights from experimental and theoretical models suggest that self-renewal not only influences the stem cell population but has a crucial impact on dynamics of non-stem cells. Increased self-renewal confers competitive advantage on cancer cell clones by leading to aggressive expansion of both stem and non-stem cancer cells. Recent models suggest that self-renewal is the key parameter to understand clonal competi- tion, selection and emergence of resistance in cancer cell populations. In this paper, in order to study the impact of self-renewal on evolution of intra-patient cancer heterogeneity, we propose a mathematical model that accounts for dynamics of healthy and cancerous white blood cell populations governed by self-renewal and differentiation processes that exhibit clonal het- erogeneity. The model describes evolution of a continuum of cell clones struc- tured by the self-renewal potential. It is an extension of the model introduced Stability of a system of equations describing clonal evolution of a cell population 23 in [8], where it was shown that experimentally observed clonal selection can be explained by heterogeneity of self-renewal capacity of leukemic stem cells. Clones with the highest stem cell self-renewal are selected due to a nonlocal competition for regulatory factors such cytokines. Our aim is to check if the clonal selection acts also in the system perturbed by mutations. The latter in- crease clonal heterogeneity in the course of disease by introducing new clones with perturbed stem cell self-renewal capacity. The model we consider is a system of integro-differential equations for the populations of primitive and mature cells, both structured by their mutation stage x. Clonal heterogeneity is modelled by assuming the self-renewal fraction to be a function of the struc- ture variable x. Heterogeneity emergence is described by an integral operator with mutation probability ε. Applying a version of the Krein-Rutmann theo- rem, we show in Section 3 existence of steady states of our selection-mutation model, which are proved, in Section 4, to converge as ε → 0 to a weighted Dirac mass at the point of maximum fitness of the pure selection model studied in [8]. Finally, we show local asymptotic stability of the steady states. The convergence holds in the sense of narrow convergence in Radon mea- sures. In general, we cannot expect the strong (norm- total variation) conver- gence of the solution to a stationary solution. For initial data in L1 (Ω), the model solutions for any finite time point are uniformly continuous with respect to the Lebesgue measure. As uε * ρ̄1 δx̄ narrowly, the distance between the two solutions T V (uε , ρ̄1 δx̄ ) ≥ 2ρ̄1 . Alternatively, the problem can be treated by considering convergence with respect to an appropriate metric, for example flat metric (bounded Lipschitz distance), for details see [8, 19] and references therein. Our results demonstrate that mutations do not prevent the clonal selection and dominance of clones with the highest stem cell self-renewal potential. The persistent clones exhibit some heterogeneity with respect to the self-renewal values located near its maximum. This range of heterogeneity becomes smaller as the mutation rate decreases, which means that, for small mutation rates, model dynamics coincides with dynamics observed in [8] with initially assumed clonal heterogeneity. Certainly, new mutations leading to an increased stem cells self-renewal in the course of disease will lead to different clonal com- position of cancer relapse comparing to the primary manifestation. However, as discussed in [8, 40] our models do not allow distinguishing between new mutations after the first therapy and persistence of preleukemic cells (small populations) that cannot be detected at the diagnosis but may dominate in relapses. In summary, the selection-mutation model proposed in this paper provides a mathematically rigorous proof that mutations do not change the selection process. Hence, our model analysis confirms the “mechanistic” expla- nation for the clonal selection phenomena observed in acute leukemias, which is based on a non-local cytokine-based regulation of stem cell self-renewal. Im- portance of this observation in the context of leukemia evolution, emergence of resistance in response to cytotoxic chemotherapy and dynamics of the disease relapses has been discussed in [8, 40] and in a broader context in the review paper [43]. In particular, model-based analysis of acute myeloid leukemia pa- 24 Jan-Erik Busse, Sı́lvia Cuadrado, Anna Marciniak-Czochra?? tients data suggested that the increased self-renewal is correlated with a poor patient prognosis [39]. Acknowledgments The authors are greatly indebted to the Referees for their helpful comments. 7 Appendix A R Lemma 4 Let ν(t) = Ω u(t, x) dx. Under the Assumptions 1, the function U (t) = ν(t) + ρ(t) is uniformly bounded on R × Ω. Proof The equation for U (t) reads, for t > 0 Z   Z   ! ∂ 1 2a(x) ν(t) a(x) ν(t) U (t) = − 1 u(t, x) dx − 2 2 1− u(t, x) dx + θ ∂t ρ(t) Ω 1 + ρ(t) ρ (t) Ω 1 + ρ(t) ρ(t) ! ≤ U (t) 2ā + θ − 2 (1 − ā) U (t) , (23) where the latter estimate holds because   2a(x) − 1 + θ ≤ 2ā + θ 1 + ρ(t) and a(x) 1− > 1 − ā. 1 + ρ(t) Since the right-hand side of equation (23) is a logistic type nonlinearity, we conclude that   2ā + θ U (t) ≤ max U (0), =: M1 ∀ (t, x) ∈ R+ × Ω. 2(1 − ā) By definition of U we can infer that ν(t) ≤ M1 ρ(t) ∀ (t, x) ∈ R+ × Ω. Lemma 5 Under the Assumptions 1, the total masses ν(t) and ρ(t) are uni- formly bounded. Proof Integrating the second equation in (4) and applying Lemma 4 we obtain ! d 2ā ν(t) ≤ − 1 ν(t) dt 1 + ν(t)M1 Stability of a system of equations describing clonal evolution of a cell population 25 from which we can conclude that ν(t) ≤ max {ν(0), (2ā − 1)M1 } =: M2 . From the first equation in (4) and using the previous inequality, we obtain d ρ(t) ≤ 2ν(t) − θρ(t) ≤ 2M2 − θρ(t) dt and thus   2M2 ρ(t) ≤ max ρ(0), . θ 8 Appendix B Here we provide the proof of Theorem 5. The proof is divided into several parts, each dealing with a different assumption of the stability theorem 4. 8.1 Convergence of the Weinstein-Aronszajn determinants The operators Sε and S0 defined in (18) and (19) are one-dimensional range operators with basis Z 2a(x)uε  2āρ¯1   2 dx  (1 + ρε )  and  (1 + ρ¯0 )2  (24) Ω  2āρ¯1 2a(x)uε − (1+ρε )2 − (1+ ρ¯0 )2 δx̄ , respectively. Therefore the Weinstein-Aronszajn determinants   ωε (λ) := det Id + Sε R(Ãε , λ)|rg(Sε )   (25) ω0 (λ) := det Id + S0 R(Ã0 , λ)|rg(S0 ) are well defined. In the next lemma we prove convergence result (20). Lemma 6 Let ωε (λ), ω0 (λ) be the Weinstein- Aronszajn determinants defined in (25). Then, ε→0 ωε (λ) −→ ω0 (λ) uniformly in λ ∈ D = {λ ∈ C |<(λ) ≥ 0, λ 6= 0 }. Both ωε (λ) and ω0 (λ) are holomorphic in D. 26 Jan-Erik Busse, Sı́lvia Cuadrado, Anna Marciniak-Czochra?? Proof In order to prove convergence, we estimate     |ωε (λ) − ω0 (λ)| ≤ det Id + Sε R(Ãε , λ)|rg(Sε ) − det Id + Sε R(Ã0 , λ)|rg(Sε )     + det Id + Sε R(Ã0 , λ)|rg(Sε ) − det Id + S0 R(Ã0 , λ)|rg(S0 ) . For the first term on the right-hand side it can be shown, in the same way as in [15, proof of Proposition 1] that ε→0 Sε (Ãε − λ)−1 −1 |R(Sε ) − Sε (Ã0 − λ)|R(Sε ) −→ 0. ∞ It remains to prove that     ε→0 det Id + Sε R(Ã0 , λ)|R(Sε ) − det Id + S0 R(Ã0 , λ)|R(S0 ) −→ 0. For this purpose we compute the determinant explicitly. The basis of rg(Sε ) is given by (24). Then, a direct computation yields Z 1 2a(x)uε (x) Sε R(Ã0 , λ)|R(Sε ) = − dx θ+λ (1 + ρε )2 Ω Z   1 a(x) 1 2a(x)uε (x) + 2 1−   dx. θ+λ 1 + ρ¯0 2a(x) − 1 − λ (1 + ρε )2 Ω 1+ρ¯0 According to Theorem 2, we know that ρε converges strongly to ρ¯0 in R and uε converges narrowly to ρ¯1 δx̄ in M+ (Ω). Hence,   ε→0 1 2āρ¯1 1 ā 2āρ¯1 Sε R(Ã0 , λ)|R(Sε ) −−−→ − − 2 1 − θ + λ (1 + ρ¯0 )2 λ(θ + λ) 1 + ρ¯0 (1 + ρ¯0 )2 = S0 R(Ã0 , λ)|R(S0 ) . By definition of the Weinstein-Aronszajn determinant, ωε (λ) and ω0 (λ) are holomorphic in D, see [25, p. 245]. 8.2 Boundedness of Sε R(λ, Ãε ) Lemma 7 There exists a constant L > 0 such that for all |λ| > L 1 Sε R(λ, Ãε ) < . ∞ 2 Proof Since supε<ε0 Ãε and supε<ε0 kSε k∞ are bounded, we obtain for ∞ |λ| > 2 Ãε ∞ ∞  X n kSε k∞ 2 kSε k∞ Sε R(λ, Ãε ) = Sε λ−1 λ−1 Ãε ≤ ≤ . ∞ n=0 λ − Ãε |λ| ∞ ∞ n o Choosing L > max 2 Ãε , 4 kSε k∞ leads to the assertion. ∞ Stability of a system of equations describing clonal evolution of a cell population 27 8.3 Proof of hypotheses (22) (excluding 0 and values with small positive real part from the spectrum) Lemma 8 For the steady state zε := (ρε , uε (x)) of system (4), it holds F (Pε Sε zε ) 6= 0. (26) Proof From the definition of F given in formula (17), it is sufficient to show that Pε Sε zε 6= 0. Since 0 is a simple strictly dominant eigenvalue of operator Ãε , we can decompose the space L1 (Ω) = hzε i Range(Ãε ) (see Theorem L A.3.1 in [13]). Hence, we have to prove that Sε zε ∈ / Range(Ãε ) what is equiv- alent to showing that D  ρ∗   E ρε ε , Sε 6= 0, (27) u∗ε uε where (ρ∗ε , u∗ε ) is the eigenfunction corresponding to the eigenvalue 0 of the adjoint operator Ã∗ε . The adjoint operator reads    −θ  0 Ãε =  2 1 − a(x) 2a(x) R 1+ρε − (1 + ε) + ε κ(y, x) · dy  1+ρε Ω and we obtain that ρ∗ε = 0. This implies that u∗εis an eigenfunction  correspond- 2a(x) R ing the the zero eigenvalue of the operator 1+ρε − (1 + ε) · +ε κ(y, x)· Ω which is the adjoint operator of Bε,ρ + Kε defined by formulas (6). This oper- ator is a generator of an irreducible positive semigroup in the Banach lattice L1 (Ω), as it is the perturbation by an irreducible operator of the generator of a positive semigroup. By Proposition 3.5 in [4] we obtain that u∗ε is strictly positive, which, together with the fact that ρ∗ε = 0, implies that D  ρ∗  2a(x)uε (x)ρε u∗ε (x)  E Z ε ρε , Sε =− 6= 0 u∗ε uε Ω (1 + ρε )2 The last step is to show Lemma 9 For the steady state zε , it holds   lim inf λF (Ãε − λ)−1 Sε zε 6= 0. (28) (ε,λ)→(0+ ,0) Proof We start by computing the resolvent operator R(λ, Ãε ),     −1 1 R a(x) 2 1− R(λ, Cε,ρε ) · dx (Ãε − λ)−1 =  θ+λ θ+λ Ω 1+ρε , 0 R(λ, Cε,ρε ) 28 Jan-Erik Busse, Sı́lvia Cuadrado, Anna Marciniak-Czochra?? where recall that Cε,ρε is defined in (6). Condition (28) reads    ρ lim inf+ λF (Ãε − λ)−1 Sε ε (ε,λ)→(0 ,0)   u ε −λ R a(x) 2a(x)uε ρε = lim inf+ θ+λ 2 1 − 1+ρ ε R(λ, C ε,ρε ) (1+ρε )2 dx (ε,λ)→(0 ,0) Ω  + 2a(x)u ε ρε R (1+ρε )2 dx 6= 0. Ω Since the limit of the second term is zero, condition (28) becomes   −1 Z a(x) 2a(x)uε ρε lim inf 2 1− λR(λ, Cε,ρε ) dx 6= 0 (ε,λ)→(0 ,0) θ + λ + 1 + ρε (1 + ρε )2 Ω Determination of the limit is difficult. Since 0 is an eigenvalue of Cε,ρε , the limiting behaviour of λR(λ, Cε,ρε ) for λ tending to zero is not obvious, because the resolvent tends to infinity (Cε,ρε tends to a multiplication operator), while λ tends to zero. However, separation of variables of the kernel κ allows an explicit derivation of the resolvent R(λ, Cε,ρε ) which facilitates the computation of the previous limit. Under this assumption, we write   Z 2a(x) Cε,ρε u = − (1 + ε) u+εκ1 (x) κ2 (y)u(y) dy = −αε (x)u+εκ1 (x)Lu, 1 + ρε Ω   2a(x) where αε (x) := 1+ε− > 0, since ρε is the first component of the 1+ρε R steady state of system (4), and Lu = κ2 (y)u(y) dy. Following the scheme Ω proposed in [11, Section 4.2] for the explicit computation of the resolvent operator R(λ, Cε,ρε ), we obtain  −1 " 1 κ1 (x) R(λ, Cε,ρε )g = L −(αε (x) + λ)−1 g (29) ελ αε (x)(αε (x) + λ) + ε(αε (x) + λ)−1 gL (αε (x) + λ)−1 κ1 (x)  (30) # − ε(αε (x) + λ)−1 κ1 (x)L (αε (x) + λ)−1 g .  (31) Let us define   2a(x)ρε a(x) βε (x) := , H ε (x) := 2 1 − , (1 + ρε )2 1 + ρε in order to shorten the notational effort. Note that   a(x) 2a(x)ρ¯0 lim Hε (x) = 2 1 − =: H(x), lim βε (x) = =: β(x). ε→0 1 + ρ¯0 ε→0 (1 + ρ¯0 )2 Stability of a system of equations describing clonal evolution of a cell population 29 We need to determine the following limiting process Z lim inf+ Hε (x)λR(λ, Cε,ρε )βε (x)uε (x) dx =: lim inf Ξ(ε, λ). (ε,λ)→(0 ,0) (ε,λ)→(0+ ,0) Ω Substituting the expression of the resolvent operator derived in (29), we obtain Z  −1 1 κ1 (y) Ξ(ε, λ) = Hε (x) L · ε αε (y)(αε (y) + λ) Ω "   κ1 (y) −(αε (x) + λ)−1 βε (x)uε (x) + ε(αε (x) + λ)−1 βε (x)uε (x)L αε (y) + λ  # βε (y)uε (y) −ε(αε (x) + λ)−1 κ1 (x)L dx. αε (y) + λ For a better distinction between the terms, let us define Z  −1 1 κ1 (y) I := − Hε (x)L (αε (x) + λ)−1 βε (x)uε (x) dx, ε αε (y)(αε (y) + λ) Ω Z  −1   −1 κ1 (y) κ1 (y) II := Hε (x)(αε (x) + λ) βε (x)uε (x)L L dx, αε (y)(αε (y) + λ) αε (y) + λ Ω Z  −1   κ1 (y) βε (y)uε (y) III := − Hε (x)(αε (x) + λ)−1 κ1 (x)L L dx, αε (y)(αε (y) + λ) αε (y) + λ Ω Ξ(ε, λ) = I + II + III. Using the steady state equation κ1 (x) uε = ε Luε , (32) αε (x) we obtain Z  −1 βε (x)κ1 (x)Luε κ1 (y) I = − Hε (x) L dx αε (x)(αε (x) + λ) αε (y)(αε (y) + λ) Ω Z  −1 βε (x) κ1 (x)κ2 (x) κ1 (y) = −Luε Hε (x) L dx. κ2 (x) αε (x)(αε (x) + λ) αε (y)(αε (y) + λ) Ω The sequence gε denoted by  −1 κ1 (x)κ2 (x) κ1 (y) gε (x) := L , αε (x)(αε (x) + λ) αε (y)(αε (y) + λ) 30 Jan-Erik Busse, Sı́lvia Cuadrado, Anna Marciniak-Czochra?? defines a Dirac sequence. The definition also guarantees that Z ∀ ε > 0 : gε (x) > 0 and gε (x) dx = 1. Ω c / Ω c and dist(x̄, Ω c ) > Let x̄ = argmaxx∈Ω a(x) and take Ω ⊂ Ω such that x̄ ∈ 0. It follows from Theorem 2 that     2a(x) ε→0 2a(x) αε (x) = 1 + ε − −→ 1 − . 1 + ρε 1 + ρ¯0 We conclude that αε κ(x)(α 1 (x)κ2 (x) ε (x)+λ) converges and is subsequently bounded on c Ω .   Then, using 1 = εL ακ1ε(x) (x) , we estimate Z Z   κ1 (x)κ2 (x) 1 κ1 (x)κ2 (x) 1 κ1 dx ≥ dx = L αε (x)(αε (x) + λ) max(αε (x) + λ) αε (x) max(αε (x) + λ) αε Ω x∈Ω Ω x∈Ω 1 ε→0 = −→ ∞. ε max(αε (x) + λ) x∈Ω This implies that Z ∀ Ω c ⊂ Ω, x̄ ∈ / Ω c , dist(x̄, Ω c ) > 0 : gε (x) dx → 0 for ε → 0. Ωc Since we additionally know by Theorem 2 that uε converges narrowly, we infer Z Luε = κ2 (y)uε (y) dy → ρ̄1 κ2 (x̄) for ε → 0. Ω Thus, we obtain ε→0 I −→ −ρ̄1 H(x̄)β(x̄) < 0. Computing the limit for II and II and using equality (32), we obtain " Hε (x) Luκε2β(x) ε (x) κ1 (x)κ2 (x) R lim inf(II + III) = lim inf ε αε (x)(αε (x)+λ) (ε,λ)→(0+ ,0) Ω #  −1      κ1 (y) L αε (y)(αε (y)+λ) · L αεκ(y)+λ 1 (y) − αε (x) βε (x) L βε (y)κ1 (y) αε (y)(αε (y)+λ) dx =0 which concludes the proof. Remark 4 Note that the separation of variables of κ is needed only, because of the explicit computation of the resolvent R(λ, Cε,ρε ). All results up to this point do not need this assumption and work for Assumption 1 alone. Stability of a system of equations describing clonal evolution of a cell population 31 9 Appendix C In this appendix we prove that the operators n−1 n R Q Tε,ρ u(x) = κ(x, yn ) κ(yi , yi+1 ) Ωn i=1 (1+ρ)n n Q u(yn ) dy1 · · · dyn . ((1+ρ)(1+ε)−2a(yi )) i=1 d n defined in the proof of Lemma 1 satisfy dρ Tε,ρ < 0 for ε small enough. The differential operator and the integral can be interchanged, because of Leibniz’ integral rule. This implies denoting by dy = dy1 · · · dyn d n dρ (Tε,ρ u)(x) = n−1 (1 + ρ)n Z d Y κ(x, yn ) κ(yi , yi+1 ) Q n u(yn ) dy = dρ Ωn i=1 ((1 + ρ)(1 + ε) − 2a(yi )) i=1 n−1 (1 + ρ)n Z Y d κ(x, yn ) κ(yi , yi+1 ) n u(yn ) dy. dρ Q Ωn i=1 ((1 + ρ)(1 + ε) − 2a(yi )) i=1 Both κ and u are positive functions, so the sign of the derivative is solely determined by the derivative of the fraction. Performing the derivative yields d (1 + ρ)n n = dρ Q ((1 + ρ)(1 + ε) − 2a(yi )) i=1 n n(1 + ρ)n−1 Q ((1 + ρ)(1 + ε) − 2a(yi )) i=1  2 − n Q ((1 + ρ)(1 + ε) − 2a(yi )) i=1 n d (1 + ρ)n dρ Q ((1 + ρ)(1 + ε) − 2a(yi )) i=1  2 . n Q ((1 + ρ)(1 + ε) − 2a(yi )) i=1 Again we see that it is sufficient to look only at a small part of this derivative to determine the sign, namely the numerator. The claim is n n Y d Y n ((1 + ρ)(1 + ε) − 2a(yi )) − (1 + ρ) ((1 + ρ)(1 + ε) − 2a(yi )) < 0 i=1 dρ i=1 and can be shown by induction over n ∈ N. Let n = 1, then ((1 + ρ)(1 + ε) − 2a(y1 )) − (1 + ρ)(1 + ε) = −2a(y1 ) < 0, 32 Jan-Erik Busse, Sı́lvia Cuadrado, Anna Marciniak-Czochra?? by Assumption 1. Let the statement be true for n ∈ N. Then have a look at the derivative for n+1 n+1 Y (n + 1) ((1 + ρ)(1 + ε) − 2a(yi )) i=1 n+1 d Y −(1 + ρ) ((1 + ρ)(1 + ε) − 2a(yi )) dρ i=1 n Y =n ((1 + ρ)(1 + ε) − 2a(yi )) · ((1 + ρ)(1 + ε) − 2a(yn+1 )) i=1 n+1 Y + ((1 + ρ)(1 + ε) − 2a(yi )) i=1 " n d Y −(1 + ρ) ((1 + ρ)(1 + ε) − 2a(yi )) · ((1 + ρ)(1 + ε) − 2a(yn+1 )) dρ i=1 n # Y + ((1 + ρ)(1 + ε) − 2a(yi )) i=1 " n Y = ((1 + ρ)(1 + ε) − 2a(yn+1 )) n ((1 + ρ)(1 + ε) − 2a(yi )) i=1 n # d Y −(1 + ρ) ((1 + ρ)(1 + ε) − 2a(yi )) dρ i=1 "n+1 Y + ((1 + ρ)(1 + ε) − 2a(yi )) i=1 n # Y −(1 + ρ) ((1 + ρ)(1 + ε) − 2a(yi )) < 0, i=1 because the first term is negative due to the induction assumption and the second term is negative because 1 + ρ > (1 + ρ)(1 + ε) − 2a(yn+1 ) for ε small enough. References 1. Alberts, B., Bray, D., Hopkin, K., Johnson, A., Lewis, J., Raff, M., Roberts, K., Walter, P.: Essential cell biology, 4 edn. Garland Science (2013) 2. Almeida, L., Bagnerini, P., Fabrini, G., Hughes, B., Lorenzi, T.: Evolution of cancer cell populations under cytotoxic therapy and treatment optimisation: insight from a phenotype-structured model. ESAIM Math. Model. Numer. Anal. 53, 1157–1190 (2019) 3. Alon, U.: An introduction to systems biology: design principles of biological circuits. CRC press (2019) Stability of a system of equations describing clonal evolution of a cell population 33 4. Arendt, W., Grabosch, A., Greiner, G., Groh, U., Lotz, H.P., Moustakas, U., Nagel, R., Neubrander, F., Schlotterbeck, U.: One-parameter semigroups of positive operators, Lecture Notes in Mathematics, vol. 1184. Springer-Verlag, Berlin (1986) 5. Brezis, H.: Functional analysis, Sobolev spaces and partial differential equations. Uni- versitext. Springer, New York (2011) 6. Bürger, R.: Perturbations of positive semigroups and applications to population genet- ics. Mathematische Zeitschrift 197(2), 259–272 (1988) 7. Bürger, R., Bomze, I.M.: Stationary distributions under mutation-selection balance: structure and properties. Adv. in Appl. Probab. 28(1), 227–251 (1996) 8. Busse, J.E., Gwiazda, P., Marciniak-Czochra, A.: Mass concentration in a nonlocal model of clonal selection. Journal of mathematical biology 73(4), 1001–1033 (2016) 9. Calsina, À., Cuadrado, S.: Small mutation rate and evolutionarily stable strategies in infinite dimensional adaptive dynamics. Journal of mathematical biology 48(2), 135–159 (2004) 10. Calsina, À., Cuadrado, S.: Stationary solutions of a selection mutation model: The pure mutation case. Mathematical Models and Methods in Applied Sciences 15(07), 1091–1117 (2005) 11. Calsina, À., Cuadrado, S.: Asymptotic stability of equilibria of selection-mutation equa- tions. Journal of mathematical biology 54(4), 489–511 (2007) 12. Capasso, V., Thieme, H.: A threshhold theorem for a reaction-diffusion epidemic system. In: Differential equations and applications, Vol. I, II (Columbus, OH, 1988), pp. 128– 133. Ohio Univ. Press, Athens, OH (1989) 13. Clément, P., Heijmans, H.J.A.M., Angenent, S., van Duijn, C.J., de Pagter, B.: One- parameter semigroups, CWI Monographs, vol. 5. North-Holland Publishing Co., Ams- terdam (1987) 14. Cuadrado, S.: Equilibria of a predator prey model of phenotype evolution. Journal of Mathematical Analysis and Applications 354(1), 286–294 (2009) 15. Cuadrado, S.: Stability of equilibria of a predator prey model of phenotype evolution. Math. Biosci. Eng 6, 701–718 (2009) 16. Daners, D., Medina, P.K.: Abstract evolution equations, periodic problems and appli- cations, vol. 279. Chapman & Hall/CRC (1992) 17. Diekmann, O., Heesterbeek, J.A.P., Metz, J.A.J.: On the definition and the computation of the basic reproduction ratio R0 in models for infectious diseases in heterogeneous populations. J. Math. Biol. 28(4), 365–382 (1990) 18. Diekmann, O., Jabin, P.E., Mischler, S., Perthame, B.: The dynamics of adaptation: an illuminating example and a hamilton–jacobi approach. Theoretical population biology 67(4), 257–271 (2005) 19. Düll, C., Gwiazda, P., Marciniak-Czochra, A., Skrzeczkowski, J.: Spaces of Measures and their Applications to Structured Population Models. Cambridge University Press (2021) 20. Eveson, S.P.: Compactness criteria for integral operators in L∞ and L1 spaces. Proc. Amer. Math. Soc. 123(12), 3709–3716 (1995). DOI 10.2307/2161898. URL https://doi.org/10.2307/2161898 21. Getto, P., Marciniak-Czochra, A., Nakata, Y., Vivanco, M.: Global dynamics of two- compartment models for cell production systems with regulatory mechanisms. Mathe- matical Biosciences 245, 258–268 (2013) 22. Gillespie, J.H.: Population genetics: a concise guide. JHU Press (2004) 23. Greene, J., Lavi, O., Gottesman, M.M., Levy, D.: The impact of cell density and mu- tations in a model of multidrug resistance in solid tumors. Bull. Math. Biol. 76(3), 627–653 (2014) 24. Greiner, G.: A typical Perron-Frobenius theorem with applications to an age-dependent population equation. In: Infinite-dimensional systems (Retzhof, 1983), Lecture Notes in Math., vol. 1076, pp. 86–100. Springer, Berlin (1984) 25. Kato, T.: Perturbation theory for linear operators. Springer Science & Business Media (1984) 26. Knauer, F., Stiehl, T., Marciniak-Czochra, A.: Oscillations in a white blood cell produc- tion model with multiple differentiation stages. J. Math. Biol. 80(3), 575–600 (2019) 34 Jan-Erik Busse, Sı́lvia Cuadrado, Anna Marciniak-Czochra?? 27. Kondo, S., Okamura, S., Asano, Y., Harada, M., Niho, Y.: Human granulocyte colony- stimulating factor receptors in acute myelogenous leukemia. European Journal of Haematology 46, 223–230 (1991) 28. Layton, J., Hockman, H., Sheridan, W., Morstyn, G.: Evidence for a novel in vivo control mechanism of granulopoiesis: mature cell-related control of a regulatory growth factor. Blood 74, 1303–1307 (1989) 29. Lorenzi, T., Chisholm, R., Clairambault, J.: Tracking the evolution of cancer cell pop- ulations through the mathematical lens of phenotype-structured equations. Biology Direct 11, 1–17 (2016) 30. Lorenzi, T., Marciniak-Czochra, A., Stiehl, T.: Mathematical modeling of leukemogen- esis and cancer stem cell dynamics. J. Math. Biol. 79, 1587–1621 (2019) 31. Lorz, A., Lorenzi, T., Hochberg, M.E., Clairambault, J., Perthame, B.: Populational adaptive evolution, chemotherapeutic resistance and multiple anti-cancer therapies. ESAIM: Mathematical Modelling and Numerical Analysis 47(2), 377–399 (2013) 32. Lutz, C., Hoang, V.T., Buss, E., Ho, A.D.: Identifying leukemia stem cells - is it feasible and does it matter? Cancer Lett 338, 10–14 (2012) 33. Marciniak-Czochra, A., Stiehl, T., Ho, A.D., Jäger, W., Wagner, W.: Modeling of asym- metric cell division in hematopoietic stem cells-regulation of self-renewal is essential for efficient repopulation. Stem cells and development 18(3), 377–386 (2009) 34. Metzeler, K., Maharry, K., Kohlschmidt, J., Volinia, S., Mrozek, K., Becker, H., Nicolet, D., Whitman, S., Mendler, J., Schwind, S., Eisfeld, A., Wu, Y., Powell, B., Carter, T., Wetzler, M., Kolitz, J., Baer, M., Carroll, A., Stone, R., Caligiuri, M., Marcucci, G., Bloomfield, C.: A stem cell-like gene expression signature associates with inferior out- comes and a distinct microRNA expression profile in adults with primary cytogenetically normal acute myeloid leukemia. Leukemia 27(10), 2023–2031 (2013) 35. Mirrahimi, S.: Adaptation and migration of a population between patches. Discrete & Continuous Dynamical Systems-Series B 18(3) (2013) 36. Pazy, A.: Semigroups of linear operators and applications to partial differential equa- tions, Applied Mathematical Sciences, vol. 44. Springer-Verlag, New York (1983) 37. Perthame, B., Barles, G.: Dirac concentrations in lotka-volterra parabolic pdes. Indiana University Mathematics Journal 57(7), 3275–3301 (2008) 38. Schaefer, H.H.: Banach lattices and positive operators. Springer-Verlag, New York- Heidelberg (1974). Die Grundlehren der mathematischen Wissenschaften, Band 215 39. Stiehl, T., Baran, N., Ho, A., Marciniak-Czochra, A.: Cell division patterns in acute myeloid leukemia stem-like cells determine clinical course: a model to predict patient survival. Cancer Research 75, 940–949 (2015) 40. Stiehl, T., Baran, N., Ho, A.D., Marciniak-Czochra, A.: Clonal selection and therapy resistance in acute leukemias: Mathematical modelling explains different proliferation patterns at diagnosis and relapse. J. Royal Society Interface 11 (2014) 41. Stiehl, T., Ho, A., Marciniak-Czochra, A.: Cytokine response of leukemic cells has im- pact on patient prognosis: Insights from mathematical modeling. Scientific Reports 8, 2809 (2018) 42. Stiehl, T., Lutz, C., Marciniak-Czochra, A.: Emergence of heterogeneity in acute leukemias. Biology Direct 11(1), 51 (2016) 43. Stiehl, T., Marciniak-Czochra, A.: Stem cell self-renewal in regeneration and cancer: Insights from mathematical modeling. Current Opinion in Systems Biology 5, 112–120 (2017) 44. Stiehl, T., Wang, W., Lutz, C., Marciniak-Czochra, A.: Mathematical modeling provides evidence for niche competition in human aml and serves as a tool to improve risk stratification. Cancer Research 80(18), 3983–3992 (2020) 45. Thieme, H.R.: Spectral bound and reproduction number for infinite-dimensional popu- lation structure and time heterogeneity. SIAM J. Appl. Math. 70(1), 188–211 (2009) 46. Van Delft, F.W., Horsley, S., Colman, S., Anderson, K., Bateman, C., Kempski, H., Zuna, J., Eckert, C., Saha, V., Kearney, L., et al.: Clonal origins of relapse in etv6- runx1 acute lymphoblastic leukemia. Blood 117, 6247–54 (2011) 47. Wang, W., Stiehl, T., Raffel, S., Hoang, V., Hoffmann, I., Poisa-Beiro, L., Saeed, B., Blume, R., Manta, L., Eckstein, V., Bochtler, T., Wuchter, P., Essers, M., Jauch, A., Trumpp, A., Marciniak-Czochra, A., Ho, A., Lutz, C.: Reduced hematopoietic stem cell Stability of a system of equations describing clonal evolution of a cell population 35 frequency predicts outcome in acute myeloid leukemia. Haematologica 102(9), 1567– 1577 (2017) 48. Webb, G.F.: Theory of nonlinear age-dependent population dynamics. CRC Press (1985) 49. Weis, L.: The stability of positive semigroups on Lp spaces. Proc. Amer. Math. Soc. 123(10), 3089–3094 (1995)