Mass concentration in a nonlocal model of clonal selection

Self-renewal is a constitutive property of stem cells. Testing the cancer stem cell hypothesis requires investigation of the impact of self-renewal on cancer expansion. To better understand this impact, we propose a mathematical model describing the dynamics of a continuum of cell clones structured by the self-renewal potential. The model is an extension of the finite multi-compartment models of interactions between normal and cancer cells in acute leukemias. It takes a form of a system of integro-differential equations with a nonlinear and nonlocal coupling which describes regulatory feedback loops of cell proliferation and differentiation. We show that this coupling leads to mass concentration in points corresponding to the maxima of the self-renewal potential and the solutions of the model tend asymptotically to Dirac measures multiplied by positive constants. Furthermore, using a Lyapunov function constructed for the finite dimensional counterpart of the model, we prove that the total mass of the solution converges to a globally stable equilibrium. Additionally, we show stability of the model in the space of positive Radon measures equipped with the flat metric (bounded Lipschitz distance). Analytical results are illustrated by numerical simulations.


Introduction
This paper is devoted to the analysis of a structured population model describing clonal evolution of acute leukemias. Leukemia is a disease of the blood production system leading to an extensive expansion of malignant cells that are non-functional and cause an impairment of blood regeneration. Recent experimental evidence indicates that cancer cell populations are composed of multiple clones consisting of genetically identical cells (Ding et al. 2012) and maintained by cells with stem-like properties (Bonnet and Dick 1997;Hope et al. 2004). Many authors have provided evidence for heterogeneity of leukemic stem cells (LSC) attempting to identify their characteristics; for review see Lutz et al. (2012). Heterogeneity is further supported by the results of gene sequencing studies (Ding et al. 2012;Ley et al. 2008). However, it was shown in these studies that a limited number of clones contribute to the total leukemic cell mass. At most 4 contributing clones were detected in the case of acute myeloid leukemia (AML) and at most 10 in the case of acute lymphoblastic leukemia (ALL) (Ding et al. 2012;Lutz et al. 2012). Moreover, in most cases of ALL, the clones dominating the relapse have already been present at the diagnosis but undetectable by the routine methods (Van Delft et al. 2011;Choi et al. 2007;Lutz et al. 2013). Due to a quiescence, a very slow cycling or other intrinsic mechanisms (Lutz et al. 2013;Choi et al. 2007), these clones may survive chemotherapy and eventually expand (Lutz et al. 2013;Choi et al. 2007). This implies that the main mechanism of relapse in ALL might be selection of existing clones and not acquisition of therapy-specific mutations (Choi et al. 2007). Similar mechanisms have been described in AML (Ding et al. 2012;Jan and Majeti 2013). Based on these findings the evolution of malignant cells can be interpreted as a selection process for properties that enable cells to survive the treatment and to expand efficiently. The mechanisms of the underlying process and its impacts on the disease dynamics and on the response of cancer cells to chemotherapy are not understood. Gene sequencing studies allow deciphering the genetic relations among different clones; nevertheless the impact of many detected mutations on cell behaviour remains unclear (Ding et al. 2012). The multifactorial nature of the underlying processes severely limits the intuitive interpretation of the experimental data.
To investigate the impact of cell properties on the multi-clonal composition of leukemias and to elucidate the possible mechanisms of the clonal selection suggested by the experimental data, a multi-compartmental model was proposed and studied numerically in Stiehl et al. (2014). It assumes the form of the following system of ordinary differential equations,  l 1 1 (t) = 2a l 1 s(t) − 1 p l 1 l 1 1 (t), d dt l 1 2 (t) = 2 1 − a l 1 s(t) p l 1 l 1 1 (t) − d l 1 2 l 1 2 (t), . . . . . . . . . d dt l n 1 (t) = 2a l n s(t) − 1 p l n l n 1 (t), d dt l n 2 (t) = 2 1 − a l n s(t) p l n l n 1 (t) − d l n 2 l n 2 (t), with nonnegative initial data. The model describes time dynamics of a healthy cell line, denoted by c j , j = 1, 2 and of n clones of leukemic cells l i j , for j = 1, 2, and i = 1, . . . , n, at time t. Each population consists of two different cell types, proliferating and non-proliferating, denoted by j = 1 and j = 2, respectively. This two-compartment model is a simplification of the more realistic model with multiple differentiation stages; see Marciniak-Czochra et al. (2009), Stiehl et al. (2013) for an introduction to the model and its application to the healthy hematopoiesis; Getto et al. (2013), Nakata et al. (2011), Stiehl and Marciniak-Czochra (2011) for its analysis; and Doumic et al. (2011) for a continuousstructure extension. This model can be viewed as a structured population model with a discrete structure describing two differentiation stages and n + 1 cell types.
Parameters p c > 0 and p l i > 0 denote the proliferation rate of the healthy cells and the cells in the leukemic clone i, respectively, and a c and a l i are the corresponding maximal fractions of self-renewal, which depend on the proportion of symmetric and asymmetric cell divisions in the respective population. More precisely, the self-renewal fractions 0 < a c < 1 and 0 < a l i < 1 are the fractions of the progeny cells that remain in the compartment of proliferating cells. Consequently, (1 − a c ) and (1 − a l i ) are fractions of the dividing cells that differentiate and become non-proliferating. By d c 2 > 0 and d l i 2 > 0 we denote the clearance rate of the non-proliferating healthy cells and the cells in the ith leukemic clone, respectively.
The model is based on the assumption that leukemic clones and their normal counterparts respond to a hematopoietic feedback signalling and compete for signalling factors (cytokines). We assume that the feedback signal, s(t), decreases if the number of non-proliferating cells increases. Derivation of such nonlinear feedback loop was proposed in Marciniak-Czochra et al. (2009). It is based on a Tikhonov-type quasistationary approximation of dynamics of the extracellular signalling molecules, such as the G-CSF cytokine, which are secreted by specialised cells at a constant rate and degraded by a receptor-mediated endocytosis. Following the evidence from clinical trials that the mature granulocytes mediate clearance of G-CSF (Layton et al. 1989), we assume that dynamics of the signalling molecules depends on the number of non-proliferating cells. This assumption has been also supported by studies of receptor expression showing that the mature cells express significantly more receptors than the cells in bone marrow (Shinjo et al. 1997). Taking into account these observations, we obtain a model with a nonlinear coupling depending on the level of non-proliferating cells.
Numerical simulations of model (1) suggest that cells with a superior self-renewal potential, i.e. a maximum value of the parameter a, reflecting the probability that a daughter cell has the same properties and fate as its parent cell, have an advantage in comparison to their competitors, which leads to the expansion of this cell subpopulation (Stiehl et al. 2014). The phenomenon was shown analytically solely in the case of two competing populations, a healthy and a cancerous cell line (Stiehl and Marciniak-Czochra 2012).
To elucidate further mechanisms of clonal selection, we propose an infinitely dimensional extension of the multi-compartment model (1). We introduce a continuous variable x ∈ Ω that represents the expression level of genes (yielding a phenotype) influencing self-renewal properties of the cells. It leads to a system of integro-differential equations describing dynamics of a structured population with the continuum of cell clones and the two-compartment differentiation structure. Cells in Population 1 (dividing cells) proliferate and may self-renew or differentiate into Population 2 cells (differentiated cells). Population 2 cells do not proliferate and die after an exponentially distributed lifetime, as depicted in Fig. 1. Cells in both populations are stratified by a structure variable x. We assume that the self-renewal parameter  (2), consisting of two compartments corresponding to undifferentiated cells (dividing cells) and mature cells (differentiated cells). Undifferentiated cells (stem cells and early progenitors) divide symmetrically or asymmetrically. Accordingly, they produce cells of the same type (self-renewal) and mature cells (differentiation). Mature cells do not divide and they die after an exponentially distributed lifetime. The cells in each compartment are heterogenous. They are stratified by a structure variable x that represents the expression level of genes (yielding a phenotype and eg. influencing the self-renewal properties of the cells). Self-renewal and differentiation of cells are regulated by a cytokine feedback which, in turn, depends on the total count of differentiated cells depends on x, i.e. the parameter a becomes a function a(x). These assumptions lead to the model ∂ ∂t Assuming s(t) = 1/(1 + K Ω u 2 (t, x)dx), we obtain a nonlocal and nonlinear coupling of the two equations.
Our approach is motivated by the theory of selection of the most fit variants in adaptive evolution. Cells with different mutational variants might have different growth properties allowing them to expand more efficiently. The phenomenon can be understood as an example of a process, which is closely related to Darwinian evolution. In our particular case, certain rare mutants may have positive growth rates and be selected in environments that otherwise result in extinction. In other words, cells with a fitness advantage expand and dominate dynamics of the population leading to extinction of the other cell clones. The model proposed belongs to the class of selection models exhibiting a mass concentration effect, similar to those presented in the books Perthame (2007) and Bürger (2000).
In the current work, we do not model mutation events. Instead, motivated by the experimental findings described earlier in Lutz et al. (2013), Choi et al. (2007), we aim to understand which aspects of the dynamics of leukemias can be explained by the selection alone. It is interesting, since the relapse caused by an expansion of a clone that could not be detected at diagnosis due to the limited sensitivity of detection methods, can be misinterpreted as a mutational event (Choi et al. 2007). A computational model of the AML with mutations was proposed in Stiehl et al. (2014). Following the biological evidence (Jan et al. 2012), it was assumed that new LSC clones were formed due to mutations occurring in LSCs or due to the influx from the so-called preleukemic cells at a rate modelled by a time inhomogeneous Poisson process. At each point of the Poisson process a new clone with random cell properties was added to the system. Simulations of that model demonstrate that leukemic cell properties at diagnosis and at relapse are comparable to the scenario without mutations. Introducing mutations to the continuous models is known to make asymptotic analysis more complicated, and therefore we do not consider this aspect in the current paper.
The mathematical angle of our study is analysis of the nonlocal effects and development of singularities in the solutions of the integro-differential equations. We show that the solutions of system (2) may tend to Dirac measures concentrated in points with the largest value of the self-renewal potential. Such dynamics can be interpreted in the terms of selection, which causes convergence of the heterogeneous initial data to a stationary solution with the mass localised on a set of measure zero. Convergence then holds in the weak * topology of Radon measures. Considering the space of positive Radon measures with a suitable metric allows formulating the result on convergence of solutions to a stationary measure in the terms of the metric instead of the weak * convergence of Radon measures. We apply the flat metric (bounded Lipschitz distance), which has proven to be useful in the analysis of a variety of transport equations models, for example to study Lipschitz dependence of solutions of nonlinear structured population models on the model parameters and initial data ; see Appendix for the definition and properties of the flat metric. Similar results have been recently shown for scalar equations including diffusion; see for instance Barles and Perthame (2008), Barles et al. (2009), Lorz et al. (2011, 2013, Desvillettes et al. (2008), and  for a model with an additional space structure. The equations studied in  and  have been also applied to address cancer heterogeneity, and the influence of the selection process on the cancer resistance to chemotherapy.
The novelty of our work lies in considering a system of two coupled equations. Difficulty of the analysis is related to the specific nonlinearities in the model, which do not allow for component-wise estimates. The proof of boundedness of mass in the scalar equations is based on existence of sub-and supersolutions. In the case of a system, we face a difficulty which appears already in the proof of boundedness of solutions of a structure-independent model. The estimates cannot be concluded directly from the equations. To tackle this problem, we investigate the dynamics of the quotients of solutions of the two variables. Systems of equations also cause additional difficulties when analysing the long-term dynamics in comparison to the scalar equations due to the lack of a rich class of entropies. Convergence to a stationary positive Radon measure has been previously studied for a scalar integro-differential equation which is linear in the nonlocal term as in Jabin and Raoul (2011). This is often referred to as the Evolutionarily Stable Distribution. To deal with model nonlinearities, we make use of a Lyapunov function established previously for a finite dimensional counterpart of the model in Getto et al. (2013) and we show that the total masses of solutions tend asymptotically to the same equilibria.
A system of two equations describing selection and mutation in a stage-structured population has been investigated in Calsina and Cuadrado (2004) and Calsina and Cuadrado (2005) in the context of adaptive dynamics. Analysis of that model is based on a specific structure of nonlinearities appearing only in the mortality terms. Using irreducibility of the mutation operator and the infinite dimensional version of the Perron-Frobenius Theorem, it has been shown that solutions of the model converge to a stationary distribution, which concentrates at the point of maximum fitness in the case of the frequency of mutations tending to zero. The nonlinearity in our model is related to the growth term, which requires a different approach to the analysis of the asymptotic behaviour of the model solutions. The difference in the structure of nonlinear feedbacks is related to a different biological definition of the described processes. While the classical juvenile-adult dynamics is based on a loop of two positive feedbacks and no self-enhancement, the model of cell differentiation involves a negative feedback and a self-enhancement of the first population. Interestingly, the two-stage structure in our model yields stabilisation of the total populations, while even in the basic juvenile-adult models, the two-stage structure may lead to multiple attractors and limit cycles; see for example Baer et al. (2006). The paper is organised as follows: in Sect. 2, the main results are stated. Analytical results are illustrated by numerical simulations. Proofs of boundedness and strict positivity of the total masses and of the exponential decay of the model solutions outside the set corresponding to the maximal value of the self-renewal parameter are presented in Sect. 3. Section 4 contains the proof of mass convergence to a globally stable equilibrium. Finally, the asymptotic dynamics of the model solutions is shown in Sect. 5. Additionally, in Sect. 6, we show how to extend the analysis of our model to the framework of positive Radon measures with a suitable metric. Finally, in Section 7 we discuss biological conclusions and ideas stemming from this work. A summary of properties of the metrics used in Sect. 5 is provided in the Appendix.

Main results
We consider the following system of integro-differential equations and Ω ⊂ R is open and bounded.
In the remainder of this work we make the following assumptions on the model parameters and initial data.
either consists of a single point or it is a set with a positive Lebesgue measure.

Remark 1
The assumption (iv) on the self-renewal fraction a(x) is made to streamline the presented analysis. If Ω a consists of several isolated points, then the solution is attracted by a finite dimensional subspace spanned by Dirac deltas located at the maximum points of a; see Fig. 3. However, in this case the exact pattern may also depend on the shape of function a(x) near its maximal points. Since analysis of this case requires stronger assumptions on regularity of the initial data and the function a(x), we consider it separately in Theorem 3.

Existence and uniqueness of a classical solution
follow by the standard theory of ordinary differential equations in Banach spaces. More delicate is the question of asymptotic behaviour of the solutions of system (3). Our goal is to show that the solution u tends asymptotically to a stationary measure, as it is observed in the numerical simulations, see Figs. 2 and 3. The phenomenon is characterised by the following Theorem.
Theorem 1 Let Assumptions 1 hold and let (u 1 , u 2 ) be a solution of system (3) with initial data (u 0 1 , u 0 2 ). Then, u 1 and u 2 converge to stationary measures with supports contained in the set Ω a defined in expression (4), as t tends to infinity. Moreover, 2 , then the solution converges to a stationary measure (Dirac measure multiplied by a positive

Convergence holds in the flat metric (bounded Lipschitz distance); see Appendix for the definition and properties of the bounded Lipschitz distance. (ii)
If Ω a is a set with positive measure andā = max x∈Ω a(x) > 1 2 , then the solution converges to a stationary L 1 -function, such that Remark 2 If a(x) ≤ 1 2 for some points x ∈ Ω, then the solutions of the model converge point-wise to zero, i.e. lim t→∞ (u 1 (t, x), u 2 (t, x)) = (0, 0) for every x ∈ Ω − := {x ∈ Ω a(x) ≤ 1 2 . This is a straightforward consequence of Eq. (3), since ρ 2 is strictly positive, as shown in Lemma 1, and hence 2a(x) 1+Kρ 2 (t) − 1 < 0 for x ∈ Ω − . Therefore, we are interested in evolution of the system for x ∈ Ω + := Ω\Ω − . Subpopulations with a(x) ≤ 1 2 may affect short-term dynamics of the system; however they have no influence on the asymptotic behaviour.
Details of the proof are presented in Sects. 3, 4 and 5. The proof is based on the following key steps: Step 1. Uniform boundedness and strict positivity of masses ρ i (t) = Ω u i (t, x) dx for i = 1, 2 (Lemma 1).
Proof of this lemma is deferred to Sect. 3.1. Parameters used in the simulation: K = 0.01, p = 1, d = 0.2 and the initial data: We observe mass concentration in the point x = arg max x∈Ω a(x) and convergence of the mass to a stable stationary value  (3) with the self-renewal function a(x) having two equal local maxima (shown in the upper panel) and the parameters the same as in Fig. 2. We observe mass concentration in two points corresponding to the maximum of the function a(x) with unequal distribution of the mass between the two points Step 2. Exponential extinction of solutions in points outside the set Ω a (Lemma 3).
We start with characterising the asymptotic behaviour of the ratios of solutions taken at different x points.
a.e. with respect to the Lebesgue measure.
The proof of this lemma is deferred to Sect. 3.2. Lemma 2 yields the following result: As a consequence of Lemma 2 we also obtain ∈ Ω a a.e. with respect to the Lebesque measure.
The corresponding proof is presented in Sect. 3.2.
Step 3. Convergence of solutions to stationary measures. Convergence to the stationary solutions follows from the property of the total masses of the solutions ( Ω u 1 (t, x)dx, Ω u 2 (t, x)dx). We show that ifā = max x∈Ω a(x) > 1 2 , then the solutions converge to the stationary state of the system withā = max x∈Ω a(x).

Corollary 2 Total masses converge to the valuesρ
Details of the proof of mass convergence are deferred to Sect. 4.
If Ω a consists of a single pointx andā = max x∈Ω a(x) > 1 2 , then the exponential decay of the solutions outside the set Ω a together with the convergence of total masses, yields convergence of the solutions to a stationary measure concentrated at x (a Dirac measure multiplied by a positive constant). In the case of Ω a having a positive Lebesgue measure, convergence of solutions together with Corollary 1 on the stationary distribution of masses among different domain points yields convergence of solutions to the stationary equilibrium. Further details of the proof of convergence of solutions to the stationary measures are given in Sect. 5.
Remark 3 In the case Ω a = {x}, the convergence holds in the weak * topology of Radon measures. In general, we cannot expect the strong (norm-total variation) convergence of the solution to a stationary solution. If the set Ω a ⊂ R has zero Lebesgue measure and consists of a single point [compare Assumptions 1 (iv)], then the model solutions for any finite time point are uniformly continuous with respect to the Lebesgue measure and u i (t, ·)L 1 → c i δx , weakly * , for i = 1, 2. Here, u i (t, ·)L 1 denotes the measure such that u is its Radon-Nikodym derivative with respect to L 1 .
Hence, the distance between the two solutions T V (u i (t, ·), c i δx ) ≥ 2c i . The problem can be solved by considering convergence with respect to a suitable metric, for example the flat metric (bounded Lipschitz distance); for details see Sect. 5.
If the support ofā is not a single point set, then the stationary distribution of masses depends on the initial conditions. If Ω a has a positive Lebesgue measure, then the distribution of masses results from Corollary 1. If Ω a consists of a discrete set of points, then the stationary solution takes the form of a linear combination of Dirac deltas; see Fig. 3. We show that in such case the limit function depends on the shape of a(x) in the neighbourhood of the concentration points.
Theorem 3 (Co-existence of different stationary solutions) Let Assumptions 1 (i)-(iii) hold and, additionally, the initial functions u 0 1 , u 0 2 ∈ C(Ω). Let the set Ω a of the maximum values of the self-renewal parameter a (as defined in expression (4)) consist of two points Ω a = {x 1 ,x 2 } and u 0 1 be strictly positive on Ω a . Then, The proof of this theorem is deferred to Sect. 5.
Remark 4 If a is an analytic function and Ω ⊂ R, then a diffeomorphism satisfying condition (6) exists if the first nonconstant nonzero terms of Taylor expansion of the function a(x) are of the same order.
This observation suggests how to construct a(x) with Ω a = {x 1 ,x 2 } such that solutions extinct at one of the points of Ω a . For example, we may define a(x) with f or x ∈ ( 5 8 , 1]. and a smooth extension of a(x) on the interval ( 3 8 , 5 8 ) satisfying 0 < a(x) < 1. We obtain Ω a = { 1 4 , 3 4 }, and a mapping Φ(x) = x − 1 4 + 3 4 satisfying condition (6) on and it is singular at x = 3 4 . Hence, the total mass concentrates at the point x = 3 4 and there is an extinction of mass at x = 1 4 .

Boundedness and strict positivity of masses
All considerations in this Section hold for x ∈ Ω a.e. with respect to the Lebesque measure.

Lemma 4 Under the assumptions of Lemma 1, the function U
Since > 1 −ā, and the right-hand side of Eq. (7) is a logistic type nonlinearity, we conclude that By definition of U , we can infer that As a straightforward consequence of Lemma 4, we deduce Corollary 3 Under the assumptions of Lemma 1, it holds Now we state another technical result in the spirit of Lemma 4.

Lemma 5 There exist constants M
Proof Calculating the derivative of the quotient of ρ 2 (t) and ρ γ 1 (t), we obtain d dt This estimate holds for arbitrary γ ∈ (0, 1), so in particular for those satisfying γ p − d < 0. Arguing as before, we deduce that, for all t ≥ 0, Equipped with Lemmas 4 and 5, we prove Lemma 1.
To show boundedness of ρ 1 , we apply inequality (8) to the first equation of system (3) Integrating this inequality over Ω yields Using a similar argument as in the proof of Lemma 4, we conclude that Boundedness of ρ 2 results from the second equation of system (3), nonnegativity of ρ 2 and the assumptions on a. It holds Integrating over Ω and using (11), we obtain Hence, we conclude that (ii) We show that masses ρ 1 and ρ 2 have a strictly positive lower bound, uniform in time.
We estimate the growth of ρ 1 using a decomposition of the domain Ω = Ω − +Ω + , where Ω − := {x ∈ Ω a(x) ≤ 1 2 and Ω + := {x ∈ Ω a(x) > 1 2 . First, we assume that the set Ω − is nonempty, i.e. Ω − u 0 1 (x) > 0. We denote Using the explicit form of the solution 1+Kρ 2 (τ ) −1 p dτ (13) and the properties of the function a(x) on the two subdomains, we obtain Combining estimates (9) and (14) yields With estimate (15) at hand, we show that ρ + 1 is strictly positive for every t ∈ R + . We estimate its dynamics The term in the brackets is strictly positive for ρ + 1 small enough, i.e. for which is a positive constant, since a > 1 2 . Hence, we obtain the estimate Consequently, we obtain the strict positivity of ρ 1 and using the second equation of (3), also the strict positivity of ρ 2 . In the case of Ω − = ∅, it holds ρ 1 = ρ + 1 and the proof is complete if we set M 5 = M 4 .

Asymptotic behaviour of the solutions
In the next step, we show that the first component of the solution of system (3) tends to zero forx / ∈ Ω a a.e. with respect to the Lebesque measure.
Proof (of Lemma 2) We choose two points x 1 , x 2 ∈ Ω such that a(x 1 ) − a(x 2 ) < 0, and calculate ∂ ∂t Solving the above differential inequality for u 1 (t,x 1 ) u 1 (t,x 2 ) , we obtain the assertion of this Lemma by the choice of x 1 and x 2 .
a.e. with respect to the Lebesque measure.
Proof We use a similar ansatz as in Lemma 2 and calculate for t > 0 ∂ ∂t Applying Lemma 2, we obtain ∂ ∂t Thus, we deduce the following bound for u 2 (t,x 1 ) Having shown the dynamics of the ratios of the values of a solution at different x points, we prove that the solutions converge to zero outside the set of points with a maximum value of the parameter a(x).

Proof (of Lemma 3) Letx be a point different fromx and assume that lim t→∞ u(t,x)
> 0. Continuity of a(x) implies that the set of x, such that a(x) > a(x), is an open nonempty set and, therefore, it has positive measure. Since Lemma 2 holds for every x,x ∈ Ω such that a(x) − a(x) > 0, we conclude that u(t, x) tends exponentially to +∞ for every x such that a(x) > a(x). This is, however, in contradiction with the uniform boundedness of the mass Ω u(t, x)dx.

Proof of convergence of the total mass
We begin the proof of Theorem 2 by showing the following lemma, which allows comparing two dynamical systems.

Lemma 7 Let t → X F (t, ·) be a one-parameter family of C
, with a single minimumū, is a strict Lyapunov functional, i.e. d dt X F (t, u)| t=0 ·∇V (u) = 0 for u =ū and d dt X F (t, u)| t=0 ·∇V (u) < 0 otherwise. Then, ifũ is a solution of where lim t→∞ sup τ ∈[t,∞) | f (τ )| = 0 and Im(ũ(·)) : Proof For arbitrary a >V , we define a truncation Since V a ∈ W 1,∞ (U ), where U is the intersection of all convex sets containing Im(ũ(·)), U = conv(Im(ũ(·))) ⊂ (0, ∞) × (0, ∞), and d dtũ ∈ L 1 (Ω), then we can define the time derivative of V a (ũ(t)) using the chain rule. ∇ u V a is defined in a classical sense only outside the set V (u) = a, but it has a Clarke derivative, i.e. a generalised subdifferential for a locally Lipschitz function (Clarke 1983), on the set V = a. In the following, ∇ũ V a (ũ) is an extension of the classical definition, involving the maximal element of the Clarke derivative, to the set where the classical derivative is not defined.
Let us define β : ImV (u) → (0, ∞) such that Since β is a continuous function defined on a compact set, it achieves a strictly positive minimum. Furthermore, for the truncation function V a , there exists a positive constantβ a such that β(V a ) ≥β a V a . Hence, we obtain Using compactness of the set U , we estimate ∇ũ V a (ũ(t)) by its L ∞ norm, which yields the following inequality, Integrating the above estimate, we obtain We show that the right-hand side of inequality (19) tends to zero for t → ∞.

Proof (of Theorem 2)
To apply Lemma 7 to system (3), we consider a finite dimensional model obtained by setting a(x) to a constant valueā Note that the above equation generates a C 1 -semiflow, which is invariant on (0, ∞) × (0, ∞). We check that the two systems (20) and (25) fulfill the assumptions of Lemma 7. Lyapunov function for system (20) has been previously constructed in Getto et al. (2013). It assumes the form where (v 1 ,v 2 ) is the stationary solution, and Lyapunov function (21) is well-defined for every Note that V 1 (v 1 ) is strictly convex and therefore ∂ ∂v 1 V 1 = 0 for v 1 =v 1 . Similar observation holds for V 2 (v 2 ). Hence (v 1 ,v 2 ) is the global minimum of the Lyapunov function.
Direct calculations, as provided in Getto et al. (2013), allow to check that for the solutions of system (20). Moreover, the equality d dt V (v 1 (t), v 2 (t)) = 0 holds only for the stationary solution (v 1 ,v 2 ).
To show that the perturbation function on the right-hand side converges to zero as t → ∞, we calculate where Ω a is defined in the expression (4). Consequently, using boundedness of ρ 1 , boundedness of a(x) as well as Lemma 3, we obtain that and hence we conclude that system (25) fulfills the assumptions of Lemma 7. Consequently, we obtain that the total mass of a solution of system (3) converges to a globally stable equilibrium, which is equal to the equilibrium of the ordinary differential equations model (20) corresponding to the maximum value of the self-renewal parameterā. Thus, we have proven the assertion of Theorem 2. Fig. 4 Using the trapezoid rule to approximate the integral of ρ 1 (t), ρ 2 (t), we observe numerically the convergence of the total mass to a constant value. The parameter set is the same as in Fig. 2

Proof of the convergence result
Finally, we obtain the main assertion.
Proof (of Theorem 1) Lemma 3 implies that the solutions of system (3) decay exponentially to zero in all points x / ∈ Ω a . We consider two cases (compare Assumptions 1 (iv)): Convergence to a stationary solution follows from the convergence of mass given by Theorem 2. Hence, the solutions converge to measures concentrated atx: where L 1 denotes a one dimensional Lebesgue measure and u i (t, ·)L 1 is the measure which Radon-Nikodym derivative with respect to L 1 is equal to u, δx is a Dirac measure localised atx and c i , i = 1, 2, are the stationary masses, i.e. c 1 =ρ 1 = d p 2ā−1 K and c 2 =ρ 2 = 2ā−1 K . The convergence result can be understood in a suitable metric on the space of positive Radon measures. We apply here the flat metric ρ F , also known as the bounded Lipschitz distance (Neunzert 1981). For completeness of presentation, the definition and basic properties of this metric are provided in Appendix.
To estimate the distance between a solution u i (t, ·) and the stationary measure c i δx , i = 1, 2, we use the following inequality for the distance of two measures μ and ν For the proof of this inequality we refer to  and Jabłoński and Marciniak-Czochra (2013). Here W 1 μ μ(Ω) , ν ν(Ω) denotes the Wasserstein distance between two probabilistic measures; see Appendix for the definition of the Wasserstein metric. We calculate, for i = 1, 2, The first term on the right hand-side of inequality (27) can be estimated using the exponential estimates of Lemma 2. To show that it converges to zero we apply the Kantorovich-Rubinstein Theorem (Villani 2003(Villani , 2006 and use the equivalent definition of the Wasserstein metric given as the cost of optimal transport with the cost function |x − y|, i.e. where γ ∈ Γ μ μ(Ω) , ν ν(Ω) is a joint distribution (probabilistic measure) with the marginal distributions μ μ(Ω) and ν ν(Ω) , and where is the family of all joint distributions with marginal distributions μ μ(Ω) and ν ν(Ω) . We estimate the difference between a normalised solution π i (t) := u i (t, ·) ρ i (t) L 1 and its limit δx , i = 1, 2. Using a joint distribution γ i = δx ⊗ π i , i = 1, 2, we obtain To show that the right-hand side of inequality (29) converges to zero, we define a set Ω a−ε = {x : a(x) >ā − ε}. For ε small enough, there existsε > 0 such that the set Ω a−ε is contained in aε−neighbourhood of Ω a , i.e. Ω a−ε ∈ [x −ε,x +ε]. By Lemma 2, π i (t) (Ω\[x −ε,x +ε]) → 0 for t → ∞. Therefore, we obtain x +ε]) +ε →ε, for t → ∞.
Since the above convergence holds for anyε > 0, we conclude that Convergence of the second term in formula (27) is due to Theorem 2. Hence, we obtain that If Ω a is a set with positive measure, no singularities emerge due to the uniform boundedness of the total mass. In this case, the solution tends to zero outside Ω a and to a positive L 1 -function on Ω a . Following Corollary 1, we conclude that the exact shape of the limit solution depends on the initial distribution. (iii) Ifā = max x∈Ω a(x) ≤ 1 2 , then the solutions converge exponentially to zero, what is a consequence of Eq. (3). We estimate where C = − 1 1+K min t∈[0,∞) ρ 2 (t) − 1 p > 0, due to Lemma 1. Hence, using the Gronwall inequality, we obtain the exponential decay to zero. Finally, convergence ρ 2 (t) → 0 as t → ∞ follows from the estimate Since the solutions (u 1 , u 2 ) are nonnegative, they converge to zero in L 1 (Ω).
Finally, we analyse the case with Ω a consisting of two points and prove the coexistence and the extinction result.
Proof (of Theorem 3) (i) We investigate dynamics of the mass of a solution of system (3) around the points of Ω a . Let us assume that there exists a diffeomorphism Φ ∈ C 1 (U 1 ), where U 1 is an open neighbourhood ofx 1 , such that Φ(x 1 ) =x 2 and a(x) = a(Φ(x)) for all x ∈ U 1 . Using the explicit form of the solution (13) and the property Φ(x 1 ) =x 2 , we obtain Changing variables on the right hand-side of (30) leads to where J Φ is Jacobian of the diffeomorphism Φ.
does not depend on time and is continuous with respect to y and since u(t, x) converges pointwise to zero outside Ω a = {x 1 ,x 2 } (see Lemma 3), we obtain Hence, the solution converges to a measure c 1,1 δx 1 + c 1,2 δx 2 with strictly positive c 1,1 and c 1,2 such that c 1,1 Since the total mass of u 1 is equal to c 1,1 + c 1,2 =ρ 1 , whereρ 1 is given in Corollary 2, the constants c 1,1 and c 1,2 are uniquely determined. Relationship (33) indicates that the mass distribution between the different concentration points depends on the shape of the function a(x) and on the initial data.

Remark 5 Continuity of
J Φ −1 (y) requires continuity of the initial data and strict positivity of u 0 1 on Ω a , which is reflected in the stronger assumptions of the theorem compared to Assumption 1.

Extension to initial data in the space of Radon measures
The phenomenon of mass concentration provides a motivation to consider the model in the space of positive Radon measures, as defined by the following equations with the initial data where μ 0 i are nonnegative Radon measures for i = 1, 2. x ∈ Ω ⊂ R n , for some n ≥ 1, denotes the state of a cell and, for every Selection-mutation models in the spaces of positive Radon measures have been studied by many authors Ackleh et al. (1999), Ackleh et al. (2005), Bürger and Bomze (1996), Bürger (2000), Caizo et al. (2013), Cleveland and Ackleh (2013), Desvillettes et al. (2008). In this context, convergence of the solutions with respect to the Prokhorov metric has been considered in Ackleh et al. (1999). For the relation between the Prokhorov metric and the Wasserstein distance used in our paper we refer to Gibbs and Su (2017).
Steps of the proof of Theorem 1 can be repeated for the measure-valued solutions with some modifications of the lemmas which rely on point-wise estimates of the quotients of solutions. Assuming that the initial data are measures such that μ 0 1 is absolutely continuous with respect to μ 0 2 , Lemma 4 can be reformulated for the model (34)-(36) by considering a Radon-Nikodym derivative instead of the point-wise quotients.
Next technical difficulty appears in Lemma 2. To show the asymptotic behaviour of the measure-valued solutions, we can apply the framework developed in Bürger and Bomze (1996). In the remainder of this section, we briefly discuss this extension.
The first equation of the model (34)-(36) can be re-defined in the terms of a probabilistic measure modelling the frequency of a certain phenotype x ∈ B in the population of mitotic cells μ 1 . It is given by the quotient where B ⊂ Ω is a Borel set, as defined before. Using the equation for μ 1 , we obtain The model can be then formulated in the framework presented in the book by Bürger (Bürger 2000). Denoting the mean fitness by and the multiplication operator A(t) by we rewrite Eq. (38) as an ordinary differential equation in the space of Radon measures However, the obtained equation is more general than the abstract equation in Bürger (2000), due to the dependence of A on time. Nevertheless, it holds

A(t) = (A(t)π(t)) (Ω).
Using the form of the operator (40), we rewrite it as a function of time α(t) = 2 p 1+ρ 2 (t) multiplied by a time independent operator (Aπ(t)) (B) = B a(x)π(t)(dx), This structure allows to follow the lines of Bürger and Bomze (1996) and focus on a differential equation given by The structure assures that the family of operators A commutes. The operator A is bounded and it generates a positive semigroup on the space of positive Radon measures M + (Ω).
Since α is a strictly positive and bounded function, due to the properties of ρ 2 shown in Lemma 1, we can rescale time, s = t 0 α(ξ )dξ , and obtain a linear autonomous differential equation Equivalence to a linear differential equation yields convergence of solutions to a solution π(t) with the support concentrated on the set of maximal value of a(x), a = sup x∈Ω∩supp(μ 0 1 ) a(x). The latter result is the extension of our Lemma 2 to the measure-valued solutions.
In summary, by adapting the framework developed in Bürger and Bomze (1996), our results can be extended to the measure-valued solutions in the case of the model of the clonal evolution without mutations. Asymptotic analysis carried out in Bürger and Bomze (1996) is based on the application of the infinite-dimensional version of the Perron-Frobenius Theorem, which is possible in models with dynamics governed by an irreducible operator. The latter is the case in models involving mutations described by an integral operator satisfying irreducibility conditions. That approach cannot be, however, directly applied to the extension of our model to the case with mutations. The difficulty is related to the estimates for the time dependent operator A defined in expression (40), which rely on the equations for the ratios of solutions in Lemma 4, or Radon-Nikodym derivatives (37), which cannot be established in the model with an additional nonlocal mutation operator. Therefore, including mutations in our model requires a different proof of the uniform boundedness and strict positivity of ρ 2 and extension of the analysis to the model with mutations remains an open question.

Discussion
In this paper, a discrete multi-compartmental model of multiple cell lineages has been extended to a model coupling a two-stage differentiation structure with a continuous structure of phenotypes. The latter allows to investigate the role of the intra-cancer heterogeneity, including competition between healthy and cancer cells and dynamics of the multi-clonal structure of the system.
Based on recent analyses of the clones consisting of mutational variants in cancer (Miller et al. 2014), it follows that the dynamics of clone distributions may in many cases consist solely of change in relative frequencies of different clones. More specifically, the clones that have been dominant in the primary tumour, are out-competed by other clones in the relapsing or metastatic tumours, which had low frequencies in the primary. The model in this paper provides a "mechanistic" explanation for these observations, which is also mathematically rigorous.
Asymptotic analysis of the proposed system of integro-differential equations suggests that the selection process may be governed by the cell's property of self-renewal that determines the fitness of each clone and ultimately leads to survival or extinction.
Theorem 1 shows that, in a well-mixed cell production system, a negative nonlinear feedback such as that the one proposed in Lander (2009), Marciniak-Czochra et al. (2009, leads to the selection of the subpopulation with the superior self-renewal potential. The assumption that the cell population is well-mixed leads to the nonlocal effect and is modelled using the integral term. This assumption reflects well the structure of the hematopoietic system. Consequently, our results suggest that the greater clonal heterogeneity observed in solid cancers than in blood cancers may be due to spatial effects of the cell-to-cell interactions. Additionally, Theorem 3 suggests some explanation of the co-existence of different clones having the same fitness. The results stress the importance of self-renewal in cancer dynamics and allow concluding that slowly proliferating cancer cells with a high self-renewal potential are able to outcompete the cells that divide faster. It suggests an explanation of the clinical dynamics such as resistance to treatment. Importance of this observation in the context of the leukemia evolution, the response to chemotherapy and the dynamics of the disease relapses has been discussed in Stiehl et al. (2014). The results obtained provide an explanation of the observed clonal selection in the acute myeloid leukemia in the course of the disease development and the relapse after chemotherapy reported by Ding et al. (2012). Recently, fitting the AML model to patients' data has suggested that an increased self-renewal is correlated with a poor patient prognosis (Stiehl et al. 2015).
The ρ F distance metrizes both weak* and narrow topologies on each tight subset of Radon measures with uniformly bounded total variation (Schwartz 1973;Ambrosio et al. 2000).
Remark 6 Every bounded Radon measure on a bounded set Ω has an integrable first moment and hence the distance ρ F is finite.