Adaptive dynamics of hematopoietic stem cells and their supporting stroma: A model and mathematical analysis

We propose a mathematical model to describe the evolution of hematopoietic stem cells (HSCs) and stromal cells in considering the bi-directional interaction between them. Cancerous cells are also taken into account in our model. HSCs are structured by a continuous phenotype characterising the population heterogeneity in a way relevant to the question at stake while stromal cells are structured by another continuous phenotype representing their capacity of support to HSCs. We then analyse the model in the framework of adaptive dynamics. More precisely, we study single Dirac mass steady states, their linear stability and we investigate the role of parameters in the model on the nature of the evolutionary stationary distributions (ESDs) such as monomorphism, dimorphism and the uniqueness properties. We also study the dominant phenotypes by an asymptotic approach and we obtain the equation for dominant phenotypes. Numerical simulations are employed to illustrate our analytical results. In particular, we represent the case of the invasion of malignant cells as well as the case of co-existence of cancerous cells and healthy HSCs.


Introduction
Hematopoietic stem cells (HSCs), developing in the bone marrow, are immature cells that are (the earliest in development) precursors of all lineages of blood cells: red blood cells, white blood cells and megacaryocytes (whose fragmentation gives rise to platelets).Blood cell formation, also called hematopoiesis, is a complex phenomenon basing on the self-renewal, differentiation and maturation of HSCs.It produces about 10 11 blood cells per day in humans and is one of the most stable biological processes in vertebrate organisms.A dysfunction in the hematopoietic process may induce blood cancer diseases (usually named malignant hemopathies) such as leukaemia where blockade of maturation and of differentiation occurs in the hematopoietic tree.As a consequence, malignant cells, resulting from an accumulation of irregular genetic events, appear and proliferate abnormally.
Many mathematical models have been proposed to understand blood cell development and blood diseases.Mackey [18], inspired by Burns and Tannock [6] and Lajtha [15], have introduced a first mathematical model of the form of a system of delay differential equations for the dynamics of HSCs where the populations are divided into two groups (proliferating cells and quiescent cells) and the time delay corresponds to the proliferating phase duration.Further improvements both in modelling and mathematical analysis are investigated by many authors; see, for example, [2,3,19,27] -models in the form of ODEs or age-structured transport equations with applications to chronic myelogenous leukaemia, [12] -a diffusion model including spatial competition between cells -, reviews [1,8,26] and the references therein.
Despite extensive studies on the dynamics of HSCs and diseases of the hematopoietic system, none of the above-mentioned models takes into account the interactions between HSCs and the hematopoietic niche which is a specific microenvironment ensuring the maintenance and regulation of HCSs locally.It is worth noting that the interactions between HSCs and their niche, of which mesenchymal stem cells (MSCs) are the most important component, play a crucial role in the formation of mature blood cells.Also, alterations in the bidirectional exchanges between HSCs and MSCs may give rise among HSCs to blood cancer stem cells, i.e., leukemic stem cells, LSCs.
Note that healthy HSCs need the close presence of stromal cells for their development but stromal cells can proliferate without HSCs.Similar to HSCs, cancer cells in the early stages need stromal cells for their development whereas in the later stages they can proliferate without support cells.In other words, the more malignant a cell is, the more independent of stromal cells it is.Here, cancer cells in earlier stages are cells with few mutation events and cancer cells in later stages stand for the ones with more mutation events.We refer to, for example, [4,7,28] for reviews of the interaction between HSCs and stromal cells, [13] for acute myeloid leukemic cells.
In the present paper, we introduce a mathematical model for the interaction between HSCs and stromal cells with the aim to better understand the nature of the dialogue between them as well as their dynamics.We also perform a mathematical analysis for the long-time behaviour of the hepatopoietic and stromal cells in the framework of adaptive dynamics.Our mathematical model and some notions in the framework of adaptive dynamics, in particular, evolutionary stable distributions (ESDs) are given in the remaining part of this section.

A mathematical model
Let n h (t, x) be the population density of hematopoietic stem cells (HSCs) and cancer cells at time t with phenotype x, continuous structure variable assumed to characterise the population heterogeneity in a way relevant to the question at stake.Here x will represent a malignancy potential of HSCs, from its minimal (representing a totally healthy state) to its maximal value (representing maximum malignancy), considered independently of their stromal support.From a biological point of view, x might represent a pathological combination of both high plasticity (i.e., ability to change phenotype; stem cells are plastic, but physiologically, they not proliferate much) and fecundity (i.e., ability to proliferate; differentiated cells are able to proliferate, but physiologically, they show little plasticity).Let n s (t, y) be their corresponding stromal cell population density -that we will sometimes call support cells -at time t and with phenotype y (here the continuous phenotype variable y will denote the supporting capacity of MSCs to HSCs).Assume for simplicity that x and y are real variables with x ∈ (a, b), y ∈ (c, d), where 0 < a < b and 0 < c < d.Totally healthy HSCs will thus have a phenotype x close to a, while aggressive leukemic HSCs (i.e., LSCs) will have a phenotype x close to b.We consider a mathematical model of the form This system is completed with initial data Here our assumptions and notations are n s (t, y) dy are the total populations of HSCs and their support cells, respectively.
ψ s (y)n s (t, y) dy denote an assumed chemical signal (Σ h ) from the hematopoietic immature stem cells (HSCs) to their supporting stroma (MSCs), i.e., "call for support" and conversely, a trophic message (Σ s ) from MSCs to HSCs.The cytokine stem cell factor (SCF) and the C-X-C motif chemokine ligand 12 (CXCL12) are typical examples for such supporting messages [7,28].The nonnegative functions ψ h , ψ s defined on (a, b) and (c, d) measure the contribution of each phenotype in the interactive messages.Assume that ψ s ≥ 0 (the higher the support phenotype in MSCs, the stronger the trophic message to HSCs); in the same way, unless otherwise specified, we shall assume that ψ h ≥ 0.
• The term r h ≥ 0 represents the intrinsic (i.e., without contribution from trophic messages from MSCs, nor limitation by the non local logistic terms ρ h and ρ s , that represent competition for space and nutrients within the whole population of cells) proliferation rate of HSCs.Assume that r h is non-decreasing (the more malignant, the more proliferative), r h (a) = 0 and r h (b) > 0.
• The term α ≥ 0, satisfying α ≤ 0 and α(b) = 0, is the sensitivity of HSCs to the trophic messages from support cells • For the term r s ≥ 0, we assume that r s (y) ≤ 0 (there is a cost in proliferation for support cells to increase their support capacity).The term β(y) ≥ 0 with β (y) ≥ 0 represents the sensitivity of the stromal cells MSCs to the (call for support) message coming from HSCs.
System (1.1) falls within the broader class of models for interacting populations where competitive, preypredator and cooperative types are typical examples of such interaction; see, for example, [22,Chapter 3].Apart from the cases mentioned in [22,Chapter 3], in the context of adaptive dynamics, the populations are often structured by phenotypical traits to take into account the heterogeneity in the population (e.g.[14]).We refer to [24] a related competitive system with healthy and cancer cells structured by a phenotypic variable related with their resistance to chemotherapy, to [25] an integro-differential Lotka-Volterra system for the interaction of N populations (N ≥ 2).In our model, besides the competition terms between cells, we introduce new terms Σ h , Σ s , α, β to represent the interacting messages between HSCs and stromal cells.The presence of these terms makes the problem difficult to study since the nature of (1.1) is unknown and may vary in time.It could be competitive, co-operative or other types depending on the sign of the terms −ρ s (t) + α(x)Σ s (t) and −ρ h (t) + β(y)Σ h (t).Note that if Σ h = Σ s = 0, our model reduces to the cases studied in [14,24].Also when ψ h = ψ s = 1, our problem becomes a particular case of [25].
Let us briefly sum up the meaning of our assumptions.From a biological point of view, the healthy HSCs cannot proliferate without support cells while cancer cells persist even without support cells.In our model n h (t, a) corresponds to healthy HSCs and n h (t, b) are leukemic cells since the intrinsic proliferation rate r h satisfies r h (a) = 0, r h (b) > 0. The monotonicity of r h implies that the higher x is, the more aggressive is a HSC (in fact, a LSC, since it is then cancerous).Also the monotonicity of α (this function stands for the sensitivity of HSCs to the trophic messages coming from MSCs) indicates that the more aggressive x is, the less sensitivity has a HSC to the trophic message sent by MSCs (i.e., the more independent it is from the surrounding stroma).Moreover, the condition α(b) = 0 shows that n(t, b) (i.e., cancer cells in the latter stages) proliferate independently of the supporting stroma.Furthermore, the monotonicity of r s , β shows that the more supporting stromal capacity MSCs have, the less they proliferate and the more sensitive to messages from HSCs they are.
As a simple case, the parameters r h , α, r h , r h , ψ h , ψ s can be chosen as linear or quadratic functions.For example, r h , α are given by r A more general model which includes the possibility of mutations has the form: Here the diffusion terms represent mutation with rates µ h , µ s and c 11 , c 12 , c 22 , c 22 measure the strength of competition between cells.Problem (1.3) reduces to (1.1) by setting µ h = µ s = 0 and c 11 (x) = c 12 (x) = 1, c 21 (y) = c 22 (y) = 1.Thus (1.1) can be considered as a good approximation of (1.3) in the regime µ h , µ s << 1 which is realistic since mutations occur rarely in physiology (to fix ideas, let us say between once every 10 6 and 10 9 cell divisions; of course, in evolved cancers, such low rates may increase).

HSCs Cancer cells
Stromal cells Interacting messages (Σ h , Σ s ) are represented by black arrow lines.The sensitivities of cells are described by red curves: HSCs are very sensitive to interacting messages (α > 0) while cancer cells (at their later stages) are independent of the surrounding stroma (α = 0).Blue curves refer to the intrinsic proliferation rates: HSCs cannot survive without supporting messages (r h = 0) while cancer cells can proliferate without supporting messages (r h > 0).

Some notions in the theory of adaptive dynamics
Hereafter, we will use the notation to denote the integrals over  As our equation arises from biological cell population dynamics, we only consider non-negative steady states and ESDs in this paper.
In the context of adaptive dynamics, when (1.5) holds, we say that the phenotypes x ∈ I, y ∈ J are living in the stationary environment given by (ρ h , ρs , Σh , Σs ).The functions G h , G s are fitness functions associated with this environment.The quantities G h (x), G s (y) are growth rates of phenotypes x, y and they tell us whether (x, y) can invade this environment: If G h (x) > 0, G s (y) > 0, (x, y) can grow and the system will reach a new equilibrium.We refer to [9,20] for more details about the framework of adaptive dynamics.

Summary of main results and organization of the remaining part of the paper
We first present mathematical results to verify biological properties concerning the independence of stromal cells on HSCs and its vital support to HSCs in Section 2.1.A uniform bound in time and a well-posedness result are given in Section 2.2.
As generally only a finite number of traits is represented in the equilibrium, we study, for simplicity, in Section 3, equilibria with only one trait, i.e., those of the form of single Dirac mass.The linear stability result of single Dirac mass steady states (Theorem 3.2) exhibits the mechanisms by which another trait can or cannot invade the stationary state produced by a given trait.
ESDs-equilibria corresponding to optimal states of the evolution-are studied in Section 4. We study the impact of the parameters on the form of ESDs.Two cases are investigated: monomorphic situation (i.e., only one trait is represented in ESDs) and dimorphic one (two traits are represented in ESDs).More precisely, we provide sufficient conditions to guarantee that all ESDs are monomorphic or dimorphic (Proposition 4.1).Also in Theorem 4.2 we obtain a result on the uniqueness of ESDs and we show that this unique ESD is monomorphic.Another result on the uniqueness of ESD (Theorem 4.5) hold for rather than general functions r h , r s and under some homogeneity assumptions of stromal cells.This theorem is concerned with more general ESDs which are not necessary to be monomorphic or dimorphic.
Section 5 is concerned with dominant traits which are the best adapted ones to the environment and favored at high population densities.These traits are represented by maximum points of the population densities and change in time because of the variation of environment.We study the movement of these traits in a long time scale, hence make the change of variable, τ = εt, with small ε > 0. We represent the dynamics of dominant traits (Theorems 5.1 and 5.2) in the regime as ε → 0 (asymptotic analysis).Also we obtain the equation for dominant traits Eq. (5.11).
In Section 6, we provide numerical simulations to illustrate our results and finally, some discussions are given in Section 7.

Hematopoietic stem cells or stromal cells without mutual interaction
In the absence of stromal cells, the system (1.1) reduces to the equation Similarly, the behaviour of stromal cells without HSCs is given by (i) For (2.1), we have n h (t, x) → r h (b)δ {x=b} weakly in the sense of measures as t → ∞.
(ii) For (2.2), we have n h (t, y) → r s (c)δ {y=c} weakly in the sense of measures as t → ∞.
These results confirm biological properties mentioned in Section 1 about the bi-directional interaction between HSCs and stromal cells.More precisely, in the absence of stromal cells for HSCs or of HSCs for stromal cells, the phenotypes of n h and n s , respectively, behave as monomorphic.Moreover, Lemma 2.1 (i) implies that healthy hematopoietic cells eventually go extinct while cancer cells (with the phenotype y = b) are selected and persist.Furthermore, stromal cells persist without HSCs and stromal cells with the lowest capacity of support will be selected (cf.Lemma 2.1 (ii)).

A well-posedness result
For a function f defined on an interval I, we set In the results below, we only need the boundedness of r h , r s , α, β, ψ h , ψ s .We do not need the monotonicity assumption of these functions.
Then the solution of (1.1)-(1.2) is non-negative and satisfies Proof.First note that the solution of (1.1) can be written in the form where Thus n h (t, x) ≥ 0, n s (t, y) ≥ 0 since n h0 (x) ≥ 0, n s0 (y) ≥ 0. This yields the first inequality of (2.4).
Integrating the two equations in (1.1) yields Summing up these two identities and using the non-negativity of n h , n s , we obtain Hence the second inequality of (2.4) follows.
Proposition (2.2) and a standard argument imply the following result.

Steady states and linear stability
Hereafter, for simplicity, we assume that ψ h (x) = x, ψ s (y) = y.Then Problem (1.1) becomes (3.1) 3.1 Single Dirac mass steady states Consider a particular case where hematopoietic and support cells evolve as single Dirac masses concentrated at x, ŷ.In other words, we focus on the behaviour of the size of the populations.In this case, n h , n s have the form We can obtain an explicit form of single Dirac mass steady states.
Proof.Single Dirac mass steady state solution yields An elementary calculation shows that (ρ h , ρs ) is a steady state of (3.2) and the corresponding Jacobian matrix is given by In view of (3.4), tr(J) < 0, det(J) > 0 so that the two eigenvalues of J are negative.Thus the linear stability of (ρ h , ρs ) follows.

Linear stability of single Dirac mass steady states
The results below are concerned with the linear stability of single Dirac mass steady states among perturbations of particular forms.
Then (n h , ns ) is linearly stable among perturbations starting by Proof.We linearize the system at (n h , ns ).For g h (t, x) := n h (t, x) − nh (x), g s (t, y) := n s (t, y) − ns (y) we obtain Note that g h , g s have the form Therefore, we obtain The corresponding matrix is given by The eigenvalues of the above matrix are and the two eigenvalues of the matrix J in (3.6).All are negative due to (3.7) and Lemma 3.1.Thus the stability of (n h , ns ) follows.
Interpretations of Theorem 3.2: With the notation (1.4), in view of (3.5), we have G h (x) = 0, G s (ŷ) = 0. Also the conditions G h (x * ) < 0, G s (y * ) < 0 show that the phenotypes x * , y * cannot invade the stable equilibrium (n h , ns ).As a consequence, no new equilibrium can be reached but a mutant invading the resident population (x, ŷ).

Sufficient conditions for monomophic and dimorphic ESDs
Recall that from the defintion 1. ) for nh mean that the blue curve must be below the red line and that the pair (α(x), r h (x)) for all x ∈ I := supp nh are the coordinates of the intersection points between the blue curve and the red line.Similarly, we have the same illustration for ns .If α (resp.r h ) is strictly monotone and if r h (α −1 ) (resp.α(r −1 h )) is concave on [0, α(a)] (resp.[0, r h (b)]).Then there is at most one intersection point satisfying the conditions (1.5)-(1.6)for nh .Thus, the strict monotonicity of α implies that I is a singleton, hence nh is monomorphic.In the case where r h (α −1 (z)) (resp.α(r −1 h )) is convex, I contains at most two points, thus nh is at most dimorphic.
Note that by Remark 1.3, we can also check if an ESD is monomorphic or dimorphic by studying the set of maximum points of the corresponding fitness functions.We state the above results in the following proposition.
Figure 2: Elements for the analysis for nh .Left: An example when the blue curve is convex and the dimorphic situation occurs.Right: An example when the blue curve is concave and the monomorphic situation occurs Proposition 4.1 (Conditions for monomorphism or dimorphism).Assume that (n h , ns ) is an ESD arbitrarily that does not vanish.Then nh is monomorphic if one of the following hypotheses is fulfilled: (ii) or r h is strictly monotone and α(r Also nh is at most dimorphic if one of the following hypotheses is fulfilled: (i) either α is strictly monotone and r h (α −1 ) is convex on [0, α(a)], (ii) or r h is strictly monotone and α(r −1 h )) is convex on [0, r h (b)], (iii) or r h , α are strictly convex.
Furthermore, the same conclusions as above hold for ns provided that similar assumptions on r s , β are supposed.
)). Suppose furthermore that the pair (n h , ns ) below is nonnegative and that the assumptions (depending on situations) mentioned below hold.Then there exists a unique (non-negative) ESD and it is monomorphic or vanishes.
(i) Under the assumptions (H 1 ) and (H 3 ), the unique ESD is given by (ii) Under the assumptions (H 2 ) and (H 3 ), the unique ESD is given by (iii) Under the assumptions (H 1 ) and (H 4 ), the unique ESD is given by (iv) Under the assumptions (H 2 ) and (H 4 ), the unique ESD is given by Proof.We only prove (i).The other cases can be proved in the same way.First note that under the assumptions (H 1 ) and (H 3 ), β = β * ≥ 1.We have G s (y) = r s (y) − ρh − ρs + β * Σh is strictly decreasing (by (H 3 )) so that it attains its global maximum only at y = 3.This, in view of Remark 1.3, yields that supp ns = {3}.As a consequence, G s (3) = 0 so that This together with the property that α ≤ 0 implies that G h (x) = r h (x) + α (x) Σs ≤ r h (x) + 3α (x)r s (3) < 0 for all x ∈ (1, 2), where we used the hypothesis (H 1 ) in the last inequality.Consequently, G h is strictly decreasing on [1,2] so that it has only one maximum point x = 1.Hence supp nh = {1}.The expression of (n h , ns ) follows from (3.3) with (x, ŷ) := (1, 3).Below, we compute all ESDs explicitly.We also see that the dimorphic situation may occurs.For simplicity, we only consider the dimorphic distribution for hematopoietic stem cells., (c, d) := (3,4).Assume that r h is strictly convex and that α is convex.Suppose furthermore that (H 3 ) is satisfied (i.e., β = β * is a positive constant and r s is strictly decreasing).Then all ESDs are given by provided that ρh ≥ 0, ρs ≥ 0 and provided that ρh ≥ 0, ρs ≥ 0 and , ρs = r h (2) 3α( 1) , provided that ρh1 ≥ 0, ρh2 ≥ 0, ρs ≥ 0.
Remark 4.4.We can also prove similar results as in Proposition 4.3 when we suppose the hypothesis (H 4 ) instead of (H 3 ).
Proof.First note that G h (x) is strictly convex.Thus G h attains its global maximum only at endpoints of the interval [1,2].Note also that G s (y) = r s (y) − ρh − ρs + β * Σh is strictly decreasing so that y = 3 is its unique maximum point.The above observations and Remark 1.3 imply that supp ns = 3 and either supp nh = {1} or supp nh = {2} or supp nh = {1, 2}.
The conditions in the definition 1.2 are equivalent to Solving this system yields the expressions of ρh1 , ρh2 , ρs .

Let us consider two concrete examples below.
Example 1: Suppose that r h , α, r s , β are given by According to Proposition 4.3, there are only two positive ESDs which are • Dimorphic distribution: Example 2: Suppose that r h , α, r s , β are given by Then, there is only one positive ESD.This ESD is dimorphic and has the form

Partial uniqueness of ESDs under homogeneity assumptions on stromal cells
We suppose that all stromal cells have some similar properties.More precisely, we consider the case where the contribution of each phenotype y in the message from stromal cells to HSCs and the sensitivity of stromal cells to the message from HSCs are the same.Mathematically, assume that the weight function ψ s (y) = 1, and that Below we state a result about the partial uniqueness of ESDs.(ii) In addition to (i), if ns does not vanish, then Proof.The definition of ESD and the non-negativity of nh yield that Similarly, it follows from the inequality Summing up (4.9) and (4.10), we obtain Similarly, we deduce from the inequality This together with (4.11) yields Equivalently, Therefore, the hypothesis (4.3) yields that the equality in (4.13) (also in all above inequalities) holds.Thus each term in (4.13) vanishes so that (4.4) follows.The identities (4.5), (4.6) follow from the fact that In the case αΣ h = βΣ s = 0, an entropy functional has been found by Jabin and Raoul [14] and used to prove the convergence of the solution to the unique ESD.In the general form of the system (1.1), we do not expect to find an entropy functional due to the complexity of the terms αΣ h , βΣ s .However, we obtain an entropy functional similar to that of Jabin and Raoul for the system (4.2) corresponding to particular choices of αΣ h and βΣ s .We obtain below a partial information about the dynamics of the solution of (4.2) as t → ∞: the entropy functional decreases monotonically on orbits.The question of the convergence of the solution to the unique ESD remains open but this functional could be an essential ingredient to solve this issue.Proposition 4.6 (Entropy functional).Let (n h , ns ) be an (non-negative) ESD of problem (4.2).Set Then E is an entropy functional for (4.2), i.e., d dt E(t) ≤ 0 provide that β is a positive constant satisfying We have In the above inequality, we have used the definition of ESDs and the non-negativity of n h (cf.Proposition 2.2).
Similarly, for E 2 := n s − ns ln(n s ), we have Therefore, as E = βE 1 + E 2 , we have where the last inequality follows from the negativity of the discriminant of this polynomial: 5 Dynamics of the fittest traits: an asymptotic point of view We are interested in the dynamics of HSCs and stromal cells with initial data close to a monomorphic state and, in particular, in tracking the movements of concentration point towards an ESD.We follows the analysis in [16] and perform the time change variable τ = tε to accelerate time and observe the dynamics.The parameter ε is also used to measure how close is the distribution from the Dirac distribution.The change of variable t → τ converts the system (3.1) to This system is completed with the initial data Rather than working on n ε h , n ε s directly, we define as usual [10,16,5] the functions u ε , v ε given by The functions u ε , v ε satisfy ( Our purpose is to study the behaviour of u ε , v ε as ε → 0 (at least with subsequences).In order to guarantee the existence of a global solution, suppose that (2.3) is fulfilled.Thus, under the assumption that ρ ε h0 + ρ ε s0 is uniformly bounded, Proposition 2.2 yields that there exists (n h , ns ) such that as ε → 0 (after extracting subsequences), (ii) As ε → 0 (after extractions of subsequences), the functions u ε and v ε converge locally uniformly to Lipschitz continuous functions u and v.Moreover, u, v satisfy xn h , max τ,x u(τ, x) ≤ 0, max τ,y v(τ, y) ≤ 0 for all τ ≥ 0. (5.5) Furthermore we have for a.e.τ , Proof.(i) First note that by Proposition 2.2, there is a constant (5.6) Differentiating the equation for u ε with respect to x yields ∂ τ u ε x (τ, x) = r h (x) + α (x) yn ε s (τ, y)dy.
(5.7) Thus, using (5.6), we obtain On the other hand, in view of (5.1), we have Similarly, the same property holds for v ε .(ii) Using the point (i) and the Arzelà-Ascoli Theorem, we may extract subsequences (u ε , v ε ) which converge as indicated in the statement.The equations for u and v are obtained by passing to the limit in (5.1).Moreover, u, v cannot take positive values.Otherwise (ρ ε h , ρ ε s ) blows up in the limit as ε vanishes and this is in contradiction with (5.6).

Monomorphic states
We provide sufficient conditions so that (n h , ns ) defined in (5.2), (5.3) is a monomorphic state.
The parameters are chosen as in the table below.Remark 6.1.The parameters in the three cases satisfy the assumptions of Theorem 4.2 (iii), (ii) and Proposition 4.3 (iii), respectively.Also they are chosen small enough such that the condition (2.3) for the global existence of the solution holds.In the first case (Figures 3 and 4), we use the homogeneous proliferation rate and the inhomogeneous sensitivity of stromal cells.Conversely, in the last two cases, the parameters correspond to the inhomogeous proliferation rate and the homogeneous sensitivity of stromal cells.Note also that in the monomorphic situatations, the corresponding fitness functions are monotone while in the dimorphic situation, the fitness function (for HSCs) is strictly convex.Figure 3 and Figure 4 display the behavior of n h and n s in the time scale t := 10 −2 t.In Figure 3, the population densities for HSCs and their support cells n h , n s are monomorphic and behave as Dirac masses.The concentration point of n h moves towards x = 1 and the one of n s moves towards the point y = 4.In this situation, the phenotype (x, y) = (1, 4) is selected.This represents a good scenario: healthy HSCs and stromal cells with the best support capacity are selected.The evolution of the corresponding dominant phenotypes are given in Figure 4.

Conclusion and perspectives
In this paper, we have introduced a mathematical model for the interaction between hematopoietic stem cells and their support cells.Leukemic stem cells are also taken into account in the model and the phenotype x, characterising the population heterogeneity in a way relevant to the question at stake, represents for both the intrinsic proliferation rate of HSCs and the malignancy potential of cancer cells (i.e., as mentioned in the introduction, a proposed pathological combination of both plasticity and fecundity, likely related to how many mutations are involved in cancer cells).Note also that the monotonicity assumption on r h means that we assumed that the malignant cells proliferate more than healthy HSCs.
We performed a study concerning the adaptive dynamics of HSCs and support cells, in particular, investigating Dirac masses (or sums of Dirac masses) that arose in the solutions of particular cases of the system.Linear stability results for single Dirac mass steady states, suggesting that another phenotype will invade the stationary environment corresponding to the steady state if the corresponding fitness function computed at that phenotype is positive.We also provided sufficient conditions to ensure that ESDs are dimorphic or monomorphic.These conditions are related to the convexity, concavity, monotonicity assumptions of the function parameters.In many cases, we could show the existence and uniqueness of ESDs as well as compute explicitly all ESDs in the case of non-uniqueness.
Applying an asymptotic approach, we showed that without extinction, the population density of HSCs and of their support cells behave as Dirac masses: Here the points of concentration x(τ ), ŷ(τ ) represent well adapted phenotypes at the time τ .Also these points are maximum points of the phase functions u(•, τ ) and v(•, τ ).The system (5.11)gives us the dynamics of x(τ ), ŷ(τ ), in other words, the adaptive process for HSCs and its support cells during their evolution.Our numerical illustrations provide the case of the existence of HSCs, or LCSs (only).Also, we illustrate the case of invasion of LCSs as well as the coexistence of HSCs and LSCs.This latter situation does not seem to be usually seen in the clinic of acute leukemias, which may be due to the fact that in reality, the competition for space and nutrients turns to the advantage of leukemic cells (The biological fact is that stromal cells change to adapt to healthy cells or malignant cells.Thus LSCs and HSCs have different hematopoietic niches so that the competitive strength of HSCs and LSCs for space and nutrients will be different).In our model, the advantage of leukemic cells in competition could be represented by a diversified non local logistic term in the equation for HSCs such as −k 1 (b − x)n h (t, x)dx − k 2 (x − a)n h (t, x)dx, with k 2 > k 1 , instead of the neutral term −ρ h , thus attributing more importance in the competition to cells close to the malignant phenotype x = b.Or else, could it be that actual biological coexistence between HSCs and LSCs could come from the fact that leukemic cells may have been reduced to a state of dormancy?Note that this perspective of dormancy has recently been investigated in a rather different modelling setting (no adaptive dynamics, no interaction with stromal cells) in [11].
Our analytic results, except in Section 4.2, hold for more general choices of the weight functions ψ h , ψ s .However, we present our results here mainly for the case ψ h = x, ψ s = y to clarify the ideas and avoid complex computations.The mathematical question related to the convergence of the solution of (1.1) to its limit (which is an ESD) remains open.The BV-method (see, for example, [23,17]) seems not amenable be applied due to the complexity of function parameters.However, we could find an entropy functional for a simplified system (4.2).This functional decreases to −∞, however it could be an essential ingredient to solve this issue.
The present model and its mathematical analysis represent to the best of our knowledge a first attempt to study the interactions between HSCs and their supporting stromal cells in the framework of adaptive dynamics.One could certainly complicate it to introduce multidimensional phenotypes related to refined cell functionalities such as fecundity, viability, plasticity, in the two cell populations, but even simple as it is, it relies on many unknown functions that should first be experimentally evaluated to go further in this modelling work, which we actually plan to do in the future.

Figure 1 :
Figure 1: An illustration for interacting messages between (healthy and malignant) HSCs and stromal cells.Interacting messages (Σ h , Σ s ) are represented by black arrow lines.The sensitivities of cells are described by red curves: HSCs are very sensitive to interacting messages (α > 0) while cancer cells (at their later stages) are independent of the surrounding stroma (α = 0).Blue curves refer to the intrinsic proliferation rates: HSCs cannot survive without supporting messages (r h = 0) while cancer cells can proliferate without supporting messages (r h > 0).

. 2 )Lemma 2 . 1 .
In view of [23, Theorem 2.1, Page 29], we have the following selection principle: Suppose that r h is bounded and strictly increasing.Suppose furthermore that r s is bounded and strictly decreasing.Assume that n h0 , n s0 ∈ L 1 are positive on [a, b] and [c, d], respectively.

Figure 3 :
Figure 3: Behaviour of HSCs (first row) and stromal cells (second row) in time.Stromal cells with best support capacity are selected and healthy HSCs persist (no LSCs).

Figure 4 :
Figure 4: Evolution of the dominant trait (horizontal axis for the distribution of traits) with time (vertical axis).Left: phenotype x of HSCs.Right: Phenotype y of stromal cells.Monomorphic states for HSCs and stromal cells.Figure3and Figure4display the behavior of n h and n s in the time scale t := 10 −2 t.In Figure3, the population densities for HSCs and their support cells n h , n s are monomorphic and behave as Dirac masses.The concentration point of n h moves towards x = 1 and the one of n s moves towards the point y = 4.In this situation, the phenotype (x, y) = (1, 4) is selected.This represents a good scenario: healthy HSCs and stromal cells with the best support capacity are selected.The evolution of the corresponding dominant phenotypes are given in Figure4.Figures5 and 6below show another monomorphic situation where stromal cells with lowest support capacity are selected.Healthy HSCs cannot survive and cancer cells are selected.
Figure 4: Evolution of the dominant trait (horizontal axis for the distribution of traits) with time (vertical axis).Left: phenotype x of HSCs.Right: Phenotype y of stromal cells.Monomorphic states for HSCs and stromal cells.Figure3and Figure4display the behavior of n h and n s in the time scale t := 10 −2 t.In Figure3, the population densities for HSCs and their support cells n h , n s are monomorphic and behave as Dirac masses.The concentration point of n h moves towards x = 1 and the one of n s moves towards the point y = 4.In this situation, the phenotype (x, y) = (1, 4) is selected.This represents a good scenario: healthy HSCs and stromal cells with the best support capacity are selected.The evolution of the corresponding dominant phenotypes are given in Figure4.Figures5 and 6below show another monomorphic situation where stromal cells with lowest support capacity are selected.Healthy HSCs cannot survive and cancer cells are selected.

Figure 5 :
Figure 5: Behaviour of HSCs (first row) and stromal cells (second row) in time.Stromal cells with lowest support capacity are selected.Healthy HSCs go extinct and LSCs persist.

Figure 6 :
Figure 6: Left: phenotype x of HSCs.Right: Phenotype y of stromal cells (phenotype space in abscissae, time in ordinates).The support is not sufficient and healthy HSCs cannot persist; only LSCs survive.One can notice an apparent fracture between the two populations around the middle of the phenotype space.

Figure 7 :
Figure 7: Behaviour of HSCs (first row) and stromal cells (second row) in time.Stromal cells with lowest support capacity are selected.HSCs and LSCs coexist (carefully note the concentration around both x = 1 and x = 2 in the first row).

Figure 8 :
Figure 8: Left: The dominant phenotypes for HSCs move towards the left and right.The phenotypes x = 1 and x = 2 are selected, i.e., both healthy and malignant HSCs persist.Right: Evolution of dominant phenotypes for stromal cells.Stromal cells with the lowest capacity of the support are selected.

Figures 7
Figures 7 and 8 represent the dimorphic situation for HSCs in the time scale t := 10 −2 t.This situation corresponds to Proposition 4.3 (iii).Starting from an initial distribution with one peak at x = 1.45, a branching process appears.There are two dominant phenotypes of HSCs.The first one moves to the left (only healthy cells selected about the time until t = 10) and the second one move towards x = 2 (and selected in a little bit later).In other words, cancerous cells invade the population of HSCs, however without occupying it in totality: healthy HSCs and cancerous cells (LSCs) coexist.

Table 1 :
Settings in the numerical simulations.