Dynamical properties of feedback signalling in B lymphopoiesis: A mathematical modelling approach
Salvador Chulián, Alvaro Martínez-Rubio, Anna Marciniak-Czochra, Thomas Stiehl, Cristina Blázquez Goñi, Juan Francisco Rodríguez Gutiérrez, Manuel Ramirez Orellana, Ana Castillo Robleda, Víctor M. Pérez-García, María Rosa
DDynamical properties of feedback signalling in Blymphopoiesis: A mathematical modelling approach
Salvador Chuli´an (cid:73) a,b , ´Alvaro Mart´ınez-Rubio (cid:73) a,b , Anna Marciniak-Czochra c ,Thomas Stiehl c , Cristina Bl´azquez Go˜ni d , Juan Francisco Rodr´ıguezGuti´errez d , Manuel Ramirez Orellana e , Ana Castillo Robleda e , V´ıctor M.P´erez-Garc´ıa f,g,h , Mar´ıa Rosa a,b a Department of Mathematics, Universidad de C´adiz, Puerto Real, C´adiz, Spain b Biomedical Research and Innovation Institute of C´adiz (INiBICA), Hospital UniversitarioPuerta del Mar, C´adiz, Spain c Institute of Applied Mathematics, BioQuant and Interdisciplinary Center of ScientificComputing (IWR), Heidelberg University, Heidelberg, Germany d Department of Paediatric Haematology and Oncology, Hospital de Jerez C´adiz, Spain e Department of Paediatric Haematology and Oncology, Hospital Infantil Universitario Ni˜noJes´us, Universidad Aut´onoma de Madrid, Madrid, Spain f Department of Mathematics, Mathematical Oncology Laboratory (MOLAB), Universidadde Castilla-La Mancha, Ciudad Real, Spain g Instituto de Matem´atica Aplicada a la Ciencia y la Ingenier´ıa (IMACI), Universidad deCastilla-La Mancha, Ciudad Real, Spain h ETSI Industriales, Universidad de Castilla-La Mancha, Ciudad Real, Spain
Abstract
Haematopoiesis is the process of generation of blood cells. Lymphopoiesis gen-erates lymphocytes, the cells in charge of the adaptive immune response. Dis-ruptions of this process are associated with diseases like leukaemia, which isespecially incident in children. The characteristics of self-regulation of this pro-cess make them suitable for a mathematical study.In this paper we develop mathematical models of lymphopoiesis using cur-rently available data. We do this by drawing inspiration from existing structuredmodels of cell lineage development and integrating them with paediatric bonemarrow data, with special focus on regulatory mechanisms. A formal analysis ofthe models is carried out, giving steady states and their stability conditions. Weuse this analysis to obtain biologically relevant regions of the parameter spaceand to understand the dynamical behaviour of B-cell renovation. Finally, we (cid:73)
These authors contributed equally to this work.
Preprint submitted to Journal of Theoretical Biology July 28, 2020 a r X i v : . [ q - b i o . T O ] J u l se numerical simulations to obtain further insight into the influence of prolif-eration and maturation rates on the reconstitution of the cells in the B line. Weconclude that a model including feedback regulation of cell proliferation repre-sents a biologically plausible depiction for B-cell reconstitution in bone marrow.Research into haematological disorders could benefit from a precise dynamicaldescription of B lymphopoiesis. Keywords:
Mathematical medicine, Haematopoiesis, Mathematical modelling,Lymphopoiesis
1. Introduction
Blood is a tissue under continuous regeneration, and its renovation is oneof the most studied developmental processes in biology [1]. It is initiatedby haematopoietic stem cells and develops through a multi-step differentia-tion cascade [2], resulting in the generation of red blood cells, platelets and cells of the immune system. Figure 1(a) shows the standard representation ofhaematopoiesis as a tree. At the top, stem cells with the potential for self-renewal give rise to respective lineage progenitors. These cells become progres-sively more specialised as they move towards the bottom of the tree. There aretwo major cell lineages: the myeloid line and the lymphoid line. The latter gen- erates lymphocytes, involved in adaptive immune response, which is responsiblefor ‘targeted’ reactions to infections.In this work, we will focus on the description of B lymphopoiesis, i.e. thedevelopment and maturation of B cells. These cells have a range of roles, beingmainly associated with the secretion of antibodies, the elements in charge of the neutralisation of foreign invaders [3]. Figure 1(b) shows a schematic repre-sentation of the route from common lymphoid progenitor to immature B cells,which eventually exit the bone marrow to complete maturation elsewhere. Al-terations in the generation of B cells are related to diseases like autoimmunereactions, immunodeficiencies or lymphoproliferative disorders like lymphomas
2r leukaemias [4]. The latter have especial incidence in children and constitutearound one third of all childhood cancer cases [5].
NKT lymphocyte B lymphocyteBasophilMonocyteThrombocyte Myeloid progenitor
Lymphoid precursor
Hematopoietic stem cell
Lymphoid progenitor
Pro-BPre-B LargePre-B SmallImmature B … (a) (b) Figure 1: (a) Schematic representation of haematopoiesis. (b) B-cell lineage starts from acommon lymphoid progenitor and then progresses towards immature B cells, which eventuallyleave the bone marrow. Hallmarks of this process are the acquisition of the immunophenotypiccell surface markers CD19 for the whole B line and CD10 for immature B cells.
Haematopoiesis is a perfect example of self-renewal and stemness in tissues[6]. Mathematical descriptions have been performed using multi-compartmental,continuously structured models [7, 8, 9]. In those models, each cell type in Fig- ure 1(a) can be thought of as a cell compartment that receives input from theupper elements and sends its output to the lower compartments [10, 11, 12,13, 14, 15, 39]. Some models have focused explicitly on the B-cell line [17, 18].Unlike full haematopoiesis, this process is sequential, simplifying its mathe-matical conceptualisation. The compartmental models described above become nonlinear when the interactions between the different cell stages are included.Indeed, this process requires some kind of negative regulatory feedback in orderto ensure steady production [19, 20]. This feature is common in the modellingof biological systems since they normally consist of a considerable number ofinteracting components [21]. In the case of B cells, their development depends on the joint action of anumber of factors that support or inhibit B-cell growth and differentiation [22].A clear description of the participants at each stage of development is lacking34], which is in part due to the difficulties in recreating in-vivo conditions inexperimental designs [23]. Mathematical models can help in elucidating which processes are more influenced by a given type of signal, as has been done for anumber of lineages and scenarios [24, 25, 26, 27, 28], and a better comprehensionof these interactions can be useful in understanding progression to malignancies[29].In this context, a complete picture of B lymphopoiesis would be a useful complement to those modelling scenarios that consider an input of B cells or thatrepresent B-cell-related phenomena [18, 30, 31, 32, 33]. Precise knowledge of thedynamical behaviour of each stage of differentiation stage can also be of helpin clinical situations where disorders are especially linked to the characteristicsof the cell of origin. In leukaemias, for instance, the phenotype of the tumour cells is an important diagnostic criterion [34]. Our aim in this paper is thus todevelop a model of B lymphopoiesis in the bone marrow. We will construct andinvestigate the properties and dynamical behaviour of a series of models. Thiswill be complemented with data from the literature, which mainly comes from in-vitro assays and immune reconstitution studies, and with clinical data from haematological patients.This paper is structured as follows: In Sec. 2 we explain basic haematologicalmodels, going from a general model to a reduced family of models more suitablefor lymphopoiesis. In Sec. 3 we perform a mathematical analysis of thesemodels, including positivity, boundedness and stability. In Sec. 4 we carry out numerical simulations, taking into account clinical data and information fromthe literature. In Sec. 5 we discuss these results and examine the potentialof these models to describe the process, concluding with the kind of researchopportunities that this analysis paves the way for.
2. Mathematical models and methods Our aim here is to describe B-cell development taking into account what thedata can tell us about the structure of the population of this haematopoietic4ine. An illustration of the representation of cell development can be found in[35, 36], where the authors proposed n maturation stages for a cell population u n = u n ( t ), with t ∈ R representing time. In this case, u would represent stem cell population, u n a mature specialised cell and u i ( i = 2 , . . . , n − p i = p i ( t ) and theso-called self renewal fraction a i = a i ( t ), for each maturation stage i = 1 , ..., n .The latter is considered to be the probability of a cell remaining in the samecell compartment after mitosis. The authors assumed that cells at stage i en-ter mitosis with a rate p i , resulting in a total number of 2 p i u i after mitosis.Then, with probability a i , they remain in the same compartment, whereas withprobability 1 − a i they go on to the next maturation stage. Therefore, each cellcompartment has an output of p i u i and an input 2 p i a i u i from their own com-partment. Consequently, the previous compartment (except at the first stage)provides the input corresponding to the number of cells that go on to the nextmaturation stage after mitosis: 2 p i − (1 − a i − ) u i − . Lastly, mature cells die ata rate d , and they are not considered to enter mitosis. If we consider n = 3stages for stem cells ( u ), intermediate cells ( u ) and specialised cells ( u ), theresult is the following system of equations: du dt = p (2 a − u , (1a) du dt = p (2 a − u + 2 p (1 − a ) u , (1b) du dt = 2 p (1 − a ) u − du . (1c)B cells are far from the haematopoietic stem cells since they are already com-mitted to the B line, thus losing part of their potential for self-renewal. We thenchoose to specify cell behaviour in each compartment in terms of proliferation and maturation, i.e. progression to the next stage, and to restrict the system tothree different cell compartments, considering the most common immunophe-notypical characterisation used in clinical practice [37]. The first compartmentwould also receive input from previous lymphoid progenitors. However, this5arly compartment is smaller and thus a constant source term contribution would be less significant [38]. Furthermore, this input is also regulated, whichwould require adding an equation from the previous compartment, and simi-larly for even earlier compartments. Our aim was to restrict the analysis to theCD19 + fraction of the B-cell line, as suggested by the data (see Appendix A).Thus we will consider three compartments accounting for the different matu-ration stages: early B cells ( C = C ( t )), intermediate B cells ( C = C ( t )), andfinally late B cells ( C = C ( t )), where t ∈ R represents time. A compartmentalmodel can then be written as dC dt = ρ C − α C , (2a) dC dt = ρ C + α C − α C , (2b) dC dt = α C − α C . (2c)Note that this formulation is equivalent to the model from Eq. (1) with u i = C i for i = 1 , , ρ i = p i , α i = 2 p i (1 − a i ) for i = 1 ,
2, and α = d .Early B cells, described here by Eq. (2a), have a proliferation rate ρ anda transition rate into intermediate B cells of α . Analogously, the proliferationrate for intermediate B cells and transition rate into the late compartment are ρ and α , as described in Eq. (2b). In this equation, a fraction of α C cellscomes from the early B compartment. This also happens with the α C cellsthat change their phenotypes from the intermediate B-cell into the late B-cellcompartment. Late B cells in Eq. (2c) are not considered to enter mitosis andthey go into the blood flow with a blood transition rate α . This compartmental model needs to be complemented with a regulatorysystem involving cell feedback signalling s ( t ). Different types of signalling andthe importance of the regulation of self-renewal in homeostatic (steady) statehave been previously studied in the literature [14, 19, 20, 24, 25, 26, 27]. In ourcase, there are different ways of specifying which cells participate in signalling and through which processes. In this paper we consider two different hypotheses.6he first is that signals are produced either by late cells (model 1) or by allcells (model 2). The second is that signals can alternatively affect either theproliferation rate (model A) or the transition rate of the model (model B).Therefore, we will consider four possibilities regarding feedback signalling. This is summarized in Figure 2. Model 1Model 2
Maturation stages
Last stage signallingAll stages signalling ρ i α i Model A
Proliferation rate af ected
Model B
Transition rate af ectedf (a) (b) f Figure 2: Representation of the different feedback signalling possibilities. (a) The signal maybe produced by the most mature cells and then affect the previous compartments (model 1);or it may be produced by all cells, influencing the whole population (model 2). (b) Signallingcan alternatively affect either the proliferation rate (model A) or the transition rate (modelB).
As stated in the introduction, the precise number and role of interactingelements is unclear. The basic immune system chemical messengers are thecytokines, proteins that control cell production. They can be generated bymicroenvironmental elements but also from cells themselves [40]. For B cells, a number of cytokines have been shown to be relevant: IL-7 for proliferation,differentiation (transition) and survival [41], and SCF and FLT3LG [42] for pro-liferation of early stages [43]. Evidence in this regard comes mainly from in-vitro and murine models and differences with humans can be significant [23]. Due tothis uncertainty and complexity we decided to follow an implicit formulation for the signalling as in [18], where instead of including each contributor explicitlywe consider the systemic action of each, and gather them together in a singlefunction s ( t ).With respect to the form of this signalling function s ( t ), we follow the de-velopment set out in [35]. First, let us consider a maximal signal ρ S , which is7elf-limited with rate α S . Consider then a number of cells N = N ( t ) that regu-lates the production of this signal in such a way that it decreases when there isa large number. Thus, a general signalling S = S ( t ) can be modelled as dSdt = ρ S − α S S − βSN, (3)where β is the inhibitory influence of cells N . By making a change of variables s ( t ) = α S /ρ S · S ( t ) and k = β/α S , we obtain a differential equation whosequasi-steady state is s ( t ) = 1 / (1 + kN ( t )). A rigorous analysis of this quasi-steady state approximation can be found in [44]. Signal concentration dependson the number of cells of type N . Following the explanation in Figure 2 we willconsider the case when only the late cells participate in signalling, i.e. s ( t ) = 11 + kC , (4)or the case where all cells participate s ( t ) = 11 + k (cid:80) i =1 C i . (5)Note that s ( t ) is a decreasing function of the cell numbers C i so it has an in-hibitory role. The parameter k measures the strength of the inhibitory feedback. To sum up, we will consider the models: dC dt = s ρ ρ C − s α α C , (6a) dC dt = s ρ ρ C + s α α C − s α α C , (6b) dC dt = s α α C − s α α C . (6c)Feedback signalling affects either proliferation (signal s ρ = s ρ ( t )) or transitionrates (signal s α = s α ( t )). The form of the signals s ρ , s α depends on the signallingsource, either s ( t ) or s ( t ). This yields four possible models, summarised inTable 1. 8ignalling 1. Late cells 2. All CellsA. Affecting proliferation Model A1 Model A2( s ρ = s , s α = 1) ( s ρ = s , s α = 1)B. Affecting transition Model B1 Model B2( s ρ = 1 , s α = s ) ( s ρ = 1 , s α = s ) Table 1: Mathematical models considered depending on signalling.
3. Theoretical results
Models A can be written as dC dt = ρ C kN − α C , (7a) dC dt = ρ C kN + α C − α C , (7b) dC dt = α C − α C , (7c)and models B have the form dC dt = ρ C − α C kN , (8a) dC dt = ρ C + α C kN − α C kN , (8b) dC dt = α C kN − α C kN . (8c)Feedback signalling can depend on cells from the late stage with N = C (modelsA1 and B1 with signal s ( t ) as in Eq. (4)) or on all cells with N = (cid:80) i =1 C i (models A2 and B2 with signal s ( t ) as in Eq. (5)). Theorem 1.
Let us consider the following set Q = { ( C , C , C ) ∈ R : C , C , C > } , (9)9 nd the initial values in QC ( t ) = C , C ( t ) = C , C ( t ) = C , (10) Then, the initial value problem for either Eqs. (7) or Eqs. (8) has a unique local-in-time solution for each t ∈ [ t − (cid:15), t + (cid:15) ] , for some value (cid:15) > . Proof.
The existence of a solution for systems from Eq. (7) and Eq. (8) isguaranteed for each ( C , C , C ) ∈ Q by continuity of the functions f ρ = f ρ ( C , C , C ) = ρ C kN − α C , (11a) f ρ = f ρ ( C , C , C ) = ρ C kN + α C − α C , (11b) f ρ = f ρ ( C , C , C ) = α C − α C , (11c)and f α = f α ( C , C , C ) = ρ C − α C kN , (12a) f α = f α ( C , C , C ) = ρ C + α C kN − α C kN , (12b) f α = f α ( C , C , C ) = α C kN − α C kN , (12c)where again N = C or N = (cid:80) i =1 C i . Boundedness of the respective partialderivatives of f αi and f ρi for i = 1 , , Henceforth we will consider that all parameters ρ i , α i and initial conditions C i are positive for i = 1 , , Theorem 2.
The solutions of Eqs. (7) and Eqs. (8) with ( C , C , C ) ∈ Q arepositive. Proof.
Let us consider the functions f ρi , f αi with i = 1 , , f ρi , f αi > − α i C i , we know that dC i dt > − α i C i . (13)10y integrating both sides of the equation from t to t we obtain C i ( t ) > C i exp( − α i t ) > . (14)And therefore all solutions C i ( t ), i = 1 , , definition. Theorem 3.
The solutions C ( t ) , C ( t ) , C ( t ) of Eqs. (7) with ( C , C , C ) ∈ Q are bounded. Proof.
In model from Eq. (7) if we consider C to be unbounded, thenlim C →∞ dC dt = ∞ , (15)which would imply that C would also be unbounded, and analogously for thecase of C . We could then writelim C i →∞ C i kC = lim C i →∞ C i k (cid:80) j =1 C j = 1 k for i = 1 , , . (16)Let us now consider the functions f ρi from Eq. (11). From Eq. (16) we get f ρ < ρ k − α C , (17)and then for all t C < Ae − α t + ρ kα , A ∈ R , (18)which yields that C ( t ) is bounded. Considering C ( t ) < M ∈ R , for all t , then f ρ < ρ k + α M − α C , (19)and integrating as above implies C is bounded. Considering C < M ∈ R , f ρ < α M − α C , (20)and thus C is also bounded.For the models B, ruled by Eqs. (8), we can sum the three equations toobtain dC T dt = ρ C + ρ C − α C kN , (21)11here C T = C + C + C . It is clear that for C and C sufficiently large the negative term makes only a small contribution and thus solutions are notbounded. This implies that models in which signalling affects only transitionrates are not appropriate for representing biological processes of this kind. Biological processes in homeostasis are stable and robust. In addition to being mathematically well posed, we need to ensure that the models have apositive stable equilibrium in which the three populations coexist. In this sectionwe study the existence of such states and their local stability. We focus onmodels A (Eq. (7)), since models B (Eq. (8)) do not lead to biologicallyrelevant dynamics. An expanded analysis of the stability conditions for certain steady states of models A can be found in Appendix B. Further analysis ofmodels B, showing that they have only unstable non-trivial positive equilibria,is presented in Appendix C.
Let us consider last stage signalling s ρ = s ( t ) as in Eq. (4) affecting theproliferation term, i.e. we study Eq. (7) with N = C . The three steady statesfor this model are P A = (0 , , , (22a) P A = (cid:18) , α ( ρ − α ) kα , ρ − α kα (cid:19) , (22b) P A = (cid:18) α ( ρ − α )( α ρ − α ρ ) kα α ρ , α ( ρ − α ) kα α , ρ − α kα (cid:19) . (22c)The Jacobian matrix of the system at any point ( C , C , C ) is J A ( C , C , C ) = ρ C k + 1 − α − C kρ ( C k + 1) α ρ C k + 1 − α − C kρ ( C k + 1) α − α . (23)12ubstituting P A = (0 , ,
0) in Eq. (23) we get the eigenvalues λ A , = − α , (24a) λ A , = ρ − α , (24b) λ A , = ρ − α . (24c)This is the trivial equilibrium that would be unstable in normal homeostaticprocesses. Instability conditions are ρ > α , or ρ > α . (25)As we consider P A >
0, it must be ρ > α , and therefore P A is unstable. For P A , we obtain λ A , = α ρ ρ − α , (26a) λ A , = − α − (cid:112) α (4 α − α ρ + α ρ )2 √ ρ , (26b) λ A , = − α (cid:112) α (4 α − α ρ + α ρ )2 √ ρ . (26c)This equilibrium point corresponds to a situation where the less differentiatedcompartment disappears and the system is maintained only by the proliferationof the second, leading to mature cells. As before, this is not a biologicallyfeasible situation, thus this equilibrium must be unstable. Since R ( λ , ) < R ( λ , ) <
0, then α ρ ρ − α >
0, which means that for P A to be unstablewe must have α ρ − α ρ > . (27)From the positivity of the non-trivial equilibrium point P A we require ρ > α , (28a) ρ ρ > α α . (28b)Conditions (28b) and (27) are identical which means that the existence of thispositive equilibrium implies the instability of P A . Stability conditions for P A are lengthy and can be found in Appendix B.13 .2.2. Model A2 Model A2 is given by Eqs. (7) with N = (cid:80) i =1 C i . The equilibria are P A = (0 , , , (29a) P A = (cid:18) , α ( ρ − α ) kα ( α + α ) , ρ − α k ( α + α ) (cid:19) , (29b) P A = (cid:18) α ( ρ − α )( α ρ − α ρ ) α kβ , α ρ ( ρ − α ) kβ , α ρ ( ρ − α ) kβ (cid:19) , (29c)where β = ( α α ρ + α ( α ρ + α ( ρ − ρ ))) . (30)The Jacobian matrix is J A = − α + ( s − − kC ) s ρ − C ks ρ − C ks ρ α − C ks ρ − α + ( s − − kC ) s ρ − C ks ρ α − α . (31)Substituting P A in Eq. (31) we obtain λ A , = − α , (32a) λ A , = ρ − α , (32b) λ A , = ρ − α . (32c)As before, for this equilibrium to be unstable we should have either ρ > α or/and ρ > α . From the positivity of P A , it must be the case that ρ > α .For P A , we obtain the eigenvalues λ A , = α ρ ρ − α , (33a) λ A , = α α − α α ρ − α ρ + h ( α , α , ρ )2 ρ ( α + α ) , (33b) λ A , = α α − α α ρ − α ρ − h ( α , α , ρ )2 ρ ( α + α ) . (33c)where h = h ( α , α , ρ ) such that h = √ α (cid:113) α α ( α − ρ ) ρ + 4 α ( α − ρ ) ρ + α ρ + α ( α + 4 ρ ) . (34)14s in the previous model, the existence of P A as a positive equilibriuminfluences the stability of P A . For the positivity of P A we have a first set ofconditions β > , (35a) ρ > α , (35b) α ρ > α ρ ; (35c)or a second set of conditions β < , (36a) ρ < α , (36b) α ρ < α ρ . (36c)Let us first consider that Eq. (35) holds. Then, in Eq. (33), the eigenvalue λ A , is positive and P A would be unstable. On the other hand, if Eq. (36) holds,then in Eq. (33) the eigenvalue λ A , <
0. Then P A would be stable whenever |R ( h ) | > α α − α α ρ − α ρ (37)for h as defined in Eq. (34).In this case, for the stability of P A , the difference with model A1 is theexistence of a denominator β (Eq. (30)). The stability conditions for P A are not shown here for reasons of space, but they are presented in Appendix B. Stability conditions related to model B, as well as a summary of all thestability conditions depending on both models, can be found in Appendix Cand Appendix D, respectively.
4. Numerical results
In the models presented above there are two parameters related to prolifer-ation ρ i ( i = 1 , α i ( i = 1 , , k .15 direct measure of the proliferation rates for the specific subsets consideredin this paper is lacking, but we can provide an estimation based on qualitative biological information. Normal B-cell development can be compared to datafrom autologous bone marrow transplantation. This type of transplantation ismore likely to reproduce developmental ontogeny, especially if there is no prioror coadjuvant anti-B-cell therapy [45]. In this case B-cell progenitors can bedetected in bone marrow as early as 1 month after transplantation [46, 47], and in blood after some delay [48]. In-vitro studies with mice show that proliferationrates are of the order of magnitude of days [40, 49, 50]. Human lymphoidcultures suggest doubling times of 1 day [51]. With respect to the relative values,current schemes for B-cell maturation indicate that upon CD19 acquisition thereis sustained proliferation that decreases as the cell matures [40, 52, 53]. Late
B cells already express B-cell receptor and a negative selection process occursprior to this [22, 54]. Based on this we can consider that the net productionrate of intermediate cells is lower than early cells ( ρ > ρ ).In order to obtain values for transition rates we make use of steady stateexpressions given by (22c) or (29c). These values can be compared to the flowcytometry data of normal bone marrow (See Appendix A). However, in oursimulations we are measuring absolute cell counts while flow cytometry is onlyable to measure relative cell proportions. The positive non-trivial equilibrium( C ∗ , C ∗ , C ∗ ) given by (22c) or (29c) satisfies the relationships C ∗ C ∗ = α α , (38a) C ∗ C ∗ = α ρ − α ρ α ρ . (38b)These quantities can be obtained from the analysis of the relative abundanceof the three different populations. Thus, for each blood transition rate α we can obtain values of α , α that agree with steady-state data and at the sametime belong in the positive stability region of parameter space. An exampleof the correlation found between transition rates is shown in Figure 3. Theorder of magnitude of α can be estimated as follows. The term α C ∗ gives16he total number of B cells per hour sent to blood by the bone marrow in homeostatic circumstances. In mice, bone marrow produces around 0 .
1% ofthe steady state population per hour [55]. In humans, the total steady-stateB-lymphocyte population can be obtained from lymphocyte proportions andbone-marrow volumes from the literature and we estimate it to be 10 cells[56, 57, 58, 59, 60]. This yields a B-cell production of 10 cells per hour, and given that C ∗ ∼ we estimate α to be of the order of magnitude of 10 − h − . Figure 3: Dependence of the transition rates over the range of stability. Lines represent α (light purple solid line), α (purple dashed line), α (dark purple dotted line). Across thisrange α j satisfies the relationship α > α > α . Values computed from Eq. (38) with ρ = 0 . − , ρ = 0 . − and steady-state cell proportions C ∗ /C ∗ = 0 . C ∗ /C ∗ = 3 . The last parameter to be determined is the inhibition constant k . Inspectionof steady states shows that, given transition and proliferation rates of the orderof magnitude of days, k is of the order of magnitude of the inverse of the totalsteady-state population k ∝ ( C ∗ + C ∗ + C ∗ ) − . The precise values for k and α are selected so that we recover steady-state values and reconstitution timescompatible with the literature cited above. Ranges of parameters in agreementwith positivity, stability and steady-state conditions are shown in Table 2.17aram. Meaning Value Units ρ Early B proliferation rate ln(2)/24 hours -1 ρ Intermediate B proliferation rate ln(2)/36 hours -1 α Transition rate: early to intermediate (0.004, 0.025) hours -1 α Transition rate: intermediate to late (0.003, 0.02) hours -1 α Transition rate: late to blood (0.01, 0.065) hours -1 k Inhibition constant (10 − , − ) cells -1 Table 2: Parameter values for Eqs. (7).
With respect to the initial state, the absolute number of mononucleatedtransplanted cells (MNC) is in the range of 10 cells [48]. From these cells, only a 1% are B early cells (10 cells) [38, 61]. These cells travel through bloodinto the bone marrow but only 10% of cells eventually reach the bone marrow[62]. Therefore, we will consider for the numerical simulations of autologoustransplantation an initial absolute number of cells of C (0) = 10 , C (0) = 0, C (0) = 0. The influence of this initial value on the dynamics of the system is described in Sec. 4.3. Typical results of simulations of models A1 and A2 (Eq. (7)) are shown inFigure 4. Recall that in model A1 (Figure 4(a)) the signal depends only on cellsin the final compartment while in model A2 (Figure 4(b)) all cells are involved.
Both models exhibit qualitatively similar behaviour. Early cells appear first,reaching a peak in population numbers slightly before day 30. Intermediatecells follow, reaching the respective peak with days of delay and with largercell numbers. Late cells appear last and stabilize in between. The systemsettles into the steady state from day 80 onwards. This behaviour agrees with the conditions expressed above for the parameter values. In particular, theproportion of population from each stage is coherent with clinical data (SeeAppendix A) and experimental data [38]. For both models we have 8.99%early cells, 70.21% intermediate cells and 20.80% late cells.18 igure 4: Comparison of numerical solutions of models A1 (a) and A2 (b) (Eq. (7)) andrelative proportions (c) and (d), respectively, during the first 100 days. Curves representearly cells C (dark blue solid line), intermediate cells C (blue dashed-dotted line) and latecells C (light blue dashed line). Both simulations have initial data C (0) = 10 , C (0) = 0, C (0) = 0 cells and parameter values ρ = 0 . − , ρ = 0 . − , α = 0 .
008 h − , α = 0 .
006 h − , α = 0 .
02 h − and k = 10 − . There are two main differences between the two models. The first relates to the magnitude of the steady state, which is larger when only late cells participatein signalling (model A1, Figure 4(a)) for the same parameter values. This canbe observed in Figure 5(a), where total cell numbers for both models are shown.Note that peak lymphocyte count, i.e. the largest cell number, occurs at day30, when intermediate cells are maximal. The second difference relates to the stage. It is interesting to relate this to absolute counts since, as explainedin Sec. 4.1, flow cytometry data only captures relative cell proportions. Weobserve very close behaviour between the models. The percentage of early cellsquickly decreases as more mature stages appear and steady-state proportionsare reached from day 80 onwards. Note that even though intermediate and mature cells have a peak in absolute cell count, this peak is not represented interms of percentage. Also, this figure shows that a decrease in cell percentagedoes not necessarily means a decrease in absolute cell count, something to takeinto account when interpreting longitudinal flow cytometry data.
Figure 5: Comparison of numerical simulations of model A1 (orange solid line) and A2 (reddashed line) from Eq. (7) for: (a) total number of cells C + C + C , with initial cellnumbers C (0) = 10 , C (0) = 0, C (0) = 0; (b) time to peak lymphocyte count as a functionof initial early B cell number. Both simulations have parameter values ρ = 0 . − , ρ = 0 . − , α = 0 .
008 h − , α = 0 .
006 h − , α = 0 .
02 h − and k = 10 − . .3. Time to peak decreases exponentially with initial value In Sec. 4.1 we described the rationale for the choice of the initial value ofearly cells. Despite this, we sought to determine how the scale of this dataimpacted reconstitution times. In Figure 5(b) we show the time to peak cellcount for a range of initial values for early cells. There exists a decreasing ex-ponential relationship between the two magnitudes, although the delay is not significant when considering that the literature on reconstitution after autolo-gous transplantation describes reconstitution times in the range of 1-2 months[46]. Multiplying initial cell numbers by 10 results in a displacement in time of5 days.
Clinical data suggests that homeostatic bone marrow displays relatively con-stant subset proportions (See Appendix A). This, together with the analysis ofthe expressions of the steady states allowed us to derive a connection betweenthe three transition rates in the model (Eq. 38). For the ranges of parametersconsidered, we observed that α > α > α (see Figure 3), which suggests that the second compartment being more numerous could be due not to a higherproliferation rate but rather to a slower maturation time. This calls for theanalysis of the influence of transition rates on the dynamics of the system.In order to do this we focused on variations of α , the rate at which latecells exit bone marrow and enter the blood flow. We select a range of variation that lies in the positive stability region and observe the qualitative differencesin the immune reconstitution. Results are shown in Figure 6(a).We observe that an increasing blood transition rate means that the systemreaches a lower number of total cells. Indeed, Eq. (22) shows that populationlevels depend on transition rates. Note that the peak during reconstitution also decreases. In order to quantify this reduction, we show, in Figure 6(b),the proportion of the height of the peak with respect to the final steady state,for the same range of values for α . For lower values of blood transition rate,the transitory population of lymphocytes can be 1 . igure 6: Numerical simulations for variable blood transition rate. (a) Total cell number C + C + C for t ∈ [0 , α = 0 .
03 h − (solid line), α = 0 .
04 h − (dashed line), α = 0 .
05 h − (dashed-dotted line), α = 0 .
06 h − (dottedline). Line colour goes from dark red to pink as blood transition rate increases. (b) Peakto steady state value ratio for the same range of blood transition rates. Both simulationsbelong to model A2 (Eq. (7)) with initial state C (0) = 10 , C (0) = 0, C (0) = 0 cells andparameter values ρ = 0 . − , ρ = 0 . − and k = 10 − . Parameters α and α vary with α according to Eq. (38). homeostatic conditions. This proportion decreases if cells increase the rate at which they join blood flow. Another consequence coming from these numericalsimulations is that the lower the transition rates, the longer it takes for thesystem to stabilise. Eq. (22) shows that steady state values depend on transition rates (see also
Figure 6) but also on the inhibition constant k . To quantify the impact of theparameter k on the dynamics of the system, we computed the total steady-statecell number as well as the proportion of peak height to steady state in a rangeof α (and thus α and α ) and k values. The results are shown in Figure 7 formodel A2. With respect to the absolute final cell count, only low levels of both k and α result in a much higher number of cells. We highlighted an area forwhich the total cell amount C T is in the range (10 , ), as estimated in Sec.22.1 from reference values. For the proportion of peak height to steady state,there is little variation in the direction of k , so signalling intensity has littleinfluence on the dynamics. As shown in Figure 6, high peak values belong in the low blood transition rate area. We conclude that transition rates primarilycause the overshoot during reconstitution, while k is mainly responsible for theexistence of stability regimes and the size of the final states. Cell number(a) Peak/steady state(b)
Figure 7: Influence of inhibition constant k and blood transition rate α parameters in modelA2 (Eq. (7) with all cell signalling) for (a) total steady state cell number, (b) peak vs steadystates ratio. The dashed highlighted area represents the range in which total cell number C T ∈ (10 , ). Each coordinate is the result of a numerical simulation of model A2 withinitial state cell numbers C (0) = 10 , C (0) = 0, C (0) = 0 and parameter values ρ =0 . − and ρ = 0 . − . Parameters α and α vary with α ∈ (0 .
01 h − , .
03 h − )as explained in Eq. (38).
5. Discussion and conclusions
Haematopoiesis is one of the most widely studied biological developmental processes [1]. Interesting questions arise related to the processes of cell lineagespecification [2, 63], the role of stem cells [6, 64] and the way cells communi-cate to regulate and ensure steady production [65, 66]. This is also true forB lymphopoiesis, the branch of haematopoiesis pertaining to B-cell formation.Specific unknowns in B-cell biology are the origins of some developmental stages, the role of senescence or the array of cytokines that regulates this process [67].23tudies of human B lymphopoiesis are encouraged, given that most of our knowl-edge about this line comes from mouse models [4]. Answering these questionsand obtaining a precise description of the dynamics of B-cell maturation arefundamental for research avenues. Some examples are the characterisation of haematological malignancies, the reconstitution after bone marrow transplantor chemotherapy or the generation of new lymphocytes from human embryonicstem cells [67].Mathematical models have the potential to integrate biological hypothesesand clinical data in order to provide an abstract representation of biological processes [11]. This representation can be useful not only for understandingthe dynamical properties of the system, but also for testing and elucidatingmore inaccessible phenomena. Our goal in this paper was to establish a biolog-ically sensible mathematical characterization of the B lymphopoiesis. Spurredby previous models of haematopoietic processes [14, 35], we designed four mod- els each with three differentiation stages. We added an implicit and systemicconsideration of cell feedback signalling resulting in four nonlinear models. Wefirst analysed these models from a theoretical perspective, addressing existence,positivity, boundedness and local stability. We collected data from the liter-ature and clinical data from haematological patients and then used numerical simulations in order to understand the role and influence of each parameter.We learn from theoretical analysis that a stable, homeostatic situation can-not arise solely by regulating the transition rate, i.e. the process of differen-tiation or maturation to the next compartment. We have focused from thereonward on feedback regulation of cell proliferation. Inspection of the steady states allowed us to use flow cytometry data to establish a relationship betweenthe three transition rates, analysing their influence by manipulating them oneby one. The relationship suggests that, for the specific cell proportions in theB line, cells transition faster from early to intermediate than from intermediateto late. Intermediate cells could then be more numerous not because they pro- liferate more, but because of their slower transition rates downstream towardsmore mature cell types. Also, numerical simulations show that cell proportion24s independent of which cells perform the signalling.In this sense, we observed that signalling coming from all stages results ina smoother reconstitution of the B cell line. Indeed, when all populations par- ticipate in signalling, their influence occurs earlier and proliferation decreasesfaster with time than when only late cells do. In this case we observed a peakthat we understood as a consequence of the delay in the reaction of the systemto overpopulation. The amplitude of this peak is also correlated with lowertransition rates, which implies lower steady states values. This is biologically understandable: the late compartment saturates due to excessive input fromprevious compartments, delaying access to stability and thus maintaining cellproduction in upstream compartments. It is important to remark here thatsubset percentage, a common metric in follow-up samples in a clinical context,can be misleading when dynamics of expansion are at play. For example, the overshoot during early reconstitution is not observed in terms of relative pro-portions. Finally, we noted that the strength of the signalling has no impacton the dynamics. The feedback loop could then be understood as a mecha-nism for the existence of stable output, while dynamical characteristics (time toreconstitution or early peak) are more dependent on intrinsic cellular processes.
The idea of this study was to determine which conditions are sufficient, froma mathematical perspective, to represent the kind of biological data that is cur-rently available for B-line development. While we obtain a behaviour that fitswith the time scales of the in-vivo process [45, 46, 47, 48, 51, 38], our study has aseries of limitations. Firstly, the choice of three compartments could be refined or expanded following a more detailed characterisation of the cells. Multidi-mensional flow cytometry data shows that surface markers, those that specifyto which stage a given cell belongs, vary continuously. A mathematical modelwhere these markers vary continuously might be able to capture this variation.Secondly, we described signalling as a systemic phenomenon. While this was enough to recapitulate known B-cell behaviour, a more detailed description in-cluding two or more types of signalling is desirable. Lastly, the model wouldbenefit from longitudinal data coming from immune reconstitution of the B-cell25ine. In this regard, flow cytometry analyses of both peripheral blood and bonemarrow in routine follow-up would allow for a more precise parameterisation and enable the hypotheses presented above to be contrasted.To conclude, we have constructed and studied several non-linear compart-mental models describing B cell lymphocyte reconstitution. These simple mod-els describe the process of B-cell generation as portrayed by bone marrow data,and we consider it a first step in a deeper exploration of the phenomena associ- ated with B-cell development. We verified mathematical and biological consis-tence, opening the door to interesting mathematical research like the existenceof bifurcations or the conditions for global stability, something that finds im-mediate application in cases of immune reconstitution. Studies of this kind canfunction as a source of hypothesis generation in biomedical research, for example when contrasting mouse versus human dynamics. Ultimately, we aim to extendthe methodology to situations of stability disruption and abnormal growth likeB-cell disorders and other haematological diseases.
Acknowledgements
The support of the Junta de Andaluc´ıa group FQM-201 is gratefully ac- knowledged. This work has been partially supported by the Fundaci´on Espa˜nolapara la Ciencia y la Tecnolog´ıa (FECYT, project PR214 from the University ofC´adiz) and the Asociaci´on Pablo Ugarte (APU). This work has been partiallysupported by the Junta de Comunidades de Castilla-La Mancha (grant numberSBPLY/17/180501/000154). ppendix A. Data Appendix A.1. Patients
Bone marrow samples from six individuals of paediatric age (1 to 13 years)were used to estimate cell subset proportions. Four patients diagnosed withIdiopathic Thrombocytopenic Purpura (1 from Jerez Hospital and 3 from Ni˜no
Jes´us Hospital) and two patients with neutropenia (Jerez Hospital). Due to thedifficulty in obtaining healthy bone marrow samples from patients of paediatricage, we selected the above as surrogate examples of normal B-cell development.Bone marrow samples were extracted from these patients in order to checkfor more severe disorders, but they were later diagnosed with B cell-unrelated diseases. Sample inspection further ensured lack of B-line affection. Instancesof this can be found in the literature [68, 69].
Appendix A.2. Flow Cytometry
Bone marrow samples were analysed by flow cytometry. This techniquemeasures the expression of the immunophenotypic markers that characterise each maturation stage. Marker expression for both hospitals’ data was acquiredon a FACSCanto II flow cytometer following the manufacturer’s specificationsBecton Dickinson (BD) for sample preparation. Samples were stained using an8-colour panel consisting of the following six fluorochrome-conjugated antibodiesprovided by BD: CD38 FITC/ CD10 PE/ CD34 PerCP-Cy5-5/ CD19 PE-Cy7/
CD20 APC/ CD45 V-450. This panel allows for the identification of B cellsubpopulations [34]. Forward (FSC) and side scatter (SSC) were also measured.Samples were preprocessed removing debris, doublets and marginal eventsas is routinely done in clinical and computational flow cytometry [70]. CD19 + cells were then gated in order to select B lymphocytes [71]. Since the model consists of three B cell populations, we performed k-means clustering on eachsample, including all B cell markers, with 3 predefined clusters. The algorithmwas initialised randomly and 100 random sets were generated, selecting the onewith lower within-cluster variation. Following standard immunophenotyping of27 cells [37], clusters were classified into early (CD45 - /CD10 + ), intermediate (CD45 + /CD10 + ) or late (CD45 + /CD10 - ) B cells. Proportions were then com-puted with respect to the total B lymphocyte count (CD19 + population), whichcorrelates with experimental data [38]. Figure A.8 shows the three stages of theprocess. All computations were carried out in RStudio using packages flowCore[72] and flowPeaks [73]. Figure A.8: (a) B Lymphocyte selection in grey (CD19 + population). Y-axis representscellular complexity, which is low for lymphocytes, and X-axis represents B Lymphocyte surfaceantigen CD19. (b) Within B population, three clusters were found corresponding to threematuration stages: early, intermediate and late populations. As B cells mature, they gainexpression of CD45 antigen and lose expression of CD10 antigen [37]. (c) Boxplots withproportions of each maturation stage from 6 patients displaying mean and 1st and 3rd quartile.Mean values with standard deviations: 0 . ± . , . ± . , . ± . Appendix B. Stability analysis for non-trivial equilibria in model A
Appendix B.1. Model A1
We recall the model from Section 3.2.1. From Eq. (22) we obtained thesteady states P A i for i = 1 , , P A . We obtain the char-acteristic equation λ + b λ + b λ + b = 0 , (B.1)28here b = α + α − α ρ ρ , (B.2a) b = α (cid:18) α − α ρ ρ (cid:19) , (B.2b) b = α α ( α − ρ )( α ρ − α ρ ) ρ . (B.2c)Using the Routh-Hurwitz Criterion, for P A to be positive and stable we musthave b b − b > , b > , b > . (B.3)The positivity conditions found in Eq. (28) yield b > , b >
0. Furthermore,the stability condition b b − b > ρ ≤ ρ . If ρ > ρ , thestability criterion is equivalent to satisfying either α ρ ≤ ρ ( α − α ρ + α ρ ) (B.4)or α ρ (cid:0) α + ρ ( α + α ) (cid:1) + α ρ α ρ > α ρ ( α T − ρ ) + α ρ ( ρ + ρ ) , (B.5)where α T = (cid:80) i =1 α i . Appendix B.2. Model A2
We recall the model from Section 3.2.2. We obtained in Eq. (29) the steadystates P A i for i = 1 , , P A . Weobtain the characteristic equation λ + b λ + b λ + b = 0 , (B.6)29here b = ( α α ( α + α ) ρ + α ρ ( α ρ + α α (3 ρ − ρ ) ρ β + (B.7a)+ α ( ρ − ρ )) − α ( α ( ρ − ρ ) ρ + α ρ ( α + ρ )) ρ β ,b = α (cid:0) α α ρ + 2 α α ρ ( α ρ + α ( ρ − ρ )) − α ρ ( ρ − ρ )( α + ρ ) (cid:1) ρ β +(B.7b)+ α α ( ρ − ρ ) ρ ρ β − α α ρ ( α ρ + α ( ρ − ρ ) ρ + α ρ ( α − ρ + 2 ρ )) ρ β ,b = α α ( α − ρ )( α ρ − α ρ ) ρ . (B.7c)Positivity conditions for P A as in Eq. (35) and in Eq. (36) yield b > b > b b − b > α α ( α + α ) ρ + α ρ ( α ρ + α α (3 ρ − ρ ) ρ β ++ α ( ρ − ρ )) − α ( α ( ρ − ρ ) ρ + α ρ ( α + ρ )) ρ β > r ρ + r ρ + r ( ρ − ρ ))( r ρ + r ρ + r ( ρ − ρ ))) β ++ α ( ρ − α ) ρ ( α ρ − α ρ ) > r = α α + 3 α α α + α α ( α + α ) ,r = − ( α α + 2 α α α ) ρ − α α α ,r = α α ( α ρ − α ρ ) ,r = α α + 2 α α + α α , (B.8c) r = − α α − α α α − α α ρ ,r = 2 α α α ρ + α ρ − ρ α ( α α + ρ ( α + α )) . ppendix C. Stability analysis for models B Appendix C.1. Model B1
Let us consider Eqs. (8) with last stage signalling s α = s ( t ) as in Eq. (4),this is, N = C . The steady states for this model are P B = (0 , , , (C.1a) P B = (cid:18) , α ( α − ρ ) kα ρ , α − ρ kρ (cid:19) , (C.1b) P B = (cid:18) α ( α − ρ )( α ρ − α ρ ) kα α ρ , α ( α − ρ ) kα ρ , α − ρ kρ (cid:19) . (C.1c)The Jacobian matrix is J B ( C , C , C ) = J B such that J B = ρ − α C k + 1 0 α C k ( C k + 1) α C k + 1 ρ − α C k + 1 α C k − α C k ( C k + 1) α C k + 1 − α C k + α ( C k + 1) . (C.2)Substituting P B i for i = 1 , , P B , we obtain the same eigenvalues as for P A , i.e. λ B , = − α , (C.3a) λ B , = ρ − α , (C.3b) λ B , = ρ − α . (C.3c)For P B we get λ B , = ρ − α ρ α , (C.4a) λ B , = − α ρ + √ α ρ √ α + α − ρ α , (C.4b) λ B , = α ρ − √ α ρ √ α + α − ρ α . (C.4c)Finally, for P B , we obtain the characteristic equation λ + b λ + b λ + b = 0 , (C.5)31here b = α ρ + α ρ − α ρ α , (C.6a) b = α ρ ( α ρ + ρ ( ρ − α )) α , (C.6b) b = α ( α − ρ ) ρ ( α ρ − α ρ ) α . (C.6c)Considering positivity conditions for P B and P B , we find that λ B ,i < i = 1 , ,
3, and therefore P B is always stable. From the positivity conditionswe also get ρ ρ > α α . (C.7a)which implies λ B , > P B is unstable.Stability of this equilibrium P B can be analysed by the Routh-HurwitzCriterion from Eq. (B.3). However, given its own positivity conditions, we get b <
0, implying P B is always unstable. Appendix C.2. Model B2
Let us now consider Eqs. (8) with signalling coming from all cellular com-partments s α = s ( t ) as in Eq. (5), this is, N = (cid:80) i =1 C i . The steady states ofthe model are P B = (0 , , , (C.8a) P B = (cid:18) , α ( α − ρ ) k ( α + α ) ρ , α ( α − ρ ) k ( α + α ) ρ (cid:19) , (C.8b) P B = (cid:18) α ( α − ρ )( α ρ − α ρ ) ρ kβ , α α ( α − ρ ) kβ , α α ( α − ρ ) kβ (cid:19) , (C.8c)where β = ( α α ρ + α ( α ρ + α ( ρ − ρ ))) . (C.9)The Jacobian matrix of Eqs. (8) with signal s given by Eq. (5) is J B = J B ( C , C , C ) such that J B = s C kα − α s + ρ s C kα C kα α + kR kR − α s + ρ s kR kR α + kR − α − kR , (C.10)32here R = C ( α + α ) + C α , (C.11a) R = C α − C α , (C.11b) R = C α − C α , (C.11c) R = C α + C ( α + α ) , (C.11d) R = C α + C ( α + α ) . (C.11e)Substituting P B i for i = 1 , , P B , we again get λ B , = − α , (C.12a) λ B , = ρ − α , (C.12b) λ B , = ρ − α . (C.12c)For P B , we get the eigenvalues λ B , = ρ − α ρ α , (C.13a) λ B , = − α ρ + α ρ + h ∗ ( α , α , ρ )2 α ( α + α ) , (C.13b) λ B , = − α ρ − α ρ + h ∗ ( α , α , ρ )2 α ( α + α ) . (C.13c)where h ∗ = h ∗ ( α , α , ρ ) such that h ∗ = √ α ρ (cid:113) α ( α + α (2 α − ρ ) + α ( α − ρ )) + α ( α − ρ ) . (C.14)Finally, for P B , we obtain the characteristic equation λ + b λ + b λ + b = 0 , (C.15)33here b = α α ρ ( α + α + ρ ) + α ρ ( α ρ + α α ( ρ − ρ ) α β + (C.16a)+ α ( ρ − ρ )) − α ( α ρ + α ( ρ − ρ )) ρ α β ,b = α ρ ( α ρ ( α ρ + α ( α + ρ )) + α ( ρ − ρ ) ρ α β + (C.16b) − α α ρ ( ρ − α ρ − ρ ρ ) − α ( α ρ + ( α + ρ )( ρ − ρ ) ρ ) α β ,b = α ( α − ρ ) ρ ( α ρ − α ρ ) α . (C.16c)Every equilibrium stability is influenced by the positivity conditions of theother points. From the positivity of P B , we get that α > ρ . (C.17)Two different scenarios arise from the positivity conditions of P B ; either β > , (C.18a) α > ρ , (C.18b) α ρ > α ρ ; (C.18c)or β < , (C.19a) α < ρ , (C.19b) α ρ < α ρ . (C.19c)Whenever Eq. (C.18) holds, equilibrium P B is stable (mainly α < ρ , asEq.(C.17) is true whenever P B > P B would alsobe stable whenever Eq. (C.19) holds and also |R ( h ∗ ) | < α ρ + α ρ (C.20)where h ∗ is defined as in Eq. (C.14). However, the stability of P B and P B is biologically uninteresting. Focusing on the non-trivial state, with the aboveconstraints Eq.(C.18) or Eq.(C.19) and the Routh-Hurtwitz criterion, we have b <
0. Therefore, P B is positive but always unstable.34 ppendix D. Summary of stability conditions We summarise in Table D.3 the conclusions of the mathematical analysisregarding stability of the non-trivial state.Steady Model A1 Model A2 Model B1 Model B2State s ρ = s , s α = 1 s ρ = s , s α = 1 s ρ = 1 , s α = s s ρ = 1 , s α = s P j Unstable Unstable Stable Conditionallystable P j Unstable Conditionally Unstable Conditionallystable stable P j Conditionally Conditionally Unstable UnstableStable Stable
Table D.3: Steady-state stability for every model from Eq. (6) under conditions of positivityof the non-trivial steady state. Index j stands for the four models considering the differentfeedback regulations: A A B B References
References [1] S. Doulatov, F. Notta, E. Laurenti, J. E. Dick, Hematopoiesis: a humanperspective, Cell stem cell 10 (2) (2012) 120–136. doi:10.1016/j.stem.2012.01.006 .[2] E. Laurenti, B. G¨ottgens, From haematopoietic stem cells to complex dif- ferentiation landscapes, Nature 553 (7689) (2018) 418–426. doi:10.1038/nature25022 .[3] B. Alberts, D. Bray, J. Lewis, M. Raff, K. Roberts, J. Watson, Molec-ular Biology of the Cell, 4th Edition, Garland, 2002. doi:10.1016/0307-4412(94)90059-0 . doi:10.1182/blood-2008-02-078071 .[5] E. Steliarova-Foucher, M. Colombet, L. A. Ries, F. Moreno, A. Dolya,F. Bray, P. Hesseling, H. Y. Shin, C. A. Stiller, S. Bouzbid, et al., Interna- tional incidence of childhood cancer, 2001–10: a population-based registrystudy, The Lancet Oncology 18 (6) (2017) 719–731. doi:10.1002/cncr.20910 .[6] T. Reya, S. J. Morrison, M. F. Clarke, I. L. Weissman, Stem cells, cancer,and cancer stem cells, Nature 414 (6859) (2001) 105–111. doi:10.1038/ .[7] G. Clapp, D. Levy, A review of mathematical models for leukemia andlymphoma, Drug Discovery Today: Disease Models 16 (2015) 1–6. doi:10.1016/j.ddmod.2014.10.002 .[8] Pujo-Menjouet, L., Blood cell dynamics: Half of a century of modelling, Math. Model. Nat. Phenom. 11 (1) (2016) 92–115. doi:10.1051/mmnp/201611106 .[9] T. Lorenzi, A. Marciniak-Czochra, T. Stiehl, A structured populationmodel of clonal selection in acute leukemias with multiple maturationstages, Journal of Mathematical Biology 79 (5) (2019) 1587–1621. doi: .[10] A. Marciniak-Czochra, T. Stiehl, W. Wagner, Modeling of replicativesenescence in hematopoietic development, Aging 1 (8) (2009) 723–732. doi:10.18632/aging.100072 .[11] I. Roeder, Quantitative stem cell biology: computational studies in the hematopoietic system, Current Opinion in Hematology 13 (4) (2006) 222–228. doi:10.1097/01.moh.0000231418.08031.48 .3612] S. Viswanathan, P. W. Zandstra, Towards predictive models of stemcell fate, Cytotechnology 41 (2/3) (2003) 75–92. doi:10.1023/a:1024866504538 . [13] V. V. Ganusov, J. Auerbach, Mathematical modeling reveals kinetics oflymphocyte recirculation in the whole organism, PLoS Computational Bi-ology 10 (5) (2014) e1003586. doi:10.1371/journal.pcbi.1003586 .[14] T. Stiehl, A. Marciniak-Czochra, Characterization of stem cells using math-ematical models of multistage cell lineages, Mathematical and Computer Modelling 53 (7-8) (2011) 1505–1517. doi:10.1016/j.mcm.2010.03.057 .[15] M. C. Mackey, R. Rudnicki, Global stability in a delayed partial differentialequation describing cellular replication, Journal of Mathematical Biology33 (1) (1994) 89–109. doi:10.1007/bf00160175 .[16] D. Dingli, F. Michor, Successful therapy must eradicate cancer stem cells, Stem Cells 24 (12) (2006) 2603–2610. doi:10.1634/stemcells.2006-0136 .URL https://doi.org/10.1634/stemcells.2006-0136 [17] G. Shahaf, S. Zisman-Rozen, D. Benhamou, D. Melamed, R. Mehr, Bcell development in the bone marrow is regulated by homeostatic feed- back exerted by mature b cells, Frontiers in immunology 7 (2016) 77. doi:10.3389/fimmu.2016.00077 .[18] S. Hu, O. A. Smirnova, F. A. Cucinotta, A biomathematical model oflymphopoiesis following severe radiation accidents-potential use for doseassessment, Health physics 102 (4) (2012) 425–436. doi:10.1097/HP. .[19] N. L. Komarova, Principles of regulation of self-renewing cell lineages, PloSone 8 (9). doi:10.1371/journal.pone.0072847 .[20] E. Manesso, J. Teles, D. Bryder, C. Peterson, Dynamical modelling ofhaematopoiesis: an integrated view over the system in homeostasis and doi:10.1098/rsif.2012.0817 .[21] D. S. Jones, M. Plank, B. D. Sleeman, Differential equations andmathematical biology, Chapman and Hall/CRC, 2009. doi:10.1201/9781420083583 . [22] K. Murphy, C. Weaver, Janeway’s immunobiology, Garland science, 2016. doi:10.1201/9781315533247 .[23] T. W. LeBien, Fates of human b-cell precursors, Blood, The Journal of theAmerican Society of Hematology 96 (1) (2000) 9–23. doi:10.1182/blood.V96.1.9 . [24] D. H. Fuertinger, F. Kappel, S. Thijssen, N. W. Levin, P. Kotanko, Amodel of erythropoiesis in adults with sufficient iron availability, Jour-nal of mathematical biology 66 (6) (2013) 1209–1240. doi:10.1007/s00285-012-0530-0 .[25] T. Walenda, T. Stiehl, H. Braun, J. Fr¨obel, A. D. Ho, T. Schroeder, T. W. Goecke, B. Rath, U. Germing, A. Marciniak-Czochra, et al., Feedback sig-nals in myelodysplastic syndromes: increased self-renewal of the malignantclone suppresses normal hematopoiesis, PLoS computational biology 10 (4). doi:10.1371/journal.pcbi.1003599 .[26] T. Stiehl, N. Baran, A. D. Ho, A. Marciniak-Czochra, Clonal selection and therapy resistance in acute leukaemias: mathematical modelling explainsdifferent proliferation patterns at diagnosis and relapse, Journal of TheRoyal Society Interface 11 (94) (2014) 20140079–20140079. doi:10.1098/rsif.2014.0079 .[27] W. Wang, T. Stiehl, S. Raffel, V. T. Hoang, I. Hoffmann, L. Poisa-Beiro,
B. R. Saeed, R. Blume, L. Manta, V. Eckstein, et al., Reduced hematopoi-etic stem cell frequency predicts outcome in acute myeloid leukemia,38aematologica 102 (9) (2017) 1567–1577. doi:10.3324/haematol.2016.163584 .[28] F. Knauer, T. Stiehl, A. Marciniak-Czochra, Oscillations in a white blood cell production model with multiple differentiation stages, Jour-nal of Mathematical Biology 80 (3) (2020) 575–600. doi:10.1007/s00285-019-01432-6 .[29] H. Jumaa, R. W. Hendriks, M. Reth, B cell signaling and tumorigen-esis, Annu. Rev. Immunol. 23 (2005) 415–445. doi:10.1146/annurev. immunol.23.021704.115606 .[30] R. Mehr, G. Shahaf, A. Sah, M. Cancro, Asynchronous differentiation mod-els explain bone marrow labeling kinetics and predict reflux between thepre-and immature b cell pools, International immunology 15 (3) (2003)301–312. doi:10.1093/intimm/dxg025 . [31] O. A. Smirnova, S. Hu, F. A. Cucinotta, Analysis of the lymphocytopoiesisdynamics in nonirradiated and irradiated humans: a modeling approach,Radiation Research 181 (3) (2014) 240–250. doi:10.1667/RR13256.1 .[32] R. Mostolizadeh, Z. Afsharnezhad, A. Marciniak-Czochra, Mathematicalmodel of chimeric anti-gene receptor (car) t cell therapy with presence of cytokine, Numerical Algebra, Control & Optimization 8 (1) (2018) 63. doi:10.3934/naco.2018004 .[33] O. Le´on-Triana, S. Sabir, G. F. Calvo, J. Belmonte-Beitia, S. Chuli´an,´A. Mart´ınez-Rubio, M. Rosa, A. P´erez-Mart´ınez, M. R. Orellana, V. M.P´erez-Garc´ıa, Car t cell therapy in b-cell acute lymphoblastic leukaemia: Insights from mathematical models, arXiv preprint. arXiv:2003.10236 .[34] J. Van Dongen, L. Lhermitte, S. B¨ottcher, J. Almeida, V. Van der Velden,J. Flores-Montero, A. Rawstron, V. Asnafi, Q. Lecrevisse, P. Lucio, et al.,Euroflow antibody panels for standardized n-dimensional flow cytomet-39ic immunophenotyping of normal, reactive and malignant leukocytes,
Leukemia 26 (9) (2012) 1908–1975. doi:10.1038/leu.2012.120 .[35] A. Marciniak-Czochra, T. Stiehl, A. D. Ho, W. Jger, W. Wagner, Modelingof asymmetric cell division in hematopoietic stem cells—regulation of self-renewal is essential for efficient repopulation, Stem Cells and Development18 (3) (2009) 377–386. doi:10.1089/scd.2008.0143 . [36] T. Stiehl, A. Marciniak-Czochra, Stem cell self-renewal in regeneration andcancer: insights from mathematical modeling, Current Opinion in SystemsBiology 5 (2017) 112–120. doi:10.1016/j.coisb.2017.09.006 .[37] E. van Lochem, V. van der Velden, H. Wind, J. te Marvelde, N. West-erdaal, J. van Dongen, Immunophenotypic differentiation patterns of nor- mal hematopoiesis in human bone marrow: Reference patterns for age-related changes and disease-induced shifts, Cytometry 60B (1) (2004) 1–13. doi:10.1002/cyto.b.20008 .[38] P. L´ucio, A. Parreira, M. van den Beemd, E. van Lochem, E. van Wering,E. Baars, A. Porwit-MacDonald, E. Bjorklund, G. Gaipa, A. Biondi, A. Or- fao, G. Janossy, J. van Dongen, J. S. Miguel, Flow cytometric analysis ofnormal b cell differentiation: a frame of reference for the detection of min-imal residual disease in precursor-b-all, Leukemia 13 (3) (1999) 419–427. doi:10.1038/sj.leu.2401279 .[39] D. Dingli, J. M. Pacheco, Modeling the architecture and dynamics of hematopoiesis, Wiley Interdisciplinary Reviews: Systems Biology andMedicine 2 (2) (2009) 235–244. doi:10.1002/wsbm.56 .[40] H. Kraus, S. Kaiser, K. Aumann, P. B¨onelt, U. Salzer, D. Vestweber,M. Erlacher, M. Kunze, M. Burger, K. Pieper, H. Sic, A. Rolink, H. Eibel,M. Rizzi, A feeder-free differentiation system identifies autonomously pro- liferating b cell precursors in human bone marrow, The Journal of Im-munology 192 (3) (2014) 1044–1054. doi:10.4049/jimmunol.1301815 .4041] R. Maddaly, G. Pai, S. Balaji, P. Sivaramakrishnan, L. Srinivasan, S. S.Sunder, S. F. Paul, Receptors and signaling mechanisms for b-lymphocyteactivation, proliferation and differentiation–insights from both in vivo and in vitro approaches, FEBS letters 584 (24) (2010) 4883–4894. doi:10.1016/j.febslet.2010.08.022 .[42] J. J. O’Shea, M. Gadina, R. M. Siegel, Cytokines and cytokine recep-tors, in: Clinical immunology, Elsevier, 2019, pp. 127–155. doi:10.1016/B978-0-7020-6896-6.00009-0 . [43] G. Petkau, M. Turner, Signalling circuits that direct early b-cell devel-opment, Biochemical Journal 476 (5) (2019) 769–778. doi:10.1042/BCJ20180565 .[44] A. Marciniak-Czochra, A. Mikeli´c, T. Stiehl, Renormalization groupsecond-order approximation for singularly perturbed nonlinear ordinary dif- ferential equations, Mathematical Methods in the Applied Sciences 41 (14)(2018) 5691–5710. doi:10.1002/mma.5107 .[45] T. Guillaume, D. B. Rubinstein, M. Symann, Immune reconstitution andimmunotherapy after autologous hematopoietic stem cell transplantation,Blood, The Journal of the American Society of Hematology 92 (5) (1998) doi:10.1182/blood.V92.5.1471 .[46] D. Leitenberg, J. M. Rappeport, B. R. Smith, B-cell precursor bone mar-row reconstitution after bone marrow transplantation, American journal ofclinical pathology 102 (2) (1994) 231–236. doi:10.1093/ajcp/102.2.231 .[47] J. E. Talmadge, E. Reed, K. Ino, A. Kessinger, C. Kuszynski, D. Heimann, M. Varney, J. Jackson, J. M. Vose, P. J. Bierman, Rapid immunologic re-constitution following transplantation with mobilized peripheral blood stemcells as compared to bone marrow, Bone marrow transplantation 19 (2)(1997) 161–172. doi:10.1038/sj.bmt.1700626 .4148] A. Parrado, S. Casares, J. Prieto, M. Carmona, A. Vaquero, J. Rodriguez-
Fernandez, Repopulation of circulating t, b and nk lymphocytes followingbone marrow and blood stem cell transplantation, Hematology and celltherapy 39 (6) (1997) 301–306. doi:10.1007/s00282-997-0301-3 .[49] G. J. Deenen, I. Van Balen, D. Opstelten, In rat b lymphocyte genesis sixtypercent is lost from the bone marrow at the transition of nondividing pre-b cell to sigm+ b lymphocyte, the stage of ig light chain gene expression,European journal of immunology 20 (3) (1990) 557–564. doi:10.1002/eji.1830200315 .[50] Y.-H. Park, D. G. Osmond, Dynamics of early b lymphocyte precursorcells in mouse bone marrow: proliferation of cells containing terminal de- oxynucleotidyl transferase, European journal of immunology 19 (11) (1989)2139–2144. doi:10.1002/eji.1830191125 .[51] H. E. Skipper, S. Perry, Kinetics of normal and leukemic leukocyte pop-ulations and relevance to chemotherapy, Cancer Research 30 (6) (1970)1883–1897. [52] A. K. Abbas, A. H. Lichtman, S. Pillai, Cellular and molecular immunology,Elsevier Health Sciences, 1994.[53] S. C. Bendall, K. L. Davis, E.-a. D. Amir, M. D. Tadmor, E. F. Simonds,T. J. Chen, D. K. Shenfeld, G. P. Nolan, D. Peer, Single-cell trajectorydetection uncovers progression and regulatory coordination in human b cell development, Cell 157 (3) (2014) 714–725. doi:10.1016/j.cell.2014.04.005 .[54] J. G. Monroe, G. Bannish, E. M. Fuentes-Panana, L. B. King, P. C.Sandel, J. Chung, R. Sater, Positive and negative selection during blymphocyte development, Immunologic research 27 (2-3) (2003) 427–442. doi:10.1385/IR:27:2-3:427 . 4255] D. G. Osmond, Population dynamics of bone marrow b lymphocytes, Im-munological reviews 93 (1) (1986) 103–124. doi:10.1111/j.1600-065x.1986.tb01504.x .[56] P. Clark, D. Normansell, D. Innes, C. Hess, Lymphocyte subsets in normal bone marrow, Blood 67 (6) (1986) 1600–1606. doi:10.1182/blood.V67.6.1600.1600 .[57] C. Andreoni, D. Rigal, M. Bonnard, J. Bernaud, Phenotypic analysis ofa large number of normal human bone marrow sample by flow cytometry,Blut 61 (5) (1990) 271–277. doi:10.1007/BF01732876 . [58] C. W. Caldwell, E. Poje, M. A. Helikson, B-cell precursors in normal pe-diatric bone marrow, American journal of clinical pathology 95 (6) (1991)816–823. doi:10.1093/ajcp/95.6.816 .[59] E. M. Rego, A. B. Garcia, S. R. Viana, R. P. Falc˜ao, Age-related changes oflymphocyte subsets in normal bone marrow biopsies, Cytometry: The Jour- nal of the International Society for Analytical Cytology 34 (1) (1998) 22–29. doi:10.1002/(SICI)1097-0320(19980215)34:1<22::AID-CYTO4>3.0.CO;2-G .[60] C. Nombela-Arrieta, M. G. Manz, Quantification and three-dimensionalmicroanatomical organization of the bone marrow, Blood Advances 1 (6) (2017) 407–416. doi:10.1182/bloodadvances.2016003194 .[61] Kleiveland, Charlotte R., Peripheral Blood Mononuclear Cells, SpringerInternational Publishing, Cham, 2015, pp. 161–167. doi:10.1007/978-3-319-16104-4_15 .[62] G. Caocci, M. Greco, G. La Nasa, Bone marrow homing and engraftment defects of human hematopoietic stem and progenitor cells, Mediterraneanjournal of hematology and infectious diseases 9 (1). doi:10.4084/MJHID.2017.032 . 4363] H. Kawamoto, H. Wada, Y. Katsura, A revised scheme for developmentalpathways of hematopoietic cells: the myeloid-based model, International immunology 22 (2) (2010) 65–70. doi:10.1093/intimm/dxp125 .[64] A. Wilson, E. Laurenti, G. Oser, R. C. van der Wath, W. Blanco-Bose,M. Jaworski, S. Offner, C. F. Dunant, L. Eshkind, E. Bockamp, et al.,Hematopoietic stem cells reversibly switch from dormancy to self-renewalduring homeostasis and repair, Cell 135 (6) (2008) 1118–1129. doi:10. .[65] L. Biasco, D. Pellin, S. Scala, F. Dionisio, L. Basso-Ricci, L. Leonardelli,S. Scaramuzza, C. Baricordi, F. Ferrua, M. P. Cicalese, et al., In vivotracking of human hematopoiesis reveals patterns of clonal dynamics duringearly and steady-state reconstitution phases, Cell stem cell 19 (1) (2016) doi:10.1016/j.stem.2016.04.016 .[66] K. Busch, K. Klapproth, M. Barile, M. Flossdorf, T. Holland-Letz, S. M.Schlenner, M. Reth, T. H¨ofer, H.-R. Rodewald, Fundamental propertiesof unperturbed haematopoiesis from stem cells in vivo, Nature 518 (7540)(2015) 542–546. doi:10.1038/nature14242 . [67] R. R. Hardy, P. W. Kincade, K. Dorshkind, The protean nature of cells inthe b lymphocyte lineage, Immunity 26 (6) (2007) 703–714. doi:10.1016/j.immuni.2007.05.013 .[68] Z. Ahmad, N. Durrani, T. Hazir, Bone marrow examination in itp in chil-dren: is it mandatory?, Journal of the College of Physicians and Surgeons– Pakistan: JCPSP 17 (6) (2007) 347–349. doi:06.2007/jcpsp.347349 .[69] H. Zafar, S. Anwar, M. Faizan, S. Riaz, Clinical features and outcome inpaediatric newly diagnosed immune thrombocytopenic purpura in a ter-tiary care centre, Pakistan journal of medical sciences 34 (5) (2018) 1195. doi:doi:10.12669/pjms.345.15687 . doi:10.1038/nri.2016.56 .[71] G. Finak, M. Langweiler, M. Jaimes, M. Malek, J. Taghiyar, Y. Korin,K. Raddassi, L. Devine, G. Obermoser, M. L. Pekalski, et al., Standard- izing flow cytometry immunophenotyping analysis from the human im-munophenotyping consortium, Scientific reports 6 (1) (2016) 1–11. doi:10.1038/srep20686 .[72] F. Hahne, N. LeMeur, R. R. Brinkman, B. Ellis, P. Haaland, D. Sarkar,J. Spidlen, E. Strain, R. Gentleman, flowcore: a bioconductor package for high throughput flow cytometry, BMC bioinformatics 10 (1) (2009) 106. doi:10.1186/1471-2105-10-106 .[73] Y. Ge, S. C. Sealfon, flowpeaks: a fast unsupervised clustering for flow cy-tometry data via k-means and density peak finding, Bioinformatics 28 (15)(2012) 2052–2058. doi:10.1093/bioinformatics/bts300 .725