Loading [MathJax]/jax/output/SVG/jax.js
Research article Special Issues

Characterisation of the salmon cystic fibrosis transmembrane conductance regulator protein for structural studies

  • The cystic fibrosis transmembrane conductance regulator protein (CFTR) is a chloride channel highly expressed in the gills of Salmo salar, with a role in osmoregulation. It shares 60% identity with the human CFTR channel, mutations to which can cause the common genetic disorder cystic fibrosis CF. The expression and localisation of salmon CFTR have been investigated, but the isolated protein has not been extensively characterised. Here we present a protocol for the purification of recombinant salmon CFTR, along with biophysical and structural characterisation of the purified protein. Salmon CFTR was overexpressed in Saccharomyces cerevisiae, solubilised in the detergent LPG-14 and chromatographically purified by nickel-affinity and size-exclusion chromatography methods. Prior to size-exclusion chromatography samples of salmon CFTR had low purity, and contained large quantities of aggregated protein. Compared to size-exclusion chromatography profiles of other orthologues of CFTR, which had less evidence of aggregation, salmon CFTR appeared to have lower intrinsic stability than human and platypus CFTR. Nonetheless, repeated size-exclusion chromatography allowed monodisperse salmon CFTR to be isolated, and multi-angle light scattering was used to determine its oligomeric state. The monodispersity of the sample and its oligomeric state were confirmed using cryo-electron microscopy and small-angle X-ray scattering (SAXS). These data were also processed to calculate a low-resolution structure of the salmon CFTR, which showed similar architecture to other ATP-binding cassette proteins.

    Citation: Naomi L. Pollock, Oscar Moran, Debora Baroni, Olga Zegarra-Moran, Robert C. Ford. Characterisation of the salmon cystic fibrosis transmembrane conductance regulator protein for structural studies[J]. AIMS Molecular Science, 2014, 1(4): 141-161. doi: 10.3934/molsci.2014.4.141

    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 cystic fibrosis transmembrane conductance regulator protein (CFTR) is a chloride channel highly expressed in the gills of Salmo salar, with a role in osmoregulation. It shares 60% identity with the human CFTR channel, mutations to which can cause the common genetic disorder cystic fibrosis CF. The expression and localisation of salmon CFTR have been investigated, but the isolated protein has not been extensively characterised. Here we present a protocol for the purification of recombinant salmon CFTR, along with biophysical and structural characterisation of the purified protein. Salmon CFTR was overexpressed in Saccharomyces cerevisiae, solubilised in the detergent LPG-14 and chromatographically purified by nickel-affinity and size-exclusion chromatography methods. Prior to size-exclusion chromatography samples of salmon CFTR had low purity, and contained large quantities of aggregated protein. Compared to size-exclusion chromatography profiles of other orthologues of CFTR, which had less evidence of aggregation, salmon CFTR appeared to have lower intrinsic stability than human and platypus CFTR. Nonetheless, repeated size-exclusion chromatography allowed monodisperse salmon CFTR to be isolated, and multi-angle light scattering was used to determine its oligomeric state. The monodispersity of the sample and its oligomeric state were confirmed using cryo-electron microscopy and small-angle X-ray scattering (SAXS). These data were also processed to calculate a low-resolution structure of the salmon CFTR, which showed similar architecture to other ATP-binding cassette proteins.


    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] Plog S, Mundhenk L, Bothe MK, et al. (2010) Tissue and cellular expression patterns of porcine CFTR: similarities to and differences from human CFTR. J Histochem Cytochem 58: 785-797. doi: 10.1369/jhc.2010.955377
    [2] Crawford I, Maloney PC, Zeitlin PL, et al. (1991) Immunocytochemical localization of the cystic fibrosis gene product CFTR. P Natl Acad Sci USA 88: 9262-9266. doi: 10.1073/pnas.88.20.9262
    [3] Gadsby DC, Nairn AC (1999) Regulation of CFTR Cl- ion channels by phosphorylation and dephosphorylation. Adv Sec Messenger Phosphoprotein Res 33: 79-106. doi: 10.1016/S1040-7952(99)80006-8
    [4] Gadsby DC, Nairn AC (1999) Control of CFTR channel gating by phosphorylation and nucleotide hydrolysis. Physiol Rev 79: S77-S107.
    [5] Kirk KL, Wang W (2011) A unified view of cystic fibrosis transmembrane conductance regulator (CFTR) gating: combining the allosterism of a ligand-gated channel with the enzymatic activity of an ATP-binding cassette (ABC) transporter. J Biol Chem 286: 12813-12819. doi: 10.1074/jbc.R111.219634
    [6] Quinton PM, Reddy MM (2000) CFTR, a rectifying, non-rectifying anion channel? J Korean Med Sci 15 Suppl: S17-20.
    [7] Goss CH, Ratjen F (2013) Update in cystic fibrosis 2012. Am J Resp Crit Care 187: 915-919. doi: 10.1164/rccm.201301-0184UP
    [8] Welsh MJ, Ramsey BW (1998) Research on cystic fibrosis: a journey from the Heart House. Am J Resp Crit Care 157: S148-154. doi: 10.1164/ajrccm.157.4.nhlbi-13
    [9] Hiroi J, McCormick SD (2012) New insights into gill ionocyte and ion transporter function in euryhaline and diadromous fish. Resp Physiol Neurobi 184: 257-268. doi: 10.1016/j.resp.2012.07.019
    [10] Christensen AK, Hiroi J, Schultz ET, et al. (2012) Branchial ionocyte organization and ion-transport protein expression in juvenile alewives acclimated to freshwater or seawater. J Exp Biol 215: 642-652. doi: 10.1242/jeb.063057
    [11] Chen JM, Cutler C, Jacques C, et al. (2001) A combined analysis of the cystic fibrosis transmembrane conductance regulator: implications for structure and disease models. Mol Biol Evol 18: 1771-1788. doi: 10.1093/oxfordjournals.molbev.a003965
    [12] Kiilerich P, Kristiansen K, Madsen SS (2007) Cortisol regulation of ion transporter mRNA in Atlantic salmon gill and the effect of salinity on the signaling pathway. J Endocrinol 194: 417-427. doi: 10.1677/JOE-07-0185
    [13] Nilsen TO, Ebbesson LO, Madsen SS, et al. (2007) Differential expression of gill Na+, K+-ATPase alpha- and beta-subunits, Na+, K+, 2Cl- cotransporter and CFTR anion channel in juvenile anadromous and landlocked Atlantic salmon Salmo salar. J Exp Biol 210: 2885-2896. doi: 10.1242/jeb.002873
    [14] Mio K, Ogura T, Mio M, et al. (2008) Three-dimensional reconstruction of human cystic fibrosis transmembrane conductance regulator chloride channel revealed an ellipsoidal structure with orifices beneath the putative transmembrane domain. J Biol Chem 283: 30300-30310. doi: 10.1074/jbc.M803185200
    [15] Rosenberg MF, O'Ryan LP, Hughes G, et al. (2011) The cystic fibrosis transmembrane conductance regulator (CFTR): three-dimensional structure and localization of a channel gate. J Biol Chem 286: 42647-42654. doi: 10.1074/jbc.M111.292268
    [16] Zhang L, Aleksandrov LA, Riordan JR, et al. (2011) Domain location within the cystic fibrosis transmembrane conductance regulator protein investigated by electron microscopy and gold labelling. BBA-Biomembranes 1808: 399-404. doi: 10.1016/j.bbamem.2010.08.012
    [17] Awayn NH, Rosenberg MF, Kamis AB, et al. (2005) Crystallographic and single-particle analyses of native- and nucleotide-bound forms of the cystic fibrosis transmembrane conductance regulator (CFTR) protein. Biochem Soc T 33: 996-999. doi: 10.1042/BST20050996
    [18] Lewis HA, Buchanan SG, Burley SK, et al. (2004) Structure of nucleotide-binding domain 1 of the cystic fibrosis transmembrane conductance regulator. EMBO J 23: 282-293. doi: 10.1038/sj.emboj.7600040
    [19] Thibodeau PH, Brautigam CA, Machius M, et al. (2005) Side chain and backbone contributions of Phe508 to CFTR folding. Nat Struct Mol Biol 12: 10-16. doi: 10.1038/nsmb881
    [20] Galeno L, Galfre E, Moran O (2011) Small-angle X-ray scattering study of the ATP modulation of the structural features of the nucleotide binding domains of the CFTR in solution. Eur Biophys J 40: 811-824. doi: 10.1007/s00249-011-0692-5
    [21] Galfre E, Galeno L, Moran O (2012) A potentiator induces conformational changes on the recombinant CFTR nucleotide binding domains in solution. Cell Mol Life Sci 69: 3701-3713. doi: 10.1007/s00018-012-1049-7
    [22] Marasini C, Galeno L, Moran O (2013) A SAXS-based ensemble model of the native and phosphorylated regulatory domain of the CFTR. Cell Mol Life Sci 70: 923-933. doi: 10.1007/s00018-012-1172-5
    [23] Hudson RP, Chong PA, Protasevich, II, et al. (2012) Conformational changes relevant to channel activity and folding within the first nucleotide binding domain of the cystic fibrosis transmembrane conductance regulator. J Biol Chem 287: 28480-28494. doi: 10.1074/jbc.M112.371138
    [24] Huang P, Liu Q, Scarborough GA (1998) Lysophosphatidylglycerol: a novel effective detergent for solubilizing and purifying the cystic fibrosis transmembrane conductance regulator. Anal biochem 259: 89-97. doi: 10.1006/abio.1998.2633
    [25] Wiener MC (2004) A pedestrian guide to membrane protein crystallization. Methods 34: 364-372. doi: 10.1016/j.ymeth.2004.03.025
    [26] Carpenter EP, Beis K, Cameron AD, et al. (2008) Overcoming the challenges of membrane protein crystallography. Curr Opin Struc Biol 18: 581-586. doi: 10.1016/j.sbi.2008.07.001
    [27] Dobrovetsky E, Menendez J, Edwards AM, et al. (2007) A robust purification strategy to accelerate membrane proteomics. Methods 41: 381-387. doi: 10.1016/j.ymeth.2006.08.009
    [28] Granseth E, Seppala S, Rapp M, et al. (2007) Membrane protein structural biology--how far can the bugs take us? Mol Membr Biol 24: 329-332. doi: 10.1080/09687680701413882
    [29] Lewinson O, Lee AT, Rees DC (2008) The funnel approach to the precrystallization production of membrane proteins. J Mol Biol 377: 62-73. doi: 10.1016/j.jmb.2007.12.059
    [30] Graeslund S (2008) Protein production and purification. Nat Meth 5: 135-146. doi: 10.1038/nmeth.f.202
    [31] Mancia F, Love J (2010) High-throughput expression and purification of membrane proteins. J Struct Biol 172: 85-93. doi: 10.1016/j.jsb.2010.03.021
    [32] Aller SG, Yu J, Ward A, et al. (2009) Structure of P-glycoprotein reveals a molecular basis for poly-specific drug binding. Science 323: 1718-1722. doi: 10.1126/science.1168750
    [33] Kawate T, Gouaux E (2006) Fluorescence-detection size-exclusion chromatography for precrystallization screening of integral membrane proteins. Structure 14: 673-681. doi: 10.1016/j.str.2006.01.013
    [34] Sonoda Y, Cameron A, Newstead S, et al. (2010) Tricks of the trade used to accelerate high-resolution structure determination of membrane proteins. FEBS Lett 584: 2539-2547. doi: 10.1016/j.febslet.2010.04.015
    [35] Sonoda Y, Newstead S, Hu NJ, et al. (2011) Benchmarking membrane protein detergent stability for improving throughput of high-resolution X-ray structures. Structure 19: 17-25. doi: 10.1016/j.str.2010.12.001
    [36] Drew D, Newstead S, Sonoda Y, et al. (2008) GFP-based optimization scheme for the overexpression and purification of eukaryotic membrane proteins in Saccharomyces cerevisiae. Nat Protoc 3: 784-798. doi: 10.1038/nprot.2008.44
    [37] Newstead S, Kim H, von Heijne G, et al. (2007) High-throughput fluorescent-based optimization of eukaryotic membrane protein overexpression and purification in Saccharomyces cerevisiae. Proc Natl Acad Sci USA 104: 13936-13941. doi: 10.1073/pnas.0704546104
    [38] Clark KM, Fedoriw N, Robinson K, et al. Purification of transmembrane proteins from Saccharomyces cerevisiae for X-ray crystallography. Protein Expres Purif 71: 207-223.
    [39] Slotboom DJ, Duurkens RH, Olieman K, et al. (2008) Static light scattering to characterize membrane proteins in detergent solution. Methods 46: 73-82. doi: 10.1016/j.ymeth.2008.06.012
    [40] Ouano AC, Kaye W (1974) Gel-permeation chromatography: X. Molecular weight detection by low-angle laser light scattering. J Polym Sci: Polym Chem Edit 12: 1151-1162.
    [41] Miller JL, Tate CG (2011) Engineering an ultra-thermostable beta(1)-adrenoceptor. J Mol Biol 413: 628-638. doi: 10.1016/j.jmb.2011.08.057
    [42] Shibata Y, White JF, Serrano-Vega MJ, et al. (2009) Thermostabilization of the neurotensin receptor NTS1. J Mol Biol 390: 262-277. doi: 10.1016/j.jmb.2009.04.068
    [43] Tate CG, Schertler GF (2009) Engineering G protein-coupled receptors to facilitate their structure determination. Curr Opin Struc Biol 19: 386-395. doi: 10.1016/j.sbi.2009.07.004
    [44] Warne T, Serrano-Vega MJ, Tate CG, et al. (2009) Development and crystallization of a minimal thermostabilised G protein-coupled receptor. Protein Expres Purif 65: 204-213. doi: 10.1016/j.pep.2009.01.014
    [45] Aleksandrov AA, Kota P, Cui L, et al. (2012) Allosteric modulation balances thermodynamic stability and restores function of DeltaF508 CFTR. J Mol Biol 419: 41-60. doi: 10.1016/j.jmb.2012.03.001
    [46] Huang P, Stroffekova K, Cuppoletti J, et al. (1996) Functional expression of the cystic fibrosis transmembrane conductance regulator in yeast. Biochim Biophys Acta 1281: 80-90. doi: 10.1016/0005-2736(96)00032-6
    [47] Bear CE, Li CH, Kartner N, et al. (1992) Purification and functional reconstitution of the cystic fibrosis transmembrane conductance regulator (CFTR). Cell 68: 809-818. doi: 10.1016/0092-8674(92)90155-6
    [48] Kogan I, Ramjeesingh M, Li C, et al. (2002) Studies of the molecular basis for cystic fibrosis using purified reconstituted CFTR protein. Method Mol Med 70: 143-157.
    [49] Bai J, Swartz DJ, Protasevich, II, et al. (2011) A gene optimization strategy that enhances production of fully functional P-glycoprotein in Pichia pastoris. PloS One 6: e22577. doi: 10.1371/journal.pone.0022577
    [50] O'Ryan L, Rimington T, Cant N, et al. (2012) Expression and purification of the cystic fibrosis transmembrane conductance regulator protein in Saccharomyces cerevisiae. J Vis Exp e3860.
    [51] Pollock N, Cant N, Rimington T, et al. (2014) Purification of the cystic fibrosis transmembrane conductance regulator protein expressed in Saccharomyces cerevisiae. J Vis Exp e51447.
    [52] Schneider CA, Rasband WS, Eliceiri KW (2012) NIH Image to ImageJ: 25 years of image analysis. Nat Method 9: 671-675. doi: 10.1038/nmeth.2089
    [53] Sievers F, Wilm A, Dineen D, et al. (2011) Fast, scalable generation of high-quality protein multiple sequence alignments using Clustal Omega. Mol Syst Biol 7: 539.
    [54] Sievers F, Higgins DG (2014) Clustal Omega, accurate alignment of very large numbers of sequences. Methods in molecular biology 1079: 105-116. doi: 10.1007/978-1-62703-646-7_6
    [55] Dawson RJP, Locher KP (2007) Structure of the multidrug ABC transporter Sav1866 from Staphylococcus aureus in complex with AMP-PNP. FEBS Lett 581: 935-938. doi: 10.1016/j.febslet.2007.01.073
    [56] Hammersley A, Svensson S, Hanfland M, et al. (1996) Two-Dimensional Detector Software: From Real Detector to Idealised Image or Two-Theta Scan. High Pressure Res 14: 325-348.
    [57] Mateu L, Luzzati V, Vargas R, et al. (1990) Order-disorder phenomena in myelinated nerve sheaths. II. The structure of myelin in native and swollen rat sciatic nerves and in the course of myelinogenesis. J Mol Biol 215: 385-402.
    [58] Luzzati V, Tardieu A (1980) Recent developments in solution x-ray scattering. Annu Rev Biophys Bioeng 9: 1-29. doi: 10.1146/annurev.bb.09.060180.000245
    [59] Petoukhov M, Svergun D (2007) Analysis of X-ray and neutron scattering from biomacromolecular solutions. Curr Opin Struc Biol 17: 562-571. doi: 10.1016/j.sbi.2007.06.009
    [60] Guinier A, Fournet G (1955) Small angle scattering of x-rays. New York: Wiley.
    [61] Feigin L, Svergun D (1987) Structure analysis by small-angle x.ray and neutron scattering. New York, London: Plenum Press.
    [62] Svergun D (1992) Determination of the regularization parameter in indirect-transform methods using perceptual criteria. J Appl Crystallogr 25: 495-503. doi: 10.1107/S0021889892001663
    [63] Dawson RJ, Locher KP (2006) Structure of a bacterial multidrug ABC transporter. Nature 443: 180-185. doi: 10.1038/nature05155
    [64] Franke D, Svergun D (2009) DAMMIF, a program for rapid ab-initio shape determination in small-angle scattering. J Appl Crystallogr 42: 342-346. doi: 10.1107/S0021889809000338
    [65] Tian C, Vanoye CG, Kang C, et al. (2007) Preparation, functional characterization, and NMR studies of human KCNE1, a voltage-gated potassium channel accessory subunit associated with deafness and long QT syndrome. Biochemistry 46: 11459-11472. doi: 10.1021/bi700705j
    [66] Oliver RC, Lipfert J, Fox DA, et al. (2013) Dependence of micelle size and shape on detergent alkyl chain length and head group. PloS One 8: e62488. doi: 10.1371/journal.pone.0062488
    [67] Yang Z, Wang C, Zhou Q, et al. (2014) Membrane protein stability can be compromised by detergent interactions with the extramembranous soluble domains. Protein Sci 23: 769-789. doi: 10.1002/pro.2460
    [68] Gulati S, Jamshad M, Knowles TJ, et al. (2014) Detergent-free purification of ABC (ATP-binding-cassette) transporters. Biochem J 461: 269-278. doi: 10.1042/BJ20131477
    [69] Lyman CP (1968) Body temperature of exhausted salmon. Copeia 1968: 631-633. doi: 10.2307/1442045
    [70] Behrisch HW (1969) Temperature and the regulation of enzyme activity in poikilotherms. Fructose diphosphatase from migrating salmon. Biochem J 115: 687-696.
    [71] Handeland SO, Berge Ö, Björnsson BT, et al. (2000) Seawater adaptation by out-of-season Atlantic salmon (Salmo salar L.) smolts at different temperatures. Aquaculture 181: 377-396.
    [72] Hsu HH, Lin LY, Tseng YC, et al. (2014) A new model for fish ion regulation: identification of ionocytes in freshwater- and seawater-acclimated medaka (Oryzias latipes). Cell Tissue Res 357: 225-243. doi: 10.1007/s00441-014-1883-z
    [73] Moorman BP, Inokuchi M, Yamaguchi Y, et al. (2014) The osmoregulatory effects of rearing Mozambique tilapia in a tidally changing salinity. Gen Comp Endocrinol [in press].
    [74] Sucre E, Bossus M, Bodinier C, et al. (2013) Osmoregulatory response to low salinities in the European sea bass embryos: a multi-site approach. J Comp Physiol B 183: 83-97. doi: 10.1007/s00360-012-0687-2
    [75] Guggino WB, Stanton BA (2006) New insights into cystic fibrosis: molecular switches that regulate CFTR. Nat Rev Mol Cell Biol 7: 426-436. doi: 10.1038/nrm1949
    [76] Venerando A, Franchin C, Cant N, et al. (2013) Detection of phospho-sites generated by protein kinase CK2 in CFTR: mechanistic aspects of Thr1471 phosphorylation. PloS One 8: e74232. doi: 10.1371/journal.pone.0074232
  • 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(6762) PDF downloads(1283) Cited by(4)

Figures and Tables

Figures(7)  /  Tables(2)

/

DownLoad:  Full-Size Img  PowerPoint
Return
Return

Catalog