Loading [MathJax]/jax/output/SVG/jax.js
Commentary

Prophylactic aspirin and public

  • The clinical use of aspirin, first synthesized over 100 years ago, entered a new phase in 1974 with the reporting of the first randomised trial showing a reduction in vascular disease deaths from low-doses. More recent evidence suggests that the medicine is effective against cancer, which makes aspirin of very considerable potential importance to public health improvement and potentially also to clinical practice. It appears that prophylactic aspirin is being increasingly used throughout the community. There is need therefore for the risks and benefits of low-dose aspirin, and its' role within healthcare and within public health, to be widely discussed not least as media reports are bringing this issue into the public domain. It also follows that policy decisions need to be taken as to whether or not its use should be actively promoted. In particular, it is important that Doctors and healthcare practitioners are well informed of the risks and benefits so that they can impart this knowledge during consultations. Furthermore, it is important that low-dose aspirin is not perceived as a substitute for a healthy lifestyle, but that it is recommended and uptake monitored alongside other protective behaviours to improve on health gain, such as smoking cessation, moderate alcohol intake, exercise and diet.

    Citation: Gareth Morgan. Prophylactic aspirin and public[J]. AIMS Public Health, 2014, 1(1): 1-8. doi: 10.3934/publichealth.2013.1.1

    Related Papers:

    [1] Qingyun Zhang, Xinhua Zhao, Liang Liu, Tengda Dai . Dynamics analysis of spatial parallel robot with rigid and flexible links. Mathematical Biosciences and Engineering, 2020, 17(6): 7101-7129. doi: 10.3934/mbe.2020365
    [2] Zhen Liu, Song Yang, Tao Ding, Ruimin Chai . Dynamic modeling and performance analysis of a lower-mobility parallel robot with a rotatable platform. Mathematical Biosciences and Engineering, 2023, 20(2): 3918-3943. doi: 10.3934/mbe.2023183
    [3] Jia Tan, ShiLong Chen, ZhengQiang Li . Robust tracking control of a flexible manipulator with limited control input based on backstepping and the Nussbaum function. Mathematical Biosciences and Engineering, 2023, 20(12): 20486-20509. doi: 10.3934/mbe.2023906
    [4] Yunxia Wei, Yuanfei Zhang, Bin Hang . Construction and management of smart campus: Anti-disturbance control of flexible manipulator based on PDE modeling. Mathematical Biosciences and Engineering, 2023, 20(8): 14327-14352. doi: 10.3934/mbe.2023641
    [5] Tongyu Wang, Yadong Chen . Event-triggered control of flexible manipulator constraint system modeled by PDE. Mathematical Biosciences and Engineering, 2023, 20(6): 10043-10062. doi: 10.3934/mbe.2023441
    [6] Dawei Ren, Xiaodong Zhang, Shaojuan Lei, Zehua Bi . Research on flexibility of production system based on hybrid modeling and simulation. Mathematical Biosciences and Engineering, 2021, 18(1): 933-949. doi: 10.3934/mbe.2021049
    [7] Xiaoyu Du, Yijun Zhou, Lingzhen Li, Cecilia Persson, Stephen J. Ferguson . The porous cantilever beam as a model for spinal implants: Experimental, analytical and finite element analysis of dynamic properties. Mathematical Biosciences and Engineering, 2023, 20(4): 6273-6293. doi: 10.3934/mbe.2023270
    [8] Dan Yang, Zhiqiang Xie, Chunting Zhang . Multi-flexible integrated scheduling algorithm for multi-flexible integrated scheduling problem with setup times. Mathematical Biosciences and Engineering, 2023, 20(6): 9781-9817. doi: 10.3934/mbe.2023429
    [9] Xinyu Shao, Zhen Liu, Baoping Jiang . Sliding-mode controller synthesis of robotic manipulator based on a new modified reaching law. Mathematical Biosciences and Engineering, 2022, 19(6): 6362-6378. doi: 10.3934/mbe.2022298
    [10] Zongchao Liu, Linhui Wu, Junwei Yang, Fangsen Cui, Pei Ho, Liping Wang, Jianghui Dong, Gongfa Chen . Thoracic aorta stent grafts design in terms of biomechanical investigations into flexibility. Mathematical Biosciences and Engineering, 2021, 18(1): 800-816. doi: 10.3934/mbe.2021042
  • The clinical use of aspirin, first synthesized over 100 years ago, entered a new phase in 1974 with the reporting of the first randomised trial showing a reduction in vascular disease deaths from low-doses. More recent evidence suggests that the medicine is effective against cancer, which makes aspirin of very considerable potential importance to public health improvement and potentially also to clinical practice. It appears that prophylactic aspirin is being increasingly used throughout the community. There is need therefore for the risks and benefits of low-dose aspirin, and its' role within healthcare and within public health, to be widely discussed not least as media reports are bringing this issue into the public domain. It also follows that policy decisions need to be taken as to whether or not its use should be actively promoted. In particular, it is important that Doctors and healthcare practitioners are well informed of the risks and benefits so that they can impart this knowledge during consultations. Furthermore, it is important that low-dose aspirin is not perceived as a substitute for a healthy lifestyle, but that it is recommended and uptake monitored alongside other protective behaviours to improve on health gain, such as smoking cessation, moderate alcohol intake, exercise and diet.


    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 [1], inspired by Burns and Tannock [2] and Lajtha [3], 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, [5,6,7,8]—models in the form of ODEs or age-structured transport equations with applications to chronic myelogenous leukaemia, [9]—a diffusion model including spatial competition between cells–, reviews [4,10,11] 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, [12,13,14] for reviews of the interaction between HSCs and stromal cells, [15] 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.

    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 \in (a, b), y \in (c, d) $, where $ 0 < a < b $ and $ 0 < c < d $. Totally healthy HSCs will thus have a phenotype $ x $ close to $ a $, while agammaessive leukemic HSCs (i.e., LSCs) will have a phenotype $ x $ close to $ b $. We consider a mathematical model of the form

    $ {tnh(t,x)=[rh(x)ρh(t)ρs(t)+α(x)Σs(t)]nh,x(a,b),t>0,tns(t,y)=[rs(y)ρh(t)ρs(t)+β(y)Σh(t)]ns,y(c,d),t>0. $ (1.1)

    This system is completed with initial data

    $ nh(0,x)=nh0(x)0,ns(0,y)=ns0(y)0. $ (1.2)

    Here our assumptions and notations are

    ● $ \rho_h(t): = \int_a^b n_h(t, x)\, dx, \rho_s(t): = \int_c^d n_s(t, y)\, dy $ are the total populations of HSCs and their support cells, respectively.

    ● The terms $ \Sigma_h(t): = \int_a^b \psi_h (x) n_h(t, x)\, dx, \Sigma_s(t): = \int_c^d \psi_s (y) n_s(t, y)\, dy $ denote an assumed chemical signal ($ \Sigma_h $) from the hematopoietic immature stem cells (HSCs) to their supporting stroma (MSCs), i.e., ``call for support'' and conversely, a trophic message ($ \Sigma_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 [13,14]. The nonnegative functions $ \psi_h, \psi_s $ defined on $ (a, b) $ and $ (c, d) $ measure the contribution of each phenotype in the interactive messages. Assume that $ \psi'_s \ge 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 $ \psi'_h \ge 0 $.

    ● The term $ r_h \ge 0 $ represents the intrinsic (i.e., without contribution from trophic messages from MSCs, nor limitation by the non local logistic terms $ \rho_h $ and $ \rho_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 $ \alpha \ge 0 $, satisfying $ \alpha' \le 0 $ and $ \alpha(b) = 0 $, is the sensitivity of HSCs to the trophic messages from support cells

    ● For the term $ r_s \ge 0 $, we assume that $ r'_s(y) \le 0 $ (there is a cost in proliferation for support cells to increase their support capacity). The term $ \beta(y) \ge 0 $ with $ \beta'(y) \ge 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, prey-predator and cooperative types are typical examples of such interaction; see, for example, [27,Chapter 3]. Apart from the cases mentioned in [27,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. [16]). We refer to [17] a related competitive system with healthy and cancer cells structured by a phenotypic variable related with their resistance to chemotherapy, to [18] an integro-differential Lotka-Volterra system for the interaction of $ N $ populations $ (N \ge 2) $. In our model, besides the competition terms between cells, we introduce new terms $ \Sigma_h, \Sigma_s, \alpha, \beta $ 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 $ -\rho_s(t)+\alpha(x)\Sigma_s(t) $ and $ -\rho_h(t)+\beta(y)\Sigma_h(t) $. Note that if $ \Sigma_h = \Sigma_s = 0 $, our model reduces to the cases studied in [16,17]. Also when $ \psi_h = \psi_s = 1 $, our problem becomes a particular case of [18].

    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 agammaessive is a HSC (in fact, a LSC, since it is then cancerous). Also the monotonicity of $ \alpha $ (this function stands for the sensitivity of HSCs to the trophic messages coming from MSCs) indicates that the more agammaessive $ 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 $ \alpha(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, \beta $ 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, \alpha, r_h, r_h, \psi_h, \psi_s $ can be chosen as linear or quadratic functions. For example, $ r_h $, $ \alpha $ are given by $ r_h = r_h^*(x-a) $ or $ r_h = r_h^*(x-a)^2 $, $ \alpha(x) = \alpha^*(b-x) $ with positive constants $ r_h^*, \alpha^* $, $ \psi_h(x) = x $, $ \psi_s(y) = y $.

    A more general model which includes the possibility of mutations has the form:

    $ {tnh(t,x)=μh(nh)xx+[rh(x)c11(x)ρh(t)c12(x)ρs(t)+α(x)Σs(t)]nh,tns(t,y)=μs(ns)yy+[rs(y)c21(y)ρh(t)c22(y)ρs(t)+β(y)Σh(t)]ns. $ (1.3)

    Here the diffusion terms represent mutation with rates $ \mu_h, \mu_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 $ \mu_h = \mu_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 $ \mu_h, \mu_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).

    Hereafter, we will use the notation $ \int $ to denote the integrals over $ [a, b] $ and $ [c, d] $, as long as there is no risk of confusion. Let $ \hat n_h, \hat n_s $ be measures defined on $ [a, b], [c, d] $, respectively. We set

    $ {\rm{supp}}\ \hat n_h = I \subset [a, b], \quad {\rm{supp}}\ \hat n_s = J \subset [c, d], $
    $ \hat \rho_h = \int \hat n_h, \quad \hat \rho_s = \int \hat n_s, \quad \hat \Sigma_h = \int \psi_h(x) \hat n_h, \quad \hat \Sigma_s = \int \psi_s(y) \hat n_s, $
    $ Gh(x):=rh(x)ˆρhˆρs+α(x)ˆΣs,Gs(y):=rs(y)ˆρhˆρs+β(y)ˆΣh. $ (1.4)

    Definition 1.1. The pair $ (\hat n_h, \hat n_s) $ is a steady state of (1.1) if

    $ Gh(x)=0,Gs(y)=0  for all  xI,yJ. $ (1.5)

    Furthermore, $ \hat n_h $ (resp. $ \hat n_s $) is said

    (ⅰ) to be monomorphic if $ {\rm{supp}}\ \hat n_h $ (resp. $ {\rm{supp}}\ \hat n_s $) is a singleton,

    (ⅱ) to be dimorphic if $ {\rm{supp}}\ \hat n_h $ (resp. $ {\rm{supp}}\ \hat n_s $) is a set of two points.

    Definition 1.2. We say that $ (\hat n_h, \hat n_s) $ is an evolutionary stable distribution (ESD) of problem (1.1) if it is a steady state and the condition below is fulfilled

    $ Gh(x)0 for all x[a,b]I,Gs(y)0 for all y[c,d]J. $ (1.6)

    Remark 1.3. It follows from Definition 1.2 that any $ x \in I $ (resp. $ y \in J $) is a maximum point of $ G_h $ (resp. $ G_s $).

    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 \in I, y \in J $ are living in the stationary environment given by $ (\hat \rho_h, \hat \rho_s, \hat \Sigma_h, \hat \Sigma_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 [20,25] for more details about the framework of adaptive dynamics.

    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, $ \tau = \varepsilon t $, with small $ \varepsilon > 0 $. We represent the dynamics of dominant traits (Theorems 5.1 and 5.2) in the regime as $ \varepsilon \to 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.

    In the absence of stromal cells, the system (1.1) reduces to the equation

    $ tnh(t,x)=[rh(x)ρh(t)]nh,x(a,b),t>0. $ (2.1)

    Similarly, the behaviour of stromal cells without HSCs is given by

    $ tns(t,y)=[rs(y)ρs(t)]ns,y(c,d),t>0. $ (2.2)

    In view of [28,Theorem 2.1,Page 29], we have the following selection principle:

    Lemma 2.1. 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} \in L^1 $ are positive on $ [a, b] $ and $ [c, d] $, respectively.

    (ⅰ) For (2.1), we have $ n_h(t, x) \rightarrow r_h(b) \delta_{\{x = b\}} $ weakly in the sense of measures as $ t \to \infty $.

    (ⅱ) For (2.2), we have $ n_h(t, y) \rightarrow r_s(c) \delta_{\{y = c\}} $ weakly in the sense of measures as $ t \to \infty $.

    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(ⅰ) 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(ⅱ)).

    For a function $ f $ defined on an interval $ I $, we set

    $ \underline f: = \min\limits_{x \in I} f (x), \quad \overline f: = \max\limits_{x \in I} f (x). $

    In the results below, we only need the boundedness of $ r_h, r_s, \alpha, \beta, \psi_h, \psi_s $. We do not need the monotonicity assumption of these functions.

    Proposition 2.2. Assume that $ r_h, r_s, \alpha, \beta, \psi_h, \psi_s $ are non-negative and bounded. Suppose furthermore that

    $ ˉαˉψs+ˉβˉψh<4. $ (2.3)

    Set

    $ \rho^M: = \max \left( \rho_{h0}+ \rho_{s0}, \dfrac{4 \max(\bar r_h, \bar r_s)}{4-(\bar \alpha \bar \psi_s+\bar \beta \bar \psi_h)}\right). $

    Then the solution of (1.1)- (1.2) is non-negative and satisfies

    $ 0ρh(t)+ρs(t)ρMfor all  t0. $ (2.4)

    Proof. First note that the solution of (1.1) can be written in the form

    $ n_h(t, x) = n_{h0}(x) \exp \left( \int_0^t A(\sigma, x) \, d \sigma \right), \quad n_s(t, y) = n_{s0}(y) \exp \left( \int_0^t B(\sigma, y) \, d \sigma\right), $

    where

    $ A(\sigma, x) = r_h(x) -\rho_h(\sigma)-\rho_s(\sigma) +\alpha(x)\Sigma_s(\sigma), \quad B(\sigma, y) = r_s(y) -\rho_h(\sigma)-\rho_s(\sigma) +\beta(y)\Sigma_h(\sigma). $

    Thus $ n_h(t, x) \ge 0, n_s(t, y) \ge 0 $ since $ n_{h0}(x) \ge 0, n_{s0}(y) \ge 0 $. This yields the first inequality of (2.4).

    Integrating the two equations in (1.1) yields

    $ ddtρh=ρ2hρhρs+rh(x)nh+Σs(t)α(x)nh,ddtρs=ρ2sρhρs+rs(y)ns+Σh(t)β(y)ns. $

    Summing up these two identities and using the non-negativity of $ n_h, n_s $, we obtain

    $ ddt(ρh+ρs)=(ρh+ρs)2+rh(x)nh+Σs(t)α(x)nh+rs(y)ns+Σh(t)β(y)ns(ρh+ρs)2+max(ˉrh,ˉrs)(ρh+ρs)+ˉαˉψsρhρs+ˉβˉψhρhρs(ρh+ρs)2+max(ˉrh,ˉrs)(ρh+ρs)+ˉαˉψs+ˉβˉψh4(ρh+ρs)2=[max(ˉrh,ˉrs)(1ˉαˉψs+ˉβˉψh4)(ρh+ρs)](ρh+ρs). $

    Hence the second inequality of (2.4) follows.

    Proposition (2.2) and a standard argument imply the following result.

    Theorem 2.3. Let $ (n_{h0}, n_{s0}) \in L^1(a, b) \times L^1(c, d) $ be non-negative. Then (1.1) possesses a unique solution $ (n_h, n_s) \in C^1([0, \infty); L^1(a, b) \times L^1(c, d)) $.

    Hereafter, for simplicity, we assume that $ \psi_h(x) = x, \psi_s(y) = y $. Then Problem (1.1) becomes

    $ {tnh(t,x)=[rh(x)ρh(t)ρs(t)+α(x)yns(t,y)dy]nh,x(a,b),t>0,tns(t,y)=[rs(y)ρh(t)ρs(t)+β(y)xnh(t,x)dx]ns,y(c,d),t>0. $ (3.1)

    Let $ \hat x \in [a, b], \hat y \in [c, d] $. Consider a particular case where hematopoietic and support cells evolve as single Dirac masses concentrated at $ \hat x, \hat y $. In other words, we focus on the behaviour of the size of the populations. In this case, $ n_h, n_s $ have the form

    $ n_h(t, x) = \rho_h(t) \delta_{\{x = \hat x\}}, \quad n_s(t, y) = \rho_s(t) \delta_{\{y = \hat y\}}, $

    where $ \rho_h(t), \rho_s(t) $ satisfy

    $ {ddtρh=(rh(ˆx)ρh(1α(ˆx)ˆy)ρs)ρh,ddtρs=(rs(ˆy)(1β(ˆy)ˆx)ρhρs)ρs. $ (3.2)

    We can obtain an explicit form of single Dirac mass steady states.

    Lemma 3.1. Let $ \hat x \in [a, b], \hat y \in [c, d] $ and assume that $ 1-(1-\alpha(\hat x) \hat y)(1-\beta(\hat y)\hat x)\! \neq\! 0 $. Set

    $ ˆρh:=rh(ˆx)rs(ˆy)(1α(ˆx)ˆy)1(1α(ˆx)ˆy)(1β(ˆy)ˆx),ˆρs:=rs(ˆy)rh(ˆx)(1β(ˆy)ˆx)1(1α(ˆx)ˆy)(1β(ˆy)ˆx). $ (3.3)

    Then $ (\hat n_h = \hat \rho_h \delta_{\{x = \hat x\}}, \hat n_s = \hat \rho_s \delta_{\{y = \hat y\}}) $ is a steady state of problem (3.1). If the three conditions below hold

    $  {1(1α(ˆx)ˆy)(1β(ˆy)ˆx)>0,rh(ˆx)rs(ˆy)(1α(ˆx)ˆy)>0,rs(ˆy)rh(ˆx)(1β(ˆy)ˆx)>0, $ (3.4)

    then $ \hat \rho_h > 0, \hat \rho_s > 0 $ and $ (\hat \rho_h, \hat \rho_s) $ is a linearly stable steady state of (3.2).

    Proof. Single Dirac mass steady state solution yields

    $ {rh(ˆx)ˆρhˆρs+α(ˆx)ˆyˆρs=0,rs(ˆy)ˆρhˆρs+β(ˆy)ˆxˆρh=0. $ (3.5)

    An elementary calculation shows that $ (\hat \rho_h, \hat \rho_s) $ is a steady state of (3.2) and the corresponding Jacobian matrix is given by

    $ J=(ˆρh(1α(ˆx)ˆy)ˆρh(1β(ˆy)ˆx)ˆρsˆρs). $ (3.6)

    In view of (3.4), $ {\rm{tr}}(\mathcal J) < 0, \det(\mathcal J) > 0 $ so that the two eigenvalues of $ \mathcal J $ are negative. Thus the linear stability of $ (\hat \rho_h, \hat \rho_s) $ follows.

    The results below are concerned with the linear stability of single Dirac mass steady states among perturbations of particular forms.

    Theorem 3.2 (Stability of monomorphic steady states). Let $ (\hat n_h = \hat \rho_h \delta_{\{x = \hat x\}}, \hat n_s = \hat \rho_s \delta_{\{y = \hat y\}}) $ be a steady state of (3.1) as in Lemma 3.1. Assume that (3.4) holds. Let $ x^*, y^* $ satisfy $ x^* \neq \hat x, y^* \neq \hat y $ and

    $ Gh(x)<0,Gs(y)<0. $ (3.7)

    Then $ (\hat n_h, \hat n_s) $ is linearly stable among perturbations starting by

    $ n_{h0} : = \varepsilon_1 \delta_{\{x = x^*\}}+(\hat \rho_h+\varepsilon_2) \delta_{\{x = \hat x\}}, n_{s0} : = \varepsilon_3 \delta_{\{y = y^*\}} + (\hat \rho_s +\varepsilon_4) \delta_{\{y = \hat y\}}. $

    Proof. We linearize the system at $ (\hat n_h, \hat n_s) $. For $ g_h(t, x): = n_h(t, x)-\hat n_h(x), g_s(t, y): = n_s(t, y)-\hat n_s(y) $ we obtain

    $ {tgh=[rh(x)ˆρhˆρs+α(x)yˆns]ghˆnhgh+ˆnh[gs+α(x)ygs],tgs=[gh+β(y)xgh]ˆns+[rs(y)ˆρhˆρs+β(y)xˆnh]gsˆnsgs. $

    Note that $ g_h, g_s $ have the form

    $ g_h(t, x) = g_h^1(t) \delta_{\{x = x^*\}}+ g_h^2(t)\delta_{\{x = \hat x\}}, \quad g_s(t, y) = g_s^1(t) \delta_{\{y = y^*\}}+ g_s^2(t)\delta_{\{y = \hat y\}}. $

    Therefore, we obtain

    $ {tg1h=[rh(x)ˆρhˆρs+α(x)ˆyˆρs]g1h,tg1s=[rs(y)ˆρhˆρs+β(y)ˆxˆρh]g1s,tg2h=ˆρhg1h+ˆρh[α(ˆx)y1]g1sˆρhg2h+ˆρh(α(ˆx)ˆy1)g2s,tg2s=ˆρs(β(ˆy)x1)g1hˆρsg1s+ˆρs(β(ˆy)ˆx1)g2hˆρsg2s. $

    The corresponding matrix is given by

    $ \left( rh(x)ˆρhˆρs+α(x)ˆyˆρs0000rs(y)ˆρhˆρs+β(y)ˆxˆρh00ˆρhˆρh[α(ˆx)y1]ˆρhˆρh(α(ˆx)ˆy1)ˆρs(β(ˆy)x1)ˆρsˆρs(β(ˆy)ˆx1)ˆρs \right). $

    The eigenvalues of the above matrix are

    $ r_h(x^*) -\hat \rho_h-\hat \rho_s+\alpha(x^*) \hat y \hat \rho_s, \quad r_s(y^*) -\hat \rho_h-\hat \rho_s+ \beta(y^*) \hat x \hat \rho_h $

    and the two eigenvalues of the matrix $ \mathcal J $ in (3.6). All are negative due to (3.7) and Lemma 3.1. Thus the stability of $ (\hat n_h, \hat n_s) $ follows.

    Interpretations of Theorem 3.2: With the notation (1.4), in view of (3.5), we have $ G_h(\hat x) = 0, G_s(\hat y) = 0 $. Also the conditions $ G_h(x^*) < 0, G_s(y^*) < 0 $ show that the phenotypes $ x^*, y^* $ cannot invade the stable equilibrium $ (\hat n_h, \hat n_s) $. As a consequence, no new equilibrium can be reached but a mutant invading the resident population $ (\hat x, \hat y) $.

    Recall that from the defintion 1.2, an ESD $ (\hat n_h, \hat n_s) $ is characterised by the conditions (1.5)– (1.6). Graphically, we plot the curve $ x\in [a, b] \mapsto (Z = \alpha(x), W = r_h(x)) $ by the blue curve and the red straight line $ Z \hat \Sigma_s + W = \hat \rho_h + \hat \rho_s $; see Figure 2. Then the conditions (1.5)– (1.6) for $ \hat n_h $ mean that the blue curve must be below the red line and that the pair $ (\alpha(x), r_h(x)) $ for all $ x \in I: = {\rm{supp}}\ \hat n_h $ are the coordinates of the intersection points between the blue curve and the red line. Similarly, we have the same illustration for $ \hat n_s $.

    Figure 1.  An illustration for interacting messages between (healthy and malignant) HSCs and stromal cells. Interacting messages $ (\Sigma_h, \Sigma_s) $ are represented by black arrow lines. The sensitivities of cells are described by red curves: HSCs are very sensitive to interacting messages ($ \alpha > 0 $) while cancer cells (at their later stages) are independent of the surrounding stroma $ (\alpha = 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 $).
    Figure 2.  Elements for the analysis for $ \hat n_h $. 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.

    If $ \alpha $ (resp. $ r_h $) is strictly monotone and if $ r_h(\alpha^{-1}) $ (resp. $ \alpha(r_h^{-1}) $) is concave on $ [0, \alpha(a)] $ (resp. $ [0, r_h(b)] $). Then there is at most one intersection point satisfying the conditions (1.5)– (1.6) for $ \hat n_h $. Thus, the strict monotonicity of $ \alpha $ implies that $ I $ is a singleton, hence $ \hat n_h $ is monomorphic. In the case where $ r_h(\alpha^{-1}(z)) $ (resp. $ \alpha(r_h^{-1}) $) is convex, $ I $ contains at most two points, thus $ \hat n_h $ 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.

    Proposition 4.1 (Conditions for monomorphism or dimorphism). Assume that $ (\hat n_h, \hat n_s) $ is an ESD arbitrarily that does not vanish. Then $ \hat n_h $ is monomorphic if one of the following hypotheses is fulfilled:

    (ⅰ) either $ \alpha $ is strictly monotone and $ r_h(\alpha^{-1}) $ is concave on $[0,\alpha (a)]$,

    (ⅱ) or $ r_h $ is strictly monotone and $ \alpha(r_h^{-1}) $) is concave on $ [0, r_h(b)] $,

    (ⅲ) or $ r_h, \alpha $ are strictly concave.

    Also $ \hat n_h $ is at most dimorphic if one of the following hypotheses is fulfilled:

    (ⅰ) either $ \alpha $ is strictly monotone and $ r_h(\alpha^{-1}) $ is convex on $ [0, \alpha(a)] $,

    (ⅱ) or $ r_h $ is strictly monotone and $ \alpha(r_h^{-1}) $) is convex on $ [0, r_h(b)] $,

    (ⅲ) or $ r_h, \alpha $ are strictly convex.

    Furthermore, the same conclusions as above hold for $ \hat n_s $ provided that similar assumptions on $ r_s, \beta $ are supposed.

    The next result is concerned with the existence and uniqueness of ESDs. We also show that the unique ESD is monomorphic and that the concentration points are endpoints of the intervals $ (a, b), (c, d) $. For simplicity, set $ (a, b): = (1, 2), (c, d): = (3, 4) $. Here we employ the two distinct sets $ (1, 2) $ and $ (3, 4) $ to insist on that fact that the two phenotypes of HSCs and of stromal cells are not the same. We will suppose two of the following assumptions:

    $ {\bf (H_1)} \ r_h'(x)+ 3 \alpha'(x) \bar r_s < 0 $ for all $ x \in (1, 2) $ and $ \beta(y)\ge 1 $ for all $ y \in (3, 4) $,

    $ {\bf (H_2)} \ r_h'(x)+ 4 \alpha'(x) \bar r_s > 0 $ for all $ x \in (1, 2) $ and $ \beta(y) \le 1/2 $ for all $ y \in (3, 4) $,

    $ {\bf (H_3)} \ \beta = \beta^* > 0 $ is constant and $ r_s $ is strictly decreasing,

    $ {\bf (H_4)} \ \beta $ is strictly increasing and $ r_s = r_s^* > 0 $ is constant.

    Theorem 4.2 (Existence and uniqueness of ESDs). Set $ (a, b): = (1, 2), (c, d): = (3, 4) $. Suppose that $ r_h, \alpha \in C([1, 2]) \cap C^1((1, 2)), r_s, \beta \in C([3, 4]) \cap C^1((3, 4)) $. Suppose furthermore that the pair $ (\hat n_h, \hat n_s) $ below is non-negative and that the assumptions (depending on situations) mentioned below hold. Then there exists a unique (non-negative) ESD and it is monomorphic or vanishes.

    (ⅰ) Under the assumptions $ \bf (H_1) $ and $ \bf (H_3) $, the unique ESD is given by

    $ \hat n_h = \frac{r_s(3)(3 \alpha(1)-1)}{1+(3 \alpha(1)-1)(1-\beta^*)} \delta_{\{x = 1\}}, \quad \hat n_s = \frac{r_s(3)}{1+(3 \alpha(1)-1)(1-\beta^*)}\delta_{\{y = 3\}}. $

    (ⅱ) Under the assumptions $ \bf (H_2) $ and $ \bf (H_3) $, the unique ESD is given by

    $ \hat n_h = \frac{r_h(2)-r_s(3)}{2\beta^*} \delta_{\{x = 2\}}, \quad \hat n_s = \frac{r_s(3)-r_h(2)(1-2\beta^*)}{2\beta^*}\delta_{\{y = 3\}}. $

    (ⅲ) Under the assumptions $ \bf (H_1) $ and $ \bf (H_4) $, the unique ESD is given by

    $ \hat n_h = \frac{r_s^*(4 \alpha(1)-1)}{1+(4 \alpha(1)-1)(1-\beta(4))} \delta_{\{x = 1\}}, \quad \hat n_s = \frac{r_s^*}{1+(4 \alpha(1)-1)(1-\beta(4))}\delta_{\{y = 4\}}. $

    (ⅳ) Under the assumptions $ \bf (H_2) $ and $ \bf (H_4) $, the unique ESD is given by

    $ \hat n_h = \frac{r_h(2)-r_s^*}{2\beta(4)} \delta_{\{x = 2\}}, \quad \hat n_s = \frac{r_s^*-r_h(2)(1-2\beta(4))}{2\beta(4)}\delta_{\{y = 4\}}. $

    Proof. We only prove (ⅰ). The other cases can be proved in the same way. First note that under the assumptions $ \bf (H_1) $ and $ \bf (H_3) $, $ \beta = \beta^* \ge 1 $. We have $ G_s(y) = r_s(y)-\hat \rho_h-\hat \rho_s + \beta^* \hat \Sigma_h $ is strictly decreasing (by $ \bf (H_3) $) so that it attains its global maximum only at $ y = 3 $. This, in view of Remark 1.3, yields that $ {\rm{supp}}\ \hat n_s = \{3\} $. As a consequence, $ G_s(3) = 0 $ so that

    $ r_s(3)-\hat \rho_s = \hat \rho_h -\beta^* \int_1^2 x \hat n_h \le \hat \rho_h -\beta^* \int_1^2 \hat n_h = \hat \rho_h -\beta^*\hat \rho_h \le 0. $

    Therefore,

    $ 3 r_s(3) \le 3 \hat \rho_s \le \int_3^4 y \hat n_s = \hat \Sigma_s. $

    This together with the property that $ \alpha' \le 0 $ implies that

    $ G_h'(x) = r_h'(x)+\alpha'(x) \hat \Sigma_s \le r_h'(x) +3 \alpha'(x)r_s(3) \lt 0 \;\;\text{ for all}\; \; x\in (1, 2), $

    where we used the hypothesis $ \bf (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 $ {\rm{supp}}\ \hat n_h = \{1\} $. The expression of $ (\hat n_h, \hat n_s) $ follows from (3.3) with $ (\hat x, \hat y): = (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.

    Proposition 4.3 (Explicit formulas of all ESDs). Set $ (a, b): = (1, 2), (c, d): = (3, 4) $. Assume that $ r_h $ is strictly convex and that $ \alpha $ is convex. Suppose furthermore that $ (\bf H_3) $ is satisfied (i.e., $ \beta = \beta^* $ is a positive constant and $ r_s $ is strictly decreasing). Then all ESDs are given by

    (ⅰ) $ (\hat n_h = \hat \rho_h \delta_{\{x = 1\}}, \hat n_s = \hat \rho_s \delta_{\{y = 3\}}) $ with

    $ \hat \rho_h = \frac{r_s(3)(3 \alpha(1)-1)}{1+(3 \alpha(1)-1)(1-\beta^*)}, \quad \hat \rho_s = \frac{r_s(3)}{1+(3 \alpha(1)-1)(1-\beta^*)}, $

    provided that $ \hat \rho_h \ge 0, \hat \rho_s \ge 0 $ and

    $ rh(2)3α(1)rs(3)1+(3α(1)1)(1β). $ (4.1)

    (ⅱ) $ (\hat n_h = \hat \rho_h \delta_{\{x = 2\}}, \hat n_s = \hat \rho_s \delta_{\{y = 3\}}) $ with

    $ \hat \rho_h = \frac{r_h(2)-r_s(3)}{2\beta^*}, \quad \hat \rho_s = \frac{r_s(3)-r_h(2)(1-2\beta^*)}{2\beta^*}, $

    provided that $ \hat \rho_h\ge 0, \hat \rho_s \ge 0 $ and

    $ r_h(2) \ge 3 \alpha(1)\frac{r_s(3)-r_h(2)(1-2\beta^*)}{2\beta^*}. $

    (ⅲ) $ (\hat n_h = \hat \rho_{h1} \delta_{\{x = 1\}}+\hat \rho_{h2} \delta_{\{x = 2\}}, \hat n_s = \hat \rho_s \delta_{\{y = 3\}}) $ with

    $ ˆρh1=2rh(2)3α(1)13α(1)rh(2)rs(3)β,ˆρh2=rh(2)rs(3)βrh(2)3α(1)13α(1),ˆρs=rh(2)3α(1), $

    provided that $ \hat \rho_{h1}\ge 0, \hat \rho_{h2}\ge 0, \hat \rho_s\ge 0 $.

    Remark 4.4. We can also prove similar results as in Proposition 4.3 when we suppose the hypothesis $ \bf (H_4) $ instead of $ \bf (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)-\hat \rho_h -\hat \rho_s +\beta^*\hat \Sigma_h $ is strictly decreasing so that $ {y = 3} $ is its unique maximum point. The above observations and Remark 1.3 imply that $ {\rm{supp}}\ \hat n_s = {3} $ and $ \text{either } {\rm{supp}}\ \hat n_h = \{1\} \text{ or } {\rm{supp}}\ \hat n_h = \{2\} \text{ or } {\rm{supp}}\ \hat n_h = \{1, 2\}. $

    (ⅰ) The case $ {\rm{supp}}\ \hat n_h = \{1\}, {\rm{supp}}\ \hat n_s = \{3\} $. The expressions of $ \hat \rho_h, \hat \rho_s $ follows from (3.3) with $ \hat x = 1, \hat y = 3 $. Because of the convexity of $ G_h $ and the decreasing monotonicity of $ G_s $, the condition (1.6) is equivalent to $ G_h(1) \ge G_h(2) $ which implies (4.1). Similarly, Item (ⅱ) — the case $ {\rm{supp}}\ \hat n_h = \{2\}, {\rm{supp}}\ \hat n_s = \{3\} $ — can be treated in the same way.

    (ⅲ) The case $ {\rm{supp}}\ \hat n_h = \{1, 2\}, {\rm{supp}}\ \hat n_s = \{3\} $. The pair $ (\hat n_h, \hat n_s) $ have the form

    $ \hat n_h = \hat \rho_{h1} \delta_{\{x = 1\}}+\hat \rho_{h2} \delta_{\{x = 2\}}, \quad \hat n_s = \hat \rho_s \delta_{\{y = 3\}}. $

    Thus we have

    $ \hat \rho_h = \hat \rho_{h1}+ \hat \rho_{h2}, \quad \hat \Sigma_h = \hat \rho_{h1}+2 \hat \rho_{h2}, \quad \hat \Sigma_s = 3 \hat \rho_{s}. $

    The conditions in the definition 1.2 are equivalent to

    $ G_h(1) = G_h(2) = 0, \quad G_s(3) = 0, $

    that is

    $ {ˆρh1ˆρh2ˆρs+3α(1)ˆρs=0,rh(2)ˆρh1ˆρh2ˆρs=0,rs(3)ρh1ρh2ˆρs+β(ˆρh1+2ˆρh2)=0. $

    Solving this system yields the expressions of $ \hat \rho_{h1}, \hat \rho_{h2}, \hat \rho_{s} $.

    Let us consider two concrete examples below.

    Example 1: Suppose that $ r_h, \alpha, r_s, \beta $ are given by

    $ r_h(x) = (x-1)^2, \quad \alpha(x) = 2-x, \quad r_s(y) = \frac{1}{2}+\frac{1}{4}(3-y), \quad \beta = 0.5. $

    According to Proposition 4.3, there are only two positive ESDs which are

    ● Monomorphic distribution:

    $ \hat n_h(x) = \frac{1}{2} \delta_{\{x = 1\}}, \quad \hat n_s(y) = \frac{1}{2} \delta_{\{y = 3\}}, $

    ● Dimorphic distribution:

    $ \hat n_h(x) = \frac{1}{3} \delta_{\{x = 1\}}+\frac{1}{3} \delta_{\{x = 2\}}, \quad \hat n_s(y) = \frac{1}{3} \delta_{\{y = 3\}}. $

    Example 2: Suppose that $ r_h, \alpha, r_s, \beta $ are given by

    $ r_h(x) = 0.75(x-1)^2, \quad \alpha(x) = 0.625(2-x), \quad r_s(y) = \frac{1}{2}+\frac{1}{4}(3-y), \quad \beta = 0.5. $

    Then, there is only one positive ESD. This ESD is dimorphic and has the form

    $ \hat n_h(x) = \frac{1}{3} \delta_{\{x = 1\}}+\frac{1}{6} \delta_{\{x = 2\}}, \quad \hat n_s(y) = \frac{2}{5} \delta_{\{y = 3\}}. $

    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 $ \psi_s(y) = 1 $, and that $ \beta > 0 $ is constant. Suppose furthermore that $ \alpha(x) = b-x, \psi_h(x) = x $. Problem (1.1) becomes

    $ {tnh(t,x)=[rh(x)ρh(t)ρs(t)+(bx)ns(t,y)dy]nh,x(a,b),t>0,tns(t,y)=[rs(y)ρh(t)ρs(t)+βxnh(t,x)dx]nsy(c,d),t>0. $ (4.2)

    Below we state a result about the partial uniqueness of ESDs.

    Theorem 4.5 (Partial uniqueness results of ESDs). Assume that $ (\hat n_h, \hat n_s) $ and $ (\tilde n_h, \tilde n_s) $ are two (non-negative) ESDs of the system (4.2). Assume further that $ \beta $ is a positive constant satisfying

    $ (β(1b)+1)2<4β. $ (4.3)

    Then

    $ ˆρh=˜ρh=:H,ˆρs=˜ρs=:S, $ (4.4)
    $ rh(x)ˆnhSxˆnh=rh(x)˜nhSx˜nh, $ (4.5)
    $ rs(y)ˆnsβSxˆnh=rs(y)˜nsβSx˜nh. $ (4.6)

    Moreover,

    (ⅰ) If $ r_s $ is strictly decreasing, then $ \hat n_s = \tilde n_s $ either is monomorphic concentrated at $ y = c $ or vanishes. Also we have

    $ rh(x)ˆnh=rh(x)˜nh. $ (4.7)

    (ⅱ) In addition to (i), if $ \hat n_s $ does not vanish, then

    $ xˆnh=x˜nh. $ (4.8)

    Proof. The definition of ESD and the non-negativity of $ \hat n_h $ yield that

    $ 0[rh(x)˜ρh˜ρs+(bx)˜ns]ˆnh={[rh(x)˜ρh˜ρs+(bx)˜ns][rh(x)ˆρhˆρs+(bx)ˆns]}ˆnh=[(ˆρh˜ρh)+(ˆρs˜ρs)]ˆρh+(˜ρsˆρs)(bx)ˆnh. $ (4.9)

    Similarly, it follows from the inequality

    $ \int \left[ r_h(x)-\hat \rho_h-\hat \rho_s +(b-x) \int \hat n_s \right] \tilde n_h \le 0, $

    that

    $ [(˜ρhˆρh)+(˜ρsˆρs)]˜ρh+(ˆρs˜ρs)(bx)˜nh0. $ (4.10)

    Summing up (4.9) and (4.10), we obtain

    $ (ˆρh˜ρh)2+(ˆρs˜ρs)(ˆρh˜ρh)+(˜ρsˆρs)(bx)(ˆnh˜nh)0. $ (4.11)

    Similarly, we deduce from the inequality

    $ [rs(y)˜ρh˜ρs+βx˜nh]ˆns+[rs(y)ˆρhˆρs+βxˆnh]˜ns0, $ (4.12)

    that

    $ (\hat \rho_s-\tilde \rho_s)^2 + (\hat \rho_s-\tilde \rho_s) (\hat \rho_h- \tilde \rho_h) + \beta \int x (\tilde n_h-\hat n_h) (\hat \rho_s- \tilde \rho_s) \le 0. $

    This together with (4.11) yields

    $ \beta (\hat \rho_h-\tilde \rho_h)^2 + [\beta (1-b)+1](\hat \rho_h-\tilde \rho_h)(\hat \rho_s-\tilde \rho_s) + (\hat \rho_s-\tilde \rho_s)^2 \le 0. $

    Equivalently,

    $ (β(ˆρh˜ρh)+β(1b)+12β(ˆρs˜ρs))2+(1(β(1b)+1)24β)(ˆρs˜ρs)20. $ (4.13)

    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

    $ \int \left[ r_h(x)-\tilde \rho_h-\tilde \rho_s +(b-x) \int \tilde n_s \right]\hat n_h = \int \left[ r_h(x)-\hat \rho_h-\hat \rho_s +(b-x) \int \hat n_s \right] \tilde n_h, $
    $ [rs(y)˜ρh˜ρs+βx˜nh]ˆns=[rs(y)ˆρhˆρs+βxˆnh]˜ns, $ (4.14)

    respectively.

    (ⅰ) If $ r_s $ is strictly decreasing on $ [c, d] $, the fitness function $ G_s $ attains its global maximum at $ y = c $. Thus $ \hat n_s $ and $ \hat n_s $ either vanish or are the Dirac mass concentrated at $ y = c $. As a consequence we have $ \hat n_s = \tilde n_s = S \delta_{\{y = c\}} $ with $ S \ge 0 $. Using the formula for $ \hat n_s $, (4.5) and (4.6), we obtain (4.7)

    (ⅱ) The case $ \hat n_s = \tilde n_s = S \delta_{\{y = c\}} $ with $ S > 0 $. Identity (4.8) follows from (4.6).

    In the case $ \alpha \Sigma_h = \beta \Sigma_s = 0 $, an entropy functional has been found by Jabin and Raoul [16] 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 $ \alpha \Sigma_h, \beta \Sigma_s $. However, we obtain an entropy functional similar to that of Jabin and Raoul for the system (4.2) corresponding to particular choices of $ \alpha \Sigma_h $ and $ \beta \Sigma_s $. We obtain below a partial information about the dynamics of the solution of (4.2) as $ t \to \infty $: 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 $ (\hat n_h, \hat n_s) $ be an (non-negative) ESD of problem (4.2). Set

    $ E(t): = \beta \int n_h - \beta \int \hat n_h \ln (n_h)+ \int n_s - \int \hat n_s \ln (n_s). $

    Then E is an entropy functional for (4.2), i.e., $ \dfrac{d}{dt}E(t) \le 0 $ provide that $ \beta $ is a positive constant satisfying

    $ (\beta(1-b)+1)^2 \le 4 \beta. $

    Proof. Set

    $ E_1: = \int n_h - \int \hat n_h \ln (n_h). $

    We have

    $ ddtE1(t)=(nh)t(nhˆnh)nh=[rh(x)ρhρs+(bx)ns](nhˆnh)=[(rh(x)ρhρs+(bx)ns)(rh(x)ˆρhˆρs+(bx)ˆns)](nhˆnh)+[rh(x)ˆρhˆρs+(bx)ˆns)](nhˆnh)=(ρhˆρh)2(ρhˆρh)(ρsˆρs)+(bx)(nhˆnh)(nsˆns)+[rh(x)ˆρhˆρs+(bx)ˆns)]nh(ρhˆρh)2(ρhˆρh)(ρsˆρs)+(bx)(nhˆnh)(nsˆns). $

    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: = \int n_s - \int \hat n_s \ln (n_s), $ we have

    $ ddtE2(t)(ρsˆρs)2(ρhˆρh)(ρsˆρs)+βx(nhˆnh)(nsˆns). $

    Therefore, as $ E = \beta E_1 +E_2 $, we have

    $ \frac{dE}{dt} \le -\beta (\rho_h-\hat \rho_h)^2 -[\beta (1-b)+1](\rho_h-\hat \rho_h)(\rho_s-\hat \rho_s) -(\rho_s-\hat \rho_s)^2 \le 0, $

    where the last inequality follows from the negativity of the discriminant of this polynomial:

    $ (\beta(1-b)+1)^2 - 4 \beta \le 0. $

    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 [23] and perform the time change variable $ \tau = t \varepsilon $ to accelerate time and observe the dynamics. The parameter $ \varepsilon $ is also used to measure how close is the distribution from the Dirac distribution.

    The change of variable $ t \mapsto \tau $ converts the system (3.1) to

    $ {τnεh(τ,x)=1ε[rh(x)ρεh(τ)ρεs(τ)+α(x)Σs(τ)]nεh,x(a,b),τ>0,τnεs(τ,y)=1ε[rs(y)ρεh(τ)ρεs(τ)+β(y)Σh(τ)]nεs,y(c,d),τ>0. $

    This system is completed with the initial data

    $ n^ \varepsilon_h(0, x) = n^ \varepsilon_{h0}(x) \gt 0, \quad n^ \varepsilon_s(0, y) = n^ \varepsilon_{s0}(y) \gt 0. $

    Rather than working on $ n_h^ \varepsilon, n_s^ \varepsilon $ directly, we define as usual [19,21,23] the functions $ u^ \varepsilon, v^ \varepsilon $ given by

    $ uε(τ,x)=εln(nεh(τ,x)),uε0(x)=εln(nεh0(x)),vε(τ,y)=εln(nεs(τ,y)),vε0(y)=εln(nεs0(y)). $

    The functions $ u^ \varepsilon, v^ \varepsilon $ satisfy

    $ {τuε(τ,x)=rh(x)ρεh(τ)ρεs(τ)+α(x)ynεs(τ,y),x(a,b),τ>0,τvε(τ,y)=rs(y)ρεh(τ)ρεs(τ)+β(y)xnεh(τ,x),y(c,d),τ>0. $ (5.1)

    Our purpose is to study the behaviour of $ u^ \varepsilon, v^ \varepsilon $ as $ \varepsilon \to 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 $ \rho_{h0}^ \varepsilon+\rho_{s0}^ \varepsilon $ is uniformly bounded, Proposition 2.2 yields that there exists $ (\hat n_h, \hat n_s) $ such that as $ \varepsilon \to 0 $ (after extracting subsequences),

    $ nεhˆnh  in  L(0,;M([a,b])), $ (5.2)
    $ nεsˆns  in  L(0,;M([c,d])). $ (5.3)

    Theorem 5.1. Assume that $ r_h, \alpha \in C^1([a, b]), r_s, \beta \in C^1([c, d]) $ and that

    $ ρεh0+ρεs0+uε0C1([a,b])+vε0C1([c,d])K0. $ (5.4)

    Then

    (ⅰ) The function $ u^ \varepsilon $ (resp. $ v^ \varepsilon $) is uniformly Lipschitz continuous on $ [0, T] \times [a, b] $ (resp. $ [0, T] \times [c, d] $) for all $ T > 0 $.

    (ⅱ) As $ \varepsilon \to 0 $ (after extractions of subsequences), the functions $ u^ \varepsilon $ and $ v^ \varepsilon $ converge locally uniformly to Lipschitz continuous functions $ u $ and $ v $. Moreover, $ u, v $ satisfy

    $ {u(τ,x)=u0(x)+rh(x)ττ0ˆnhτ0ˆns+α(x)τ0yˆns,v(τ,y)=v0(y)+rs(y)ττ0ˆnhτ0ˆns+β(y)τ0xˆnh,maxτ,xu(τ,x)0,maxτ,yv(τ,y)0 for all τ0. $ (5.5)

    Furthermore we have for a.e. $ \tau $,

    $ {\rm{supp}}\ \hat n_h(\tau, \cdot) \subset \{u(\tau, \cdot) = 0\}, \quad {\rm{supp}}\ \hat n_s(\tau, \cdot) \subset \{v(\tau, \cdot) = 0\}. $

    Proof. (ⅰ) First note that by Proposition 2.2, there is a constant $ K_1 > 0 $ such that

    $ nεhL(0,;L1(a,b))+nεsL(0,;L1(c,d))K1. $ (5.6)

    Differentiating the equation for $ u^ \varepsilon $ with respect to $ x $ yields

    $ τuεx(τ,x)=rh(x)+α(x)ynεs(τ,y)dy. $ (5.7)

    Thus, using (5.6), we obtain $ |\partial_\tau u_x^ \varepsilon(\tau, x)| \le |\overline{r'}_h| +\overline{|\alpha'|} d K_1 $ so that

    $ |u_x^ \varepsilon(\tau, x)| \le K_0 + (|\overline{r'}_h| +\overline{|\alpha'|} d K_1)\tau. $

    On the other hand, in view of (5.1), we have

    $ |\partial_\tau u^ \varepsilon(\tau, x)| \le \bar r_h +2K_1+\alpha K_1 d. $

    Hence $ u^ \varepsilon $ is uniformly Lipschitz continuous on $ [0, T] \times [a, b] $. Similarly, the same property holds for $ v^ \varepsilon $.

    (ⅱ) Using the point (ⅰ) and the Arzelà–Ascoli Theorem, we may extract subsequences $ (u^ \varepsilon, v^ \varepsilon) $ 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 $ (\rho^ \varepsilon_h, \rho^ \varepsilon_s) $ blows up in the limit as $ \varepsilon $ vanishes and this is in contradiction with (5.6).

    We provide sufficient conditions so that $ (\hat n_h, \hat n_s) $ defined in (5.2), (5.3) is a monomorphic state.

    Theorem 5.2. Let all hypotheses as in Theorem 5.1 hold. Suppose furthermore that

    $ {uε0xxK,vε0yyK,rh(x)0,α(x)0,rs(y)0,β(y)0. $ (5.8)

    Then, in the distributional sense

    $ u_{xx} \le -K^*, \quad v_{yy} \le -K^*. $

    Thus for all $ \tau $, the functions $ u(\tau, \cdot), v(\tau, \cdot) $ are concave so that they have a unique maximum point. As a consequence, $ \hat n_h, \hat n_s $ have the form

    $ \hat n_h(\tau, x) = \hat \rho_h(\tau) \delta_{\{x = \hat x(\tau)\}}, \quad \hat n_s(\tau, y) = \hat \rho_s(\tau) \delta_{\{y = \hat y(\tau)\}}. $

    Moreover,

    $ If  ˆρh(τ)>0,  then  maxxu(τ,x)=u(τ,ˆx(τ))=0,If  ˆρs(τ)>0,  then  maxyv(τ,y)=v(τ,ˆy(τ))=0. $

    Proof. Differentiating twice the equation of $ u^ \varepsilon $, we obtain

    $ \partial_\tau u_{xx}^ \varepsilon(\tau, x) = r''_h(x) +\alpha''(x) \int y n^ \varepsilon_s(\tau, y) \le 0. $

    Thus $ u_{xx}^ \varepsilon(\tau, x) \le u_{xx}^ \varepsilon(\tau = 0, x) \le -K^* $. Therefore $ u_{xx} \le -K^* $. Similarly, $ v_{yy} \le -K^* $. Hence the theorem follows.

    In this section we derive the equations for the concentration point $ \hat x(\tau), \hat y(\tau) $. Our equations are valid until the time $ T^* $ where $ \hat \rho_h(\tau) > 0, \hat \rho_s(\tau) > 0 $ and that $ \hat x(\tau), \hat y(\tau) $ do not touch the boundary and that $ \dot{\hat x}, \dot{\hat y} $ are smooth enough (see Remark 5.3 below for the regularity of $ \hat x, \hat y $). For all $ \tau \in (0, T^*) $ we have $ \hat x(\tau) \in (a, b) $ is the maximum point of $ u(\tau, \cdot) $ on $ [a, b] $. It follows that $ u_x(\tau, \hat x(\tau)) = 0 $ so that

    $ u_{x\tau} (\tau, \hat x(\tau))+u_{xx}(\tau, \hat x(\tau)) \dot{\hat x} = 0 \ \text{ with } \ \dot{\hat x}: = d \hat x/d\tau. $

    This implies

    $ ˙ˆx(τ)=1uxx(τ,ˆx(τ))[rh(ˆx(τ))+α(ˆx(τ))yˆns]=1uxx(τ,ˆx(τ))[rh(ˆx(τ))+α(ˆx(τ))ˆy(τ)ˆρs(τ)]. $ (5.9)

    Similarly, we have

    $ ˙ˆy(τ)=1vyy(τ,ˆy(τ))[rs(ˆy(τ))+β(ˆy(τ))ˆx(τ)ˆρh(τ)]. $ (5.10)

    The equations (5.9) and (5.10) describe the dynamics of $ \hat x(\tau), \hat y(\tau) $. We also can obtain a more explicit form of (5.9) and (5.10) by representing $ \hat \rho_h, \hat \rho_s $ in terms of $ \hat x, \hat y $. To that purpose, we first notice that $ u (\tau, \hat x(\tau)) = 0 $ for $ \tau \in [0, T^*) $ so that

    $ 0 = \dfrac{d u (\tau, \hat x(\tau))}{d\tau} = u_\tau (\tau, \hat x(\tau)) = r_h(\hat x (\tau))-\hat \rho_h(\tau)-\hat \rho_s(\tau)+ \alpha(\hat x(\tau)) \hat y(\tau) \hat \rho_s (\tau). $

    Similarly,

    $ r_s(\hat y (\tau))-\hat \rho_h(\tau)-\hat \rho_s(\tau)+ \beta(\hat y(\tau) \hat x(\tau) \hat \rho_h (\tau) = 0. $

    Therefore, under the assumption that $ 1-(1-\alpha(\hat x(\tau)) \hat y(\tau))(1-\beta(\hat y(\tau))\hat x(\tau)) \neq 0 $, we have

    $ \hat \rho_h(\tau) = \dfrac{r_h(\hat x(\tau))-r_s(\hat y(\tau))(1-\alpha(\hat x(\tau)) \hat y(\tau))}{1-(1-\alpha(\hat x(\tau)) \hat y(\tau))(1-\beta(\hat y(\tau))\hat x(\tau))}, $
    $ \hat \rho_s(\tau) = \dfrac{r_s(\hat y(\tau))-r_h(\hat x(\tau))(1-\beta(\hat y(\tau))\hat x(\tau))}{1-(1-\alpha(\hat x(\tau)) \hat y(\tau))(1-\beta(\hat y(\tau))\hat x(\tau))}. $

    Substituting these expressions of $ \hat \rho_h(\tau), \hat \rho_s(\tau) $ in (5.9) and (5.10) we obtain

    Canonical equations

    $ {˙ˆx=1uxx(τ,ˆx)[rh(ˆx)+α(ˆx)ˆyrs(ˆy)rh(ˆx)(1β(ˆy)ˆx)1(1α(ˆx)ˆy)(1β(ˆy)ˆx)],˙ˆy=1vyy(τ,ˆy)[rs(ˆy)+β(ˆy)ˆxrh(ˆx)rs(ˆy)(1α(ˆx)ˆy)1(1α(ˆx)ˆy)(1β(ˆy)ˆx)]. $ (5.11)

    Remark 5.3. In view of the two first equations in (5.5), if $ u_0, v_0, r_h, \alpha, r_s, \beta $ are smooth enough (e.g., $ C^2 $ functions), then $ u_{xx} $ is continuous. Thus the right-hand-sides of (5.11) are continuous as functions of $ \hat x, \hat y, \tau $. Therefore, the standard ordinary differential equation theory implies that the solution of (5.11) $ (\hat x, \hat y) $ is a $ C^1 $ function of time. We refer to [26], for a results about the regularity of $ u, v $ in the case of Hamilton–Jacobi equations.

    Let us illustrate numerically the convergence of the solution of (3.1) towards an ESD as well as the movement of concentration points. We employ the two distinct sets $ (a, b): = (1, 2) $ and $ (c, d): = (3, 4) $ to insist on that fact that the two phenotypes of HSCs and of stromal cells are not the same. The space and time steps are given by

    $ \delta x = \delta y = 0.005, \qquad \delta t = 0.01. $

    Define for $ 0 \le k \le 200 $ and $ p = 0, 1, 2, ... $,

    $ x_k: = 1+k \delta x, \quad y_k: = 3+ k \delta y, \quad t_p: = p \delta t, \quad (n_h)^p_k: = n_h(t_p, x_k), \quad (n_s)^p_k: = n_h(t_p, y_k). $

    Our numerical simulations are performed in MATLAB and based on the implicit-explicit scheme below:

    $ {(nh)p+1k(nh)pkδt=max(0,Rh(xk,(ρh)p,(ρs)p,(Σs)p)(nh)pkmax(0,Rh(xk,(ρh)p,(ρs)p,(Σs)p)(nh)p+1k,(ns)p+1k(ns)pkδt=max(0,Rs(xk,(ρh)p,(ρs)p,(Σs)p)(ns)pkmax(0,Rs(yk,(ρh)p,(ρs)p,(Σh)p)(ns)p+1k,Rh(xk,(ρh)p,(ρs)p,(Σs)p):=rh(xk)(ρh)p(ρs)p+α(xk)(Σs)p,Rs(yk,(ρh)p,(ρs)p,(Σs)p):=rs(yk)(ρh)p(ρs)p+β(yk)(Σh)p. $ (6.1)

    Here $ (\rho_h)^p, (\rho_s)^p, (\Sigma_h)^p, (\Sigma_s)^p $ are the approximation of $ \rho_h, \rho_s, \Sigma_h, \Sigma_s $ at the time $ p $. We use the initial conditions

    $ n_{h0} = \exp(-(x-1.5)^2/0.01), \quad n_{s0} = \exp(-(y-3.4).^2/0.01) $

    for Figures 3, 4, 5 and 6. For the last figures, we choose

    $ n_{h0} = \exp(-(x-1.45)^2/0.002), \quad n_{s0} = \exp(-(y-3.58).^2/0.002). $
    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.  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.
    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.  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.

    The parameters are chosen as in the table below.

    Table 1.  Settings in the numerical simulations.
    Parameters Figures 3 and 4 Figures 5 and 6 Figures 7 and 8
    (monomorphic situation) (monomorphic situation) (dimorphic situation)
    healthy case leukemic case co-existence case
    $ r_h $ $ 0.1(x-1) $ $ x-1 $ $ 0.75(x-1)^2 $
    $ \alpha $ $ 0.1(2-x) $ $ 0.4(2-x) $ $ 0.625(2-x) $
    $ r_s $ $ 0.1 $ $ 0.6(4.5-y) $ $ 0.5+0.25(3-y) $
    $ \beta $ $ 0.2(y-3)+1 $ $ 0.1 $ $ 0.5 $

     | Show Table
    DownLoad: CSV

    Remark 6.1. The parameters in the three cases satisfy the assumptions of Theorem 4.2(ⅲ), (ⅱ) and Proposition 4.3(ⅲ), 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.

    Figures 3 and 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. Figures 5 and 6 below show another monomorphic situation where stromal cells with lowest support capacity are selected. Healthy HSCs cannot survive and cancer cells are selected.

    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(ⅲ). 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.

    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.  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.

    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:

    $ nεh(τ,x)ˆρh(τ)δ{x=ˆx(τ)} withˆρh(τ)>0maxxu(τ,x)=0=u(τ,ˆx(τ))nεs(τ,y)ˆρs(τ)δ{y=ˆy(τ)} withˆρs(τ)>0maxyv(τ,x)=0=v(τ,ˆy(τ)). $

    Here the points of concentration $ \hat x(\tau), \hat y(\tau) $ represent well adapted phenotypes at the time $ \tau $. Also these points are maximum points of the phase functions $ u(\cdot, \tau) $ and $ v(\cdot, \tau) $. The system (5.11) gives us the dynamics of $ \hat x(\tau), \hat y(\tau) $, 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\int\!(b-x)n_h(t, x)dx - k_2\int\!(x-a)n_h(t, x)dx $, with $ k_2 > k_1 $, instead of the neutral term $ - \rho_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 [22].

    Our analytic results, except in Section 4.2, hold for more general choices of the weight functions $ \psi_h, \psi_s $. However, we present our results here mainly for the case $ \psi_h = x, \psi_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, [24,28]) 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 $ -\infty $, 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.

    B.P. has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (grant agreement No 740623). The authors acknowledge funding by Plan Cancer HTE programm EcoAML.

    The authors declare there is no conflict of interest.

    [1] Elwood, P. ; Longley, M. (2010) My health: whose responsibility?: A Jury decides. J Epidemiol Commun H 64: 761-764. doi: 10.1136/jech.2009.087767
    [2] Antithrombotic, T. (2009) Aspirin in the primary and secondary prevention of vascular disease: collaborative meta-analysis of individual participant data from randomised trials. Lancet 373:1849-1860. doi: 10.1016/S0140-6736(09)60503-1
    [3] Elwood, P. ; Morgan, G. ; Fone, D. et al. (2011) Aspirin use in a south-Wales county. Br J Cardiol 18: 238-240.
    [4] Morgan, G. (2012) Rapid health impact assessment of aspirin promotion for the secondary prophylaxis of vascular events in Wales. Qual Prim Care 20: 299-301.
    [5] Elwood, P. ; Gallagher, A. M. ; Duthie, G. G. et al. (2008) Aspirin, salicylates and cancer. Lancet373: 1301-1309.
    [6] Rothwell, P. M; Wilson, M. ; Elwin, C. E. et al. (2010) Long-term effect of aspirin on colorectal cancer incidence and mortality: 20-year follow-up of five randomised trials. Lancet 376: 1741-1750. doi: 10.1016/S0140-6736(10)61543-7
    [7] Rothwell, P. M; Fowkes, F. G. R; Belch, J. F. F. et al. (2011) Effect of daily aspirin on long-term risk of death due to cancer: analysis of individual patient data from randomised trials. Lancet377: 31-41.
    [8] Burn, J. ; Gerdes, A. M. ; Macrae, F. et al. (2011) Long-term effect of aspirin on cancer risk in carriers of hereditary colorectal cancer: an analysis from the CAPP2 RCT. Lancet 378: 281-287.
    [9] Phillips, I. ; Langley, R. ; Gilbert, D. ; et al. (2013) Aspirin as a treatment for cancer. Clin Oncol (R Coll Radiol) 25: 333-335. doi: 10.1016/j.clon.2013.03.001
    [10] Morgan, G. ; Elwood, P. (2009) Could recommendations about aspirin prophylaxis enhance colorectal cancer screening programmes? Eur J Public Health 19: 576-577. doi: 10.1093/eurpub/ckp062
    [11] Morgan, G. (2009) Aspirin for the prevention of vascular events. Public Health 123: 787-788. doi: 10.1016/j.puhe.2009.10.007
    [12] Gorelick, P. B. ; Weisman, S. M. (2005) Risk of haemorrhagic stroke with aspirin use: an update. Stroke 36: 1801-1807. doi: 10.1161/01.STR.0000174189.81153.85
    [13] Hansson, L. ; Zanchetti, A. ; Carruthers, S. G. et al. (1998) Effects of intensive blood-pressure lowering and low-dose aspirin in patients with hypertension: principal results of the Hypertension Optimal Treatment (HOT) randomised trial. HOT Study Group. Lancet 351:1755-1762.
    [14] Elwood, P. ; Morgan, G. ; Brown, G. ; et al. (2005) Aspirin for all over 50? For. British Medical Journal 330: 1440-1441. doi: 10.1136/bmj.330.7505.1440
    [15] Baigent, C. ; (2005) Aspirin for all over 50? Against. British Medical Journal 330: 1441-1442.
    [16] Baron, J. A. (2003) A randomized trial of aspirin to prevent colorectal adenomas. N Engl J Med348:891-899.
    [17] Chan, A. T. ; Giconnucci, E. L. ; Schernhammer, E. S. et al. (2004) A prospective study of aspirin use and the risk for colorectal adenoma. Ann Intern Med 140: 157-166. doi: 10.7326/0003-4819-140-3-200402030-00006
    [18] Brenner, H. ; Tao, S. ; Haug, U. (2010) Low-dose aspirin use and performance of immunochemical fecal occult blood tests. JAMA 304: 2513-2520. doi: 10.1001/jama.2010.1773
    [19] Ajani, U. A. ; Ford, E. S. ; Greenland, K. J. (2006) Aspirin use among US adults: Behavioural Risk Factor Surveillance System. Amer J Prev Med 30: 74-77. doi: 10.1016/j.amepre.2005.08.042
    [20] Sung, J. J. ; Lau, J. Y. ; Ching, J. Y. et al. (2010) Continuation of low-dose aspirin therapy in peptic ulcer bleeding: a randomized trial. Ann Intern Med 152: 1-9. doi: 10.7326/0003-4819-152-1-201001050-00179
    [21] US Preventive Services Task Force. (2009) Aspirin for the prevention of cardiovascular disease: US preventive services task force recommendation statement. Ann Intern Med 150: 396-404. doi: 10.7326/0003-4819-150-6-200903170-00008
  • This article has been cited by:

    1. Guoning Si, Wenkai Li, Hanjing Lu, Zhuo Zhang, Xuping Zhang, Vibration Modeling and Analysis of a Flexible 3-PRR Planar Parallel Manipulator Based on Transfer Matrix Method for Multibody System, 2023, 11, 2075-1702, 505, 10.3390/machines11050505
    2. Yunxia Wei, Yuanfei Zhang, Bin Hang, Construction and management of smart campus: Anti-disturbance control of flexible manipulator based on PDE modeling, 2023, 20, 1551-0018, 14327, 10.3934/mbe.2023641
    3. Decio de M. Rinaldi, Tarcisio A. Hess-Coelho, Modular modeling of parallel mechanisms with actuation redundancy, 2025, 47, 1678-5878, 10.1007/s40430-025-05671-1
    4. Mikhail Drapalyuk, Aleksey Platonov, Petr Popikov, Modelling of the working process of the hydraulic pulsator of forest machine manipulator, 2025, 15, 2222-7962, 191, 10.34220/issn.2222-7962/2025.2/12
  • Reader Comments
  • © 2014 the Author(s), licensee AIMS Press. This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0)
通讯作者: 陈斌, bchen63@163.com
  • 1. 

    沈阳化工大学材料科学与工程学院 沈阳 110142

  1. 本站搜索
  2. 百度学术搜索
  3. 万方数据库搜索
  4. CNKI搜索

Metrics

Article views(4378) PDF downloads(889) Cited by(1)

Other Articles By Authors

/

DownLoad:  Full-Size Img  PowerPoint
Return
Return

Catalog