Modelling Erythroblastic Islands: Using a Hybrid Model to Assess the Function of Central Macrophage
Stephan Fischer, Polina Kurbatova, Nikolai Bessonov, Olivier Gandrillon, Vitaly Volpert, Fabien Crauste
aa r X i v : . [ q - b i o . Q M ] M a y Modelling Erythroblastic Islands: Using a Hybrid Modelto Assess the Function of Central Macrophage
S. Fischer , † , P. Kurbatova , , † , N. Bessonov , O. Gandrillon , ,V. Volpert , , F. Crauste , , ∗ September 18, 2018 ∗ Corresponding authorTel: +33 472 448 516; Fax: +33 472 431 687; Email: [email protected] Universit´e de Lyon, CNRS, INRIA, INSA-Lyon, LIRIS Combining, UMR5205F-69621, France Universit´e de Lyon, Universit´e Lyon 1, CNRS UMR 5208, Institut Camille Jordan43 blvd du 11 novembre 1918, F-69622 Villeurbanne-Cedex, France INRIA Team Dracula, INRIA Grenoble Rhˆone-Alpes Center Institute of Mechanical Engineering Problems, 199178 Saint Petersburg, Russia Universit´e de Lyon, Universit´e Lyon 1, CNRS UMR 5534, Centre de G´en´etique et dePhysiologie Mol´eculaire et Cellulaire, F-69622 Villeurbanne-Cedex, France † Both authors contributed equally to this workShort title: Modelling Erythroblastic Islands
Funding Statement:
This work has been partially supported by the ANR grant ANR-09-JCJC-0100-01, project ProCell, and the ANR Project Bimod. The funders had no rolein study design, data collection and analysis, decision to publish, or preparation of themanuscript.
Author Contributions:
Conceived and designed the experiments: SF, PK, VV. Performed theexperiments: SF, PK. Analyzed the data: SF, PK, OG, VV, FC. Conceived and designedthe mathematical models: SF, PK, OG, VV, FC. Designed and implemented algorithms: SF,PK, NB. Wrote the paper: SF, PK, OG, VV, FC.1 bstract
The production and regulation of red blood cells, erythropoiesis, occurs in thebone marrow where erythroid cells proliferate and differentiate within particular struc-tures, called erythroblastic islands. A typical structure of these islands consists in amacrophage (white cell) surrounded by immature erythroid cells (progenitors), withmore mature cells on the periphery of the island, ready to leave the bone marrow andenter the bloodstream. A hybrid model, coupling a continuous model (ordinary dif-ferential equations) describing intracellular regulation through competition of two keyproteins, to a discrete spatial model describing cell-cell interactions, with growth factordiffusion in the medium described by a continuous model (partial differential equations),is proposed to investigate the role of the central macrophage in normal erythropoiesis.Intracellular competition of the two proteins leads the erythroid cell to either prolifer-ation, differentiation, or death by apoptosis. This approach allows considering spatialaspects of erythropoiesis, involved for instance in the occurrence of cellular interactionsor the access to external factors, as well as dynamics of intracellular and extracellularscales of this complex cellular process, accounting for stochasticity in cell cycle durationsand orientation of the mitotic spindle. The analysis of the model shows a strong effect ofthe central macrophage on the stability of an erythroblastic island, when assuming themacrophage releases pro-survival cytokines. Even though it is not clear whether or noterythroblastic island stability must be required, investigation of the model concludesthat stability improves responsiveness of the model, hence stressing out the potentialrelevance of the central macrophage in normal erythropoiesis. uthor Summary Red blood cells are produced in the bone marrow and released into the bloodstream. Dif-ferent cell populations involved in their production (from the very immature ones to dif-ferentiated ones) are usually kept in a steady state through various and complex negativefeedback controls. Inside the bone marrow, immature red blood cells mature and differen-tiate within particular structures, called erythroblastic islands, consisting of groups of cellssurrounding a macrophage (a big white cell). The role of the macrophage as well as thedifferent feedback mechanisms involved in red blood cell production (erythropoiesis) are notcompletely understood. We propose a new model of red blood cell production consideringboth mechanisms at the intracellular scale (protein regulation, cell fate decision) and at thecell population scale (cell-cell interactions, influence of feedback controls), and allowing thenumerical simulation of an erythroblastic island. Our results bring new information on therole of the macrophage and the erythroblastic island in the production and regulation of redblood cells, the macrophage being shown to play a key role in the stability and robustness ofthe production process, and will improve our understanding of in vivo erythropoiesis.
Hematopoiesis, the process of production and regulation of all blood cell types, has beenthe topic of modelling works for decades. Dynamics of hematopoietic stem cells have beendescribed by Mackey’s early works [34, 35], and later developed by Mackey and co-authors[36, 37] as well as Loeffler and co-workers [49, 63, 64]. The white blood cell production process(leukopoiesis) has been modelled in order to illustrate some pathological cases, such as effectsof radiations [63, 64], and bring some information on hematological diseases, mainly leukemia[10, 16, 24, 39, 40, 41, 57] and cyclical neutropenia [9, 17, 30]. The platelet production process(thrombopoiesis) attracted less attention through years [23, 62], with however some modellingeffort in the 2000’s [5, 52]. The red blood cell production process (erythropoiesis) has recentlybeen the focus of modelling in hematopoiesis. Early works by Wichmann et al. [65], Wulffet al. [68], B´elair et al. [8], Mahaffy et al. [38], on the regulation of cell population kinetics(mainly erythroid progenitors and mature red blood cells) have been completed by recentworks [1, 2, 6]. Particular attention was paid to the relevance of erythroid progenitor self-renewal in response to stress and considering multiscale approaches [11, 12, 18, 19], that isincluding both cell population kinetics and intracellular regulatory networks dynamics in themodels, in order to give insight in the mechanisms involved in erythropoiesis. A case studyof erythropoiesis is anemia [53], which allows observing the system in a stress situation thatcan be easily induced (at least in mice). A model of all hematopoietic cell lineages has beenproposed by Colijn and Mackey [16, 17], including dynamics of hematopoietic stem cells andwhite cell, red blood cell and platelet lineages.Almost all above-mentioned contributions to hematopoiesis modelling are based on theo-retical studies of deterministic models, and are mainly qualitative rather than quantitative.Other approaches, developed by Veng-Pedersen and co-workers [14, 25, 61, 66, 67], focused onpharmacodynamics/pharmacokinetics (PK/PD) models related to erythropoiesis modelling.Such approaches are usually centered on parameter estimation and model fitting to data (sta-tistical approaches). Nevertheless, none of the previously mentioned approaches did considerspatial aspects of hematopoiesis. Models describe cell population kinetics, either in the bonemarrow or the spleen (where cell production and maturation occurs) or in the bloodstream(where differentiated and mature cells ultimately end up). Consequently, cellular regulationby cell-cell interaction was neither considered in these models.We propose a new multiscale model for cell proliferation [11, 12], applied to erythropoiesisin the bone marrow, based on hybrid modelling. This approach, taking into account both in-teractions at the cell population level and regulation at the intracellular level, allows studying3ell proliferation at different scales. Moreover, the ‘hybrid’ modelling consists in consideringa continuous model at the intracellular scale (that is, deterministic or stochastic differentialequations), where protein competition occurs, and a discrete model to describe cell evolution(every cell is a single object evolving off-lattice and interacting with surrounding cells andthe medium), hence allowing considering small population effects as well as random effects.Hybrid multiscale models have been recently developed and used mainly, but not exclusively,to model solid tumor growth [4, 31, 32, 43, 45, 46, 51, 55].Erythropoiesis consists in the differentiation and maturation of very immature bloodcells, called progenitors, produced by differentiation of hematopoietic stem cells in the bonemarrow. These erythroid progenitors differentiate through successive divisions and under theaction of external signals into more and more mature cells, up to a stage called reticulocytes,which are almost differentiated red blood cells. After nuclei ejection, reticulocytes enterthe bloodstream and finish their differentiation process to become erythrocytes (mature redblood cells). At every step of this differentiation process, erythroid cells can die by apoptosis(programmed cell death), and erythroid progenitors have been shown to be able to self-renewunder stress conditions [7, 26, 28, 44]. Moreover, cell fate is partially controlled by externalfeedback controls. For instance, death by apoptosis has been proved to be mainly negativelyregulated by erythropoietin (Epo), a growth factor released by the kidneys when the organismlacks red blood cells [33]. Self-renewal is induced by glucocorticoids [7, 26, 28, 44], but alsoby Epo [50, 56] and some intracellular autocrine loops [28, 54].In the hybrid model, cells are individual objects supposed to be able to either self-renew,differentiate or die by apoptosis, according to the state of the environment and the cellitself. Indeed, global feedback control mediated by the population of mature red blood cells(through erythropoietin release) is supposed to influence erythroid progenitor proliferation[50, 56] and inhibit their death by apoptosis [33]. Local feedback control, mediated byreticulocytes and based upon Fas-ligand activity, is supposed to induce both differentiationand death by apoptosis [22]. Global and local feedback controls modify the activity ofintracellular proteins. A set of two proteins, previously considered by the authors [19], Erkand Fas, has been shown [50] to be involved in an antagonist loop where Erk, from the MAPKfamily, inhibits Fas (a TNF family member) and self-activates, high levels of Erk inducingcell proliferation depending on Epo levels, whereas Fas inhibits Erk and induces apoptosisand differentiation. Competition between these two proteins sets a relevant frame to observe,within a single cell, all three possible cell fates depending on the level of each protein.In addition to intracellular regulation of cell fate and feedback induced either by maturered blood cells or reticulocytes on erythroid progenitors, an important aspect of erythro-poiesis lays in the structure of the bone marrow. Erythroid progenitors are indeed producedin specific micro-environments, called niches [60], where erythroblastic islands develop. Anerythroblastic island consists in a central macrophage surrounded by erythroid progenitors[15]. The macrophage appears to act on surrounding cells, by affecting their proliferation anddifferentiation programs. For years, however, erythropoiesis has been studied under the in-fluence of erythropoietin, which may induce differentiation and proliferation in vitro withoutthe presence of the macrophage. Hence, the roles of the macrophage and the erythroblasticisland have been more or less neglected. Consequently, spatial aspects of erythropoiesis (andhematopoiesis, in general) have usually not been considered when modelling cell populationevolution or hematological diseases appearance and treatment.We use a hybrid model to take into account both intracellular and extracellular regulationof erythropoiesis, but also mainly to study the importance of spatial structure of erythroblas-tic islands in the regulation of erythropoiesis. To our knowledge, this is the first attempt tomodel erythropoiesis by taking these aspects into account. We focus in particular on the roleof the macrophage, by analyzing different situations, with and without macrophage at thecenter of the island, and actions of the feedback controls, especially regarding the feedback ofthe macrophage on erythroid cell proliferation and differentiation. The ‘classical’ structure4f the erythroblastic island, with the macrophage at the center and immature erythroid cellssurrounding it, will be shown to make the hybrid model stable (so the island can be sustainedfor an unlimited number of cell cycles) and robust to perturbations due to global and localfeedback controls.The next section is devoted to the presentation of the multiscale hybrid model. First,we focus on the intracellular scale, detailing regulatory networks considered in this work,described by nonlinear ordinary differential equations, and investigating their properties,in particular regarding the existence and stability of steady states associated to cell self-renewal, differentiation and death by apoptosis. Second, we concentrate on the extracellularscale, which consists in erythroid cell populations, possibly a macrophage, bone marrowmedium, and is described by a computational model. Interactions between cells and with thesurrounding medium are presented. Finally, we justify the coupling of both scales, throughglobal and local feedback controls, and we discuss parameters of the model. In Section 3we present the analysis of the hybrid model. This analysis consists of two parts: first, weconsider an erythroblastic island without a macrophage in its center, and investigate thestability of the island and the roles of feedback controls in the stability. Next we studyan island consisting of immature erythroid cells surrounding a macrophage, and the roles offeedback controls in its stability. This analysis concludes to the central role of the macrophagein the stability of the erythroblastic island, and consequently of the relevance of consideringspatial aspects when modeling erythropoiesis.
We introduce, in this section, the hybrid model for erythropoiesis used later for in silicoexperiments. It consists in the coupling of two models, with different space and time scales.The first model describes intracellular dynamics, represented by regulatory networks basedon specific protein competition. This intracellular dynamical level can be easily modeledwith a continuous approach, namely ordinary differential equations (ODEs) and/or partialdifferential equations (PDEs). The second model is at cell population level, where discretecells and events are computed, so that stochasticity due to random events (cell cycle duration,orientation of the mitotic spindle at division) and small population effects plays an importantrole. Extracellular regulation by growth factors is partially due to continuous models (PDEs).Such discrete-continuous models are usually called hybrid models [12, 31, 43, 45, 46, 51, 55].We hereafter present how the model is defined at each scale and how the different levelsinteract.
The intracellular scale describes how each erythroid progenitor cell chooses between self-renewal, differentiation and apoptosis. This choice depends on the level of expression ofsome proteins, involved in a regulatory network.Ordinary differential equations are used to describe a simplified regulatory network, basedon two competing proteins, Erk and Fas. Concentrations of the proteins evolve according toprotein-related mechanisms as well as extracellular factor concentrations. A brief analysis ofthe dynamics of the system (2.1)–(2.2) is performed in order to identify key factors in cellfate choice.
Precise intracellular regulatory mechanisms involved in erythroid progenitor cell fate arelargely unknown. Based on the current knowledge, we decided to focus on a simplifiedregulatory network based on two proteins, Erk and Fas, responsible respectively for cell5elf-renewal and proliferation and cell differentiation and apoptosis [18, 21, 22, 50]. Theseproteins are antagonists: they inhibit each other’s expression. They are also subject toexternal regulation, through feedback loops.One of these feedback loops is based upon the population of reticulocytes (differentiatederythroid cells). They produce Fas-ligand which is fixed to their exterior cell membrane.Fas-ligand activates Fas, a transmembrane protein, and influences progenitor differentiationand apoptosis. Another feedback control is related to erythrocytes in bloodstream. Theirquantity determines the release of erythropoietin and other hormones, called growth fac-tors. Erythropoietin is known to inhibit erythroid progenitor apoptosis [33] and to stimulateimmature erythroid progenitor self-renewal [50, 56]. Other hormones, like glucocorticoids[7, 26, 28, 44], also increase erythroid progenitor self-renewal by activating Erk.The simplified system of Erk-Fas interactions considered as the main regulatory networkfor erythroid progenitor fate is then [19] dEdt = ( α ( Epo, GF ) + βE k )(1 − E ) − aE − bEF, (2.1) dFdt = γ ( F L )(1 − F ) − cEF − dF, (2.2)where E and F denote intracellular normalized levels of Erk and Fas. Equation (2.1) describeshow Erk level evolves toward maximal value 1 by activation through hormones (function α of erythropoeitin, denoted by Epo , and other growth factors, denoted by GF ) and self-activation (parameters β and k ). In the meantime, Erk is linearly degraded with a rate a and is inhibited by Fas with a rate bF .Equation (2.2) is very similar, only there is no proof for Fas self-activation. Fas is howeveractivitated by Fas-ligand, denoted by F L , through the function γ ( F L ), it is degraded with arate d and inhibited by Erk with a rate cE .This model simply describes a competition between two proteins (Erk and Fas) throughcross-inhibition, self-activation (only for Erk), and activation by extracellular proteins.In order to explore the parameter space efficiently and observe the possible behaviors ofthe cell, it is relevant to study the mathematical interpretation of parameters. For fixed values of
Epo , GF and F L , (2.1)-(2.2) is a closed system of ordinary differentialequations. In order to analyze the impact of each parameter on the system’s behavior, wecan derive analytical forms of the zero-lines associated with each equation. The zero-line f associated with equation (2.1) is f ( E ) = 1 b (cid:16) αE − βE k − (1 − E ) − ( a + α ) (cid:17) . (2.3)The zero-line g associated with equation (2.2) is g ( E ) = γcγ + dc + E . (2.4)The shape of the zero-lines depends only on 4 and 2 parameters respectively ( α/b , β/b , k and a/b for f , γ/c and d/c for g ), the 2 remaining parameters ( b and c , related to cross-inhibition) only affecting the absolute and relative speeds of the dynamics along each dimen-sion. The existence of steady states, defined as intersections of the two zero-lines, thus onlydepends on the former 6 parameters. 6et us stress out that α represents the influence of external factors (erythropoietin, growthfactors) on Erk activation, β is related to Erk self-activation, and γ describes the influenceof Fas-ligand on Fas activation.By considering that f can be separated in 3 components, the impact of each parametercan easily be observed. The same can be done for g . Figure 1 shows that the variations ofthe zero-lines are quite constrained, with a central role of α and β for f and γ for g . Thiswill be studied in further detail in the following.One may note that the effect of k is not shown here, since k = 2 in all experiments. It ishowever noticeable that increasing k simply decreases the maximum value of the polynomialpart of f toward 0 while increasing the value for which it is reached toward 1.Simple considerations show that the zero-lines can only cross in the [0 , × [0 ,
1] domain.Furthermore, they intersect at least once, so there is at least one steady state, and there canbe only up to three steady states. When there is only one steady state, it is stable. Whenthere are three steady states, the central point is unstable and the two surrounding pointsare stable [19]. The case with two steady states is singular so it is not really relevant.Hence, the system exists in two configurations that we will call monostable (one stablestate) and bistable (two stable states). The next step is to determine how the system cango from one state to the other one. Thanks to Figure 2, two kind of bifurcations can beimagined to achieve this goal.Let first start with the case where β is large compared to α , this means Erk self-activationrate is stronger than external activators (Epo for instance). In this case, f has two localextrema (Figure 2.A). If γ , which is associated with Fas-ligand, varies then a double saddle-point bifurcation (or fold bifurcation) occurs. For low and high values of Fas-ligand (thatis low and high values of γ , respectively), there is only one stable steady state, whereas formedium values of Fas ligand two stable points coexist, a separatrix determining the attractivedomain of each point (Figure 2.B).Now consider the case where α is large compared to β . Then f becomes strictly decreasingand, if α is sufficiently large, there will be only one steady state for every value of γ . Thesystem is thus strictly monostable (Figure 2.C).Finally, the system (2.3)–(2.4) exists in two different cases: when α is low it is locallybistable, depending on the value of γ , and when α is high, it is strictly monostable. A kind offork bifurcation occurs when α increases (approximately around α = β/
27, but exact valuedepends on other parameters and is quite difficult to obtain).
Let us focus on the possible dynamics of the system (2.1)–(2.2). This system must be able todivide the erythroid progenitor cell population into 3 different types: self-renewing, apoptoticand differentiating cells. Differentiation corresponds to low values of Fas and Erk at the endof the cycle, apoptosis to high values of Fas and self-renewal to high values of Erk andreasonable values of Fas.A simple idea to achieve such a behavior is to introduce two threshold values. Thefirst one, F cr , induces apoptosis when reached. The second one, E cr , determines the choicebetween self-renewal and differentiation for cells reaching the end of their cycle. So the phaseplane of system (2.1)–(2.2) can be divided in three parts: cells having low Erk and Fas valuesdifferentiate at the end of their cycle, cells having Erk values greater than E cr self-renewat the end of their cycle and cells reaching Fas values greater than F cr die immediately byapoptosis.In order to get three robust cell subpopulations, it is necessary to determine how to cor-rectly place these thresholds. Therefore, cases where the system is either strictly monostableor locally bistable have to be considered.Let us recall that the position of the f zero-line strongly depends on parameters α and7 , whereas the position of the g zero-line strongly depends on the parameter γ . Since β denotes the self-activation rate of Erk and α represents a feedback control by global factors(this will be developed later in Section 3), then all erythroid progenitor cells will obey thesame evolution rule for Erk. On the contrary, γ represents a feedback control by a localfactor (Fas-ligand emitted by reticulocytes and acting at short-range, see Section 2.2) so allcell positions on the phase-plane will depend on their own exposition to Fas-ligand.If the system is strictly monostable, then, independently from its starting position in the( E, F )-plane, a cell will converge toward a unique stable state, located on the f zero-lineand defined only by its exposition to Fas-ligand (the value of γ ). So, at steady-state, if cellsare exposed to a continuum of Fas-ligand values, they will be continuously located alongthe f zero-line, and their level of Erk will not change: there is no real subpopulation. It isthus quite inefficient to work at steady-state (lack of robustness) to get three subpopulationscorresponding to the three possible cell fates.On the contrary, if the system is locally bistable, cells will still be asymptotically locatedalong the f zero-line, yet a part of it is unstable (roughly speaking the ascending part, seeFigure 1) so cells will be separated into two subpopulations. In this case, the starting pointis important: depending on where it is compared to the separatrix in the bistable domain, acell will converge either toward the upper or the lower stable branch of f .Based on these observations, we considered two different initial conditions and dynamicsto get the three subpopulations.The first idea is that, if working at steady-state does not allow to obtain distinct subpop-ulations, then working out of steady-state may be relevant. Consequently, in a first scenariocalled the ‘out of equilibrium case’, all cells are arbitrarily put at the origin (that is E = 0and F = 0) at the beginning of their cycle (see Figure 5, left) and, before they have reachedsteady-state, they are forced to take a decision based on Erk and Fas levels.In the case where the system is strictly monostable, this solution is nevertheless not sat-isfactory (no subpopulation appears because cells remain very close to each other). However,when the system is locally bistable, cells will converge either toward the lower stable branchor the upper stable branch of the f zero-line. For a given value of Fas-ligand, the separatrixwill cut the origin, so there is a clear threshold between the two states. Even if there is stilla continuum of cells on the phase plane for continuous values of exposition to Fas-ligand,two denser groups appear. In fact, cells exposed to high Fas-ligand values quickly convergetoward the top, cells exposed to low Fas-ligand values toward the bottom-right and thoseexposed to medium values remain trapped near the origin in a ‘slow’ zone around the sepa-ratrix. The latter cells keep low Erk and Fas values, so by introducing the thresholds, threesubpopulations are easily obtained (Figure 3).The second idea is that working at steady-state is more robust. When the system is lo-cally bistable, one notices that cells on the lower stable branch are those which will self-renew(high values of Erk, reasonable values of Fas). If quantities of Erk and Fas are equally dividedbetween the two daughter cells at birth, those exposed to medium values of Fas-ligand will ap-pear near the separatrix and be divided between the upper branch (differentiation/apoptosis)and the lower branch (self-renewal). We thus have a system where cells rapidly take the de-cision between self-renewal or differentiation/apoptosis. Once the decision is made, cells aretrapped on a hysteresis cycle controlled by Fas-ligand, so the decision is hardly reversible.Fas values have to become very high in order to get from the lower branch to the upperbranch and conversely, very low in the opposite direction (Figure 4).For this solution, called ‘hysteresis cycle case’, introducing thresholds in order to deter-mine cell fate is quite simple. The Erk threshold can be situated anywhere between the twobranches and the Fas threshold somewhere on the upper branch, cutting into two the cellson this branch and so easily controlling the differentiation/apoptosis ratio.The ‘out of equilibrium’ and the ‘hysteresis cycle’ scenarii are schematically described inFigure 5. 8 .1.4 Feedback Control Role A key point in system (2.1)–(2.2) is how it reacts to feedback controls. The analysis performedin the previous subsections dealt with the influence of each parameter on the system dynamics,mainly parameters related to feedback controls ( α and γ ). Let now study how parameters areaffected by feedback controls. Erythropoietin (through function α ) is expected to decreaseapoptosis and enhance self-renewal/differentiation. We also know that, at a smaller spacescale, Fas-ligand induces apoptosis (through function γ ) and other growth factors, on theopposite, stimulate self-renewal (through the function α ).In order to validate possible dynamics, the effect of each factor must be checked in bothcases, the ‘out of equilibrium’ and the ‘hysteresis cycle’ scenarii.The control exerted by Fas-ligand, which increases the value of Fas reached at steady-state, has already been largely discussed. Simply by existence of the threshold F cr , Fas-ligandplays the expected role. Fas-ligand is linked to the parameter γ , which is actually a functionof F L .Epo is known to have a direct action on apoptosis [33], independent of Erk pathway[13, 42, 59]. Therefore, a simple idea is to increase F cr with the level of Epo. This clearlyreduces apoptosis in favor of differentiation and self-renewal. However, it is also knownthat Epo, and other growth factors, activate Erk. So erythropoietin and growth factors aresupposed to influence the variations of α . It is noticeable that at a local scale, many growthfactors are supposed to have similar mechanisms, so α should in fact be a function of manyglobal and local growth factors.The value of α controls the monostable/bistable bifurcation (Figure 1). When the systemis locally bistable, α mainly controls the value of the local minimum of the f zero-line. In thecase where we work out of equilibrium, α changes the position of the separatrix and of theupper steady states, so when α is increased the separatrix cuts the origin for a higher value of γ ( F L ) and cells tend to be displaced toward steady-states with lower values of Fas and highervalues of Erk. In the case where the hysteresis cycle is used, increasing the local minimumof f reduces the size of the cycle. Cells located on the upper stable branch instantly ‘fall’ onthe lower stable branch, so self-renewal is immediately enhanced and differentiation reduced.The effect of α is not the same in both dynamics, but finally it is clearly possible todecrease apoptosis in favor of self-renewal/differentiation and the trends for each protein isroughly what is expected biologically. However, it is still necessary to study whether thefeedback is as efficient as observed experimentally, when the intracellular scale is coupled tothe extracellular one (see Section 3). In the previous section, we focused on the intracellular scale of the hybrid model, based ona system of ordinary differential equations describing competition between two key proteins,Erk and Fas. Concentrations of Erk and Fas induce either cell self-renewal, or differentiation,or death by apoptosis. We now focus on an other layer of the hybrid model, the discretepart, which describes how cells interact and how they influence intracellular parameters whichdepend on extracellular molecules.In order to describe the evolution of an erythroid cell population, we chose to use anindividual-based model. All cells and their interactions are then numerically computed.The objective being to represent erythroblastic islands, which have a limited size, such anapproach takes automatically into account small population effects as well as random effects(direction of division, cell cycle length).A population of cells is numerically simulated in a 2D computational domain which isa rectangle. Each cell is a discrete object, an elastic ball, considered to be circular andcomposed of two parts: a compressible part at the border and a hardly compressible part atthe center. All newborn cells have the same radius r and linearly increase in size until the9nd of their cycle, when they reach twice the initial radius. When a cell divides, it gives birthto two small cells side by side, the direction of division being randomly chosen. The cell cycleduration is itself variable. From a biological point of view, cell cycle proceeds through G /G , S , G and M phases. Duration of the first G /G phase and transitions between phases arecontrolled by various intracellular and extracellular mechanisms, inducing stochasticity incell cycle durations [29, 47]. We assumed the duration of G /G phase is a random variablewith a uniform distribution in some given interval, other phase durations were supposed tobe constant.Under the assumption of small deformations, the force acting between cells can be ex-pressed as a function of the distance between their centers. The force between two particleswith centers at x i and x j is given by a function D ( d ij ) of the distance d ij between the centers.This function is zero if the distance is greater than the sum of their radii. To describe themotion of a particle, one should determine the forces acting on it from all other particlesand possibly from the surrounding medium. The motion of each cell is then described by thedisplacement of its center by Newton’s second law: m ¨ x + µm ˙ x − X j = i D ( d ij ) = 0 , (2.5)where m is the mass of the particle, the second term in left-hand side describes the frictionby the surrounding medium, the third term is the potential force between cells. The forcebetween two spherical particles is considered in the form D ( d ij ) = K d − d ij d ij − d + 2 H , d ij < d , , d ij ≥ d , (2.6)where d is the sum of cell radii, K is a positive parameter, and H accounts for the compress-ible part of each cell. The force between the particles tends to infinity when d ij decreases to d − H . On the other hand, this force equals zero if d ij ≥ d .Three cell types are computed. First, erythroblasts, which are immature erythroid cells,also known as erythroid progenitors. They follow the growth rules explained above and theirfate is determined as described in Section 2.1. They either self-renew and give two cellsof the same type, or differentiate and give two differentiated cells (reticulocytes), or die byapoptosis, depending on their exposition to growth factors and Fas-ligand. It is noticeablethat no asymmetric division is considered in this model.Second, reticulocytes. From a biological point of view, reticulocytes are almost mature redblood cells that leave the bone marrow and enter the bloodstream after ejecting their nuclei.In this individual-based model, they are differentiated cells which stay in the bone marrow alittle while after being produced, and leave the bone marrow (computational domain) after atime equal to one cell cycle duration. Contrary to erythroblasts, reticulocytes do not have achoice to make, they only express Fas-ligand, thus influencing the development of surroundingerythroblasts.Third, macrophages. They are big white blood cells located at the center of erythroblasticislands. They very probably have an active role in the development of blood cells surround-ing them. In this model, they express growth factors, driving nearby erythroblasts towarddifferentiation and self-renewal.Macrophages express growth factors that are supposed to diffuse in their neighborhood.Reticulocytes express Fas-ligand on their surfaces which does not diffuse in a real bone mar-row, yet the expression of Fas-ligand is modeled by short-diffusion. This allows consideringseveral phenomena that are not taken into account in the hybrid model: first, cell motion inthe bone marrow is actually more important than computed and, second, there are normallyseveral maturing erythroid subpopulations which increasingly express Fas-ligand, creating a10radient of exposition. Moreover, cells are not circular in reality and this short-range diffu-sion compensates their inexact representation. Both Fas-ligand ( F L ) and growth factor ( GF )concentration evolutions are described with the following reaction-diffusion equations, ∂F L ∂t = D F L ∆ F L + W F L − σ F L F L , (2.7) ∂GF∂t = D GF ∆ GF + W GF − σ GF GF, (2.8)where W F L = k F L C ret is a source term depending on the number of reticulocytes ( C ret denotes the relative volume of reticulocytes in the computational domain) and W GF is aconstant source term for growth factor (released by the macrophage) concentration, σ F L and σ GF are degradation rates, and D F L and D GF are diffusion rates. If the diffusion coefficient D F L is sufficiently small, then Fas-ligand is concentrated in a small vicinity of reticulocytes.In this case, Fas-ligand influences erythroid progenitors when they are sufficiently close toreticulocytes. This is the short-diffusion illustrated in this work.Figure 6 shows an example of erythroblastic island in the hybrid model. The big greencell in the center is a macrophage, the green substance surrounding it is the above men-tioned growth factors released by the macrophage. At the border, blue cells are reticulocytesemitting Fas-ligand (whose concentration is represented in red). Between them are the ery-throblasts, growing in size until they die or divide into two new cells, depending on how theywere influenced by the two diffusing proteins.A fourth cell type, erythrocytes, is considered in this model. Biologically, erythrocytesare mature reticulocytes, circulating in blood and transporting oxygen to the organs. Inthe model, erythrocytes are reticulocytes that have spent a time equal to one cell cycle inthe computational domain and then left it. Hence, erythrocytes are only counted as cellsoutside the domain (they are not represented on the computational domain) which act,through feedback loops, on the fate of immature cells within the domain. We will assimilateerythrocytes and red blood cells (RBC) in the following, by considering that RBCs are onlyerythrocytes (in reality, some reticulocytes usually manage to leave the bone marrow andenter the bloodstream where they end their maturation, so they are also red blood cells).Erythrocytes will be supposed to survive (arbitrarily) during a time corresponding to 3 cellcycle durations (although RBCs do not cycle). This is short, since in mice RBCs have alifespan of 40 days, yet it allows performing a first feedback test and it can be changed whencomparing the model to experimental data. Now that both modeling scales have been described and a general picture of the model isaccessible, it becomes necessary to detail how the two scales can be linked together. At thehighest level, there is a small erythroblastic island formed with two populations of devel-opping red blood cells and, optionally, a central macrophage. The island has a particulartopology: the macrophage is at the center, differentiated cells (reticulocytes) at the borderand immature cells in between. Both macrophage and differentiated cells emit growth factorscontinuously controlling intracellular protein concentrations: immature cells at the center arerather exposed to growth factors produced by the macrophage whereas those at the borderare rather exposed to Fas-ligand. In addition, immature cells are subject to a feedbackcontrol mediated by mature red blood cells circulating in the bloodstream, representing theaction of erythropoietin. Concentration of Epo in the computational domain is supposed tobe uniform, so all cells are similarly influenced by Epo. This assumption holds true for allexternal substances acting on cells at a global scale, yet other more sophisticated choices arepossible.Erythroid progenitors are then exposed to a continuum of erythropoietin
Epo , growthfactors GF and Fas-ligand F L . Regarding the intracellular scale, this means that α and γ α drives erythroblaststo self-renewal, high γ to death, and intermediate values to differentiation. Hence, thecell population scale influences the intracellular scale through Epo, growth factors and Fas-ligand concentrations, that is through functions α and γ and the critical value F cr of Fasconcentration (Epo increases F cr in order to inhibit cell apoptosis).One can notice that when there is no macrophage, there is only one value of α for allerythroblasts (this value depends on the concentration of Epo which is related to the numberof erythrocytes, and is supposed to be the same for all immature cells), so cells are simplysituated along a unique f zero-line. It is then convenient to investigate how the value of γ will determine a cell’s choice (see Section 2.1.4).As just mentioned, Epo-mediated feedback control simply corrects the critical value forapoptosis F cr and values of α uniformly, so its impact can be easily understood individually byreferring to the analysis in Section 2.1.4. However, in order to compute this feedback control,the number of erythrocytes in blood circulation must be estimated. This is performed fromthe reticulocyte population as mentioned at the end of the previous section. In the previous section, the hybrid model has been presented in details. The applicability ofsuch a modeling will be investigated in this part. To do so, we focus on finding conditions(parameter values, that are partially discussed in Appendix A) for the system leading to stableerythroblastic islands (this means the number of cells in the island stays almost constant fora certain time when not facing a perturbation), reacting sensibly to global feedback variables(for example in the case of anemia) and proving to be robust to parameter values variations.We begin by investigating the case of one erythroblastic island without macrophage in itscenter. The case of the island with a macrophage is then considered to see how the stabilityof the island is affected. For these analyses, simple numerical stability and feedback testswere performed to check the validity of hypotheses.
Consider an island initially composed of 98 erythroid progenitors and 64 reticulocytes sur-rounding them. Parameters are carefully chosen to optimize stability and response to feed-back (see Appendix A).
Erythroblastic islands are likely to be stable. Indeed, one can expect that, in order to bebiologically responsive, an island survives for several cell cycles. It is relevant to determinewhether it is mathematically possible to achieve long-lasting stability with the hybrid modelor not.A wide range of parameters has been tested for the two possible dynamics, that is eitherout of equilibrium or with the hysteresis cycle. At first, the best results look like quite similar(Figure 7). In both cases, the island remains approximately stable during some cycles but,due to stochastic variations (cell cycle duration, orientation of the mitotic spindle at division,small size populations), it suddenly dies out or excessively proliferates (it must be noticedthat saturation is only due to space limitations). Two questions arise from this fact: whydoes the system behave this way and is it still biologically relevant?Let’s focus first on the mechanisms leading to an island’s stability. In the simulations,during each cycle there is a turnover of reticulocytes, in order to replace those which entered12loodstream. These reticulocytes then induce apoptosis and differentiation in the erythrob-last population. Within this population, in order to remain stable, the proportion of self-renewing erythroblasts has to be constant. From a geometrical point of view, self-renewingerythroblasts are located at the center of the erythroblastic island, differentiating and dyingerythroblasts at the border (Figure 6). When a random variation occurs in the size of theisland, the amount of the differentiating and dying cells varies like the perimeter of the circle(occupied by reticulocytes which act at constant range), while the amount of self-renewingcells varies like the area of the circle. As a result, the ratio ‘differentiating and dying cells’over ‘self-renewing cells’ does not remain constant. When it derives too far from the initialvalue, one subpopulation rapidly overgrows the other.It appears then that it is impossible to find parameter values for which an island wouldremain at steady-state on a long term. However, this does not imply that results are bi-ologically irrelevant. The fact that islands die in this model is not really surprising since,biologically, there might be a turnover of these structures and erythroid progenitors perma-nently arise from the differentiation of hematopoietic stem cells, which has been neglectedhere. Also, proliferation is not as dramatic as it might seem. First, it is possible to turn itdown by varying parameters (but sacrificing instead the average lifetime of an island) andalso proliferation may not be that simple in a realistic environment. Indeed, there is noobstacle in the computational domain and some particular geometries could clearly reducevariations in the proportion of self-renewing cells. Another point is that when several islandsare side-by-side, proliferating islands will collide and be in contact with more reticulocytesthan before, which may block their growth (data not shown). Nevertheless, unstable islandsrevealed hard to control and stability of islands should be expected, as previously mentioned.Regarding the two possible dynamics, either out of equilibrium or hysteresis cycle, stabil-ity of the system at equilibrium (hysteresis cycle) is slightly better (Figure 7 right). Moreover,a significant fact arose during experiments when distribution of cell cycle lengths was varied:the dynamics at equilibrium is very robust in this case, whereas the other dynamics becomesclearly biased. Indeed, when cells are placed at the origin at the beginning of their cycle,they start in the differentiating area. When their cell cycle length is short, they do not haveenough time to escape this area and differentiate. When cell cycle duration is long, they can-not differentiate. Consequently, some cell fates are only decided because of cell cycle length,independently of the position in the erythroblastic island, and the island structure is quicklylost. In conclusion, since cell cycle length variation may be large [29, 47], it was decided toabandon this ‘out of equilibrium’ dynamics and to perform further tests with the dynamicsusing the hysteresis cycle (Figure 5, right panel).
The above-described hybrid system should be clearly influenced by feedback control loops,as expected biologically, and, besides, feedback loops may play a role in the stabilization oferythroblastic islands.As mentioned in Section 2.1, different feedback controls are considered in this model: alocal control through Fas-ligand, activating Fas and influencing erythroblast differentiationand death by apoptosis; two global feedback controls by erythropoietin, one activating Erkthrough the function α , the second one varying the critical value of Fas, namely F cr , inducingapoptosis. We focus here on the global feedback control mediated by Epo.It is known that increasing Epo levels in the bloodstream, and consequently in the bonemarrow, increases Erk activity and decreases erythroblast apoptosis rate. Hence, an increaseof Epo levels induces an increase in the values of both α and F cr . The exact influence of Epoon these two quantities is however not known. Therefore, we began by simply monitoringthe behaviour of one island under constant values of these two parameters around the defaultvalue (Figure 8). 13y comparing the proportions of different erythroblast subpopulations (self-renewing, dif-ferentiating and apoptotic populations) during the first 5 cycles, it appeared that an increaseof F cr (that is, an increase of Epo levels) decreases apoptosis and self-renewal in favor ofdifferentiation. For low values of F cr , there is no differentiation (erythroblasts at the centerof the island survive, others die), so the island proliferates by absence of death factors. Forhigher values, effect on self-renewal can be tuned by modifying parameters controlling thehysteresis cycle. At some point the self-renewing population is not affected anymore and thereis simply a transfer from the apoptotic subpopulation to the differentiating subpopulation(the differentiation part of the differentiation/apoptosis branch increases).Increasing α , on the contrary, mainly acts on the self-renewing population. An in-crease of α sends all cells on the self-renewing branch of the hysteresis cycle and differen-tiation/apoptosis disappears. Reducing α increases the size of the apoptosis/differentiationbranch. For a constant value of F cr , only the differentiation part of this branch is increased.One thus observes a transfer from the self-renewing branch to the differentiation part of theother branch.In conclusion, the system’s response is qualitatively what is expected. Epo mainly de-creases apoptosis in favor of differentiation and stimulates self-renewal. However, the statis-tics computed above only indicate trends and it is important to see how islands actually reactto this feedback.The next step is then to explicitly include feedback loops in the system’s equations to seewhether the response is quantitatively correct or not, and if it is possible, through feedbackcontrols, to achieve a better stability than before. Therefore, we assumed a linear variationof both α and F cr with deficiency of red blood cells (RBC), hence implicitly describing thedependency on Epo. Denoting by N RBC the number of circulating red blood cells, and sincea decrease of RBC count induces an increase in Epo levels, one gets α := α ( N RBC ) = α + k α ( N target − N RBC ) , (3.9) F cr := F cr ( N RBC ) = F + k F ( N target − N RBC ) . (3.10)We remind that N RBC is estimated by the number of reticulocytes which have left bonemarrow and are supposed to survive during a time equivalent to 3 cell cycles. N target is thenumber of RBCs in circulation during a typical run before proliferation/extinction of theisland. Parameters α and F are ‘typical’ values of α and F cr (as given in Appendix A).Parameters k α and k F are to be estimated.It must be mentioned that these two feedback controls are not exact from a biologicalpoint of view, since one would expect the functions to saturate for low and high values of N RBC , yet they are good approximations when the number of red blood cells is not too faraway from the target (for instance, in normal erythropoiesis).Functions α and F cr in (3.9) and (3.10) were used together with the previous initialconditions (hysteresis cycle case), and it was impossible to achieve better stability thanpreviously mentioned when varying the values of k α and k F (data not shown). To be moreprecise, it was possible with strong feedback controls to stop proliferation, instead the islandsystematically died out. Other tests, based on simulations of a constant number of red bloodcells ( N RBC ), showed that the effect of feedback controls was not fast enough to deal withthe quick proliferation and extinction of erythroblastic islands. The island still lived around10–15 cycles and, when it began to expand or die out, feedback response came too late or wastoo strong to bring it back to its ideal size. Therefore, what was observed was only the usualdestabilization of the island, feedback being only able to turn proliferation into extinctionafter some cycles.As a result, it is very difficult to test feedback reliably on unstable islands. However, thestudy showed that global feedback loops are inefficient when it comes to controlling localstructures, since proliferation or extinction events occur quickly in the model.14n the next section, we consider the addition of a macrophage in the center of the islandon the stability of the island and the relevance of feedback controls.
In this case, an erythroblastic island is initially composed of one macrophage, 80 erythroidprogenitors and 64 reticulocytes (see Figure 6, right). As done in the previous section,stability of the island and influences of feedback controls are investigated.In addition to feedback controls exerted, in the previous section, on erythroblast dynamicsby Epo and Fas-ligand, the macrophage in the center of the island is supposed to releasegrowth factors, whose concentration is denoted by GF , which positively act on Erk activation.A number of proteins associated with the proliferation of erythroid progenitors and that couldfulfill such a function have been described. This includes Stem Cell factor (SCF), the c-kitligand, Ephrin-2, the EphB4 ligand, and BMP4, a member of the TGF β family ([48] andreferences therein). All those growth factors are assumed to diffuse around the macrophageas developed in Section 2.2 (Equation (2.8)). They trigger erythroblast self-renewal throughthe function α : the more growth factors, the higher α . Function α now becomes α := α ( N RBC , GF ) = α + k α ( N target − N RBC ) + k GF GF. (3.11)Meaning of parameters α , k α and N target is identical to the ones defined for (3.9). Valuesare also the same, except for α that had to be decreased to compensate the addition ofthe term k GF GF . Parameter k GF was set accordingly (influence of these parameters arediscussed below). Default values are α = 0 h − and k GF = 3 h − . molecule − . Stability of the hybrid system is first investigated without considering global feedback controls(that is, k α = 0 in (3.11) and k F = 0 in (3.10)). Results are illustrated in Figure 9.For a large range of parameter values, the system quickly converges toward a steady-state. The number of self-renewing and differentiating cells (and thus circulating RBC) canbe precisely controled by parameter values, as will be detailed in the next subsection. Thereis a first phase where the size of subpopulations oscillates (because of lack or abundance ofdeath factors) yet steady-state is always quickly reached (within a few cycles), even thoughinherent stochastic oscillations are always present.The macrophage completely controls the island. If α is sufficiently low, a cell will notbe able to self-renew without external growth factors. Therefore, the default behaviour of anisolated cell is differentiation. On the other hand, a cell in contact with the macrophage willbe exposed to a very high concentration of growth factors and will always self-renew, sinceit is generally not exposed to death factors. In the vicinity of the macrophage, there is acompetition between growth factors and death factors, which will eventually determine thesize of the island. If α or k GF gets higher, growth factors will reach further erythroblasts,thus the size of the final island will increase. Now that erythroblastic islands are stable, due to the presence of the macrophage, it is notnecessary to be limited to a period of 5 cycles to compare the subpopulations of the island,as was done in the previous case (Section 3.1.2). Rather, it becomes relevant to see how thesize of the island and the production of RBC could evolve on a long term.Similarly to what has been done in the previous case, we once again simulated the expectedconsequences (in terms of steady state value) of feedback control by running simulations underconstant values of α and F cr . Results are shown in Figure 10.15hen F cr increases, both the number of reticulocytes and RBC increases, as the differen-tiation part of the differentiation/apoptosis branch increases. For high values of F cr , there isno apoptosis anymore. The number of erythroid progenitors seems to be rather independentof the value of F cr .When α increases, which can be seen as an increase of the ground level of global ‘survivingfactors’ (here Epo and macrophage-emitted growth factors), the size of the island increasesmainly due to the increase of the self-renewing population. As a result, reticulocytes, situatedall around the island, also increase in number, as does the amount of RBC. However, for veryhigh values of α , surviving signals are so strong that almost all erythroblasts self-renew(the macrophage has no effect anymore) and the island is made only of proliferating cells,reticulocyte and RBC counts vanish. Such high values of α may only be reached in extremestress situations and mainly create a pool of cells which may afterwards differentiate andsubstantially increase the number of RBC at once (but this implies that another mechanismhas to lower the value of α at some point). In any case, these values of α may else be avoidedby adding a saturation effect on the feedback loop.In each case, a substantial increase of RBC production (about 3-fold and 4-fold respec-tively) is observed. As the two mechanisms are largely independent, they can be combinedand thus the overall production rate can be greatly increased, in magnitudes comparableto what can be observed biologically [15, 48, 58]. This seems to indicate that the reactionof the system with macrophage to feedback controls is more efficient, relevant and easy toobserve compared to the case without macrophage, again in good agreement with biologicalevidences [48]. We proposed a new model for erythropoiesis modeling based on the coupling of two relevantscales, the intracellular scale consisting of protein regulatory networks and the extracellularscale focusing on cell movement and interactions as well as growth factor distribution in themedium. Description of the competition between two proteins within erythroid progenitorshas been performed with continuous models (ordinary differential equations) whereas cellshave been studied as discrete objects on an off-lattice model, so the whole system can benamed ‘hybrid model’, for it consists in the coupling of continuous and discrete systems.Applied to normal erythropoiesis in the bone marrow, in particular to the regulation oferythroblastic islands, the model suggests an important role of macrophages in the stability ofislands. In the absence of macrophages, erythroblastic islands very quickly lose their stability,meaning they either die out or abnormally proliferate (with overproduction of erythroidprogenitors). They survive only for the equivalent of about ten cell cycles. On the contrary,when the erythroblastic island is built around a central macrophage, assumed to releasegrowth factors sustaining erythroid progenitor survival and proliferation (like SCF, Ephrin-2 or BMP-4 [48]), then the island is a very robust structure, able to produce continuouslyerythrocytes while keeping a steady structure (that is an almost constant ratio of reticulocytesover progenitors).This control of the production of differentiated red blood cells by the macrophage at thelevel of the erythroblastic island may first appear in contradiction with usual knowledge oferythropoiesis. It is indeed well known that during a response to a stress, erythropoiesis isa very intense process, large amounts of cells being produced in a short time, so the wholeprocess is expected to possess the ability to overcome its usual production, in particularthe ratio of reticulocytes over progenitors should not stay constant. The hypothesis wepropose is that the macrophage, inside the erythroblastic island, controls the explosiveness oferythropoiesis and that, contrary to other hematopoietic lineages (white cells, platelets), thiscontrol is necessary because stress erythropoiesis must deal with very large amounts of cells16nd otherwise the stability of the process could be lost. It is moreover noticeable that gettingmore stability with the central macrophage does not decrease responsiveness of the model.For instance, the island reacts to perturbations simulating anemia (sudden loss of maturered blood cells), returning quickly to ‘normal’ values of the different erythroid populations,progenitors, reticulocytes and erythrocytes (not shown here). Hence the gain of stabilityis not compensated by a loss in the system’s reactivity. Simulation of stress erythropoiesis(anemia) is actually a work in progress which should confirm that the macrophage plays arelevant role in the stability of erythropoiesis.Some comments may however be made on the behavior of an island without macrophagein its center. Although it is clear that this system cannot lead to a stable island (this can beexplained by considering that fluctuations automatically lead to proliferation or extinction),it is nevertheless possible to imagine a system working with such islands. For instance, ifislands are put side-by-side, proliferation is made difficult (islands cannot expand due toconfinement) and extinction can be compensated by the birth of new islands. Moreover,other assumptions, based on a better knowledge of erythroblastic islands, could lead to morestability (geometry of bone marrow and islands,. . . ), even though all scenarii cannot andwere not tested. However, the lack of stability makes it hard to have a clear view on howfeedback controls act on one single island, and as far as we know of, studying erythropoiesismodels (or hematopoiesis models, in a wide sense) without considering feedback controls doesnot make sense. Hence, searching for island stability appeared relevant in this study. At aglobal scale, it is still possible that island stability is not necessary and can be compensatedbetween proliferation and extinction. When looking at one single island however, the factthat feedback controls do not work correctly on an unstable island is somewhat compelling.Stability of the erythroblastic island including a central macrophage was obtained byconsidering two scenarii (see Figure 5) for the distribution of Erk and Fas quantities indaughter cells, following mitosis and mother cell’s division. In the first one, initial values ofErk and Fas in newborn cells were set to zero, and cell fate was decided depending on theevolution of concentrations through the cell cycle duration (this is the ‘out of equilibrium’assumption). This scenario was shown to lead to results strongly dependent on the length ofcell cycle. The second scenario, called the ‘hysteresis cycle case’, which is more biologicallyrelevant, consisted in dividing values of Erk and Fas of the mother cell at mitosis into two, soboth daughter cells are located on a hysteresis cycle in the Erk-Fas plane at birth. Combiningthe hysteresis cycle scenario with a structure including a local feedback on the island (themacrophage) gave relevant results and appears very promising. Parameter values which wereobtained through mathematical analysis appear to be very robust and the system exhibitsqualitatively the same behavior for a wide range of parameters. The stability of the islandallows a better control of its size and its subpopulations and, as a result, gives much moresatisfying results on feedback tests.Regarding parameter values used to study the model and check different scenarii (with orwithout macrophage in the center of the island, distribution of initial quantities of Erk andFas), one must note that most of these values are mainly not accessible through experiments.Some others, such as for instance the time of survival of RBC, are however rather quitearbitrary and not necessarily in agreement with usual knowledge on erythropoiesis. Yet,results show that the behaviour may be satisfying from a biological point of view. Sincea wide range of parameters is allowed, it is expected that this system can be tuned to fitexperimental data and numerically reproduce certain stress behavior, which will be the nextstep for modeling erythropoiesis.Feedback functions considered in this study may also be questioned. Indeed, both erythro-poietin’s influence on erythroid progenitor self-renewal capacity and on apoptosis protectionwere taken as linear functions of the number of circulating red blood cells. Saturation effectsshould be envisaged when confronting the model to strong perturbations (for instance, a se-vere anemia), in order to be closer to the reality where responses of the organism to a loss of17ed blood cells have shown to be nonlinear [18, 19]. In order to perform tests presented in thecurrent study, in which the focus was on normal erythropoiesis, linear dependencies of thefeedback controls on the red blood cell population were nevertheless probably not limiting.The question of the relevance of the protein regulatory network considered here may alsobe asked. Erythroid progenitor fate does not depend only on two proteins, Erk and Fas.However, considering a very ‘complete’ regulatory network, with several proteins, does notappear appropriate, at least in order to perform first tests on the hybrid model. Complexitythat would arise from such networks would certainly hide roles and influences of other actors,such as growth factors, cell-cell interaction, feedback controls, etc. The regulatory networkcould however be improved by considering some key proteins or receptors, such as GATA-1,involved in erythroid progenitor differentiation, or EpoR, the erythropoietin receptor, whichplays a crucial role for erythroid progenitor proliferation and inhibition of apoptosis [3, 58].One last point that has not been developed in this study deals with the geometry of thebone marrow. First, particular geometries, with obstacles or describing the porous aspectof the bone marrow, have not been included in this study. The hybrid model can deal withvarious geometries, however the choice was made to focus on the behavior of one single ery-throblastic island and on feedback controls at global and local scales. Second, the hybridmodel was numerically simulated in two dimensions. The case of simulations in three dimen-sions has not been developed. Qualitatively, such a change is not expected to modify theresults, it should nevertheless be tested.Although several aspects of erythropoiesis have not been considered in this model, webelieve the analysis of feedback control role and relevance as well as the investigation ofthe role of a central macrophage in the center of an erythroblastic island improve previousmodels of erythropoiesis and give some insight in the way erythroid cell proliferation anddifferentiation can be controlled, in vitro and in vivo. Further developments are howevernecessary, and confronting the hybrid model with experimental data on stress situations is anext step to validate this new and certainly promising approach.
Aknowledgements
Authors thank Prof. Mark Koury and all members of the INRIA Team DRACULA(http://dracula.univ-lyon1.fr) for fruitful discussions.
A Parameter Estimation
We discuss here the parameter values used when investigating stability of the hybrid modelin Section 3. Some parameters are clearly related to the discrete macroscopic model whilesome others are related to intracellular dynamics.Choosing individual-based modelling creates many implicit parameters. Most of themcan be chosen quite easily (cell radius, parameters of motion, etc.) and are not essential inthe results of simulations. However, other parameters clearly change the overall cell behaviorand have to be set carefully.Let first consider parameters associated with extracellular dynamics (mainly growth fac-tor release and action). Exact quantities of Fas-ligand, Epo and growth factors released bythe macrophage are not fundamental, because their effect will be normalized at the intra-cellular scale (intracellular parameters controlling the sensitivity to these molecules). Moreimportant are the diffusion coefficients, for both Fas-ligand and growth factors released bythe macrophage. These parameters are chosen in order to describe a short range actionof Fas-ligand, and to locate growth factors around the macrophage. It is important thatthe range of diffusion is set according to the size of the cells in order to have a realisticshort-diffusion for instance. The distribution of cell cycle lengths is also really important:18t changes how the island expands, how cell cycles are locally correlated and plays a role inthe fate of immature cells when we work out of equilibrium. Parameter values are shown inTable 1.Intracellular parameters are chosen in order to optimize robustness of the system. It isnoticeable that no experimental data is available to set values of these parameters, only theexpected behavior is accessible. Parameters b and c of system (2.1)–(2.2) control the relativeand absolute speeds of the reactions. They are set first. The magnitude of the hysteresiscycle can then be tuned by choosing carefully α/b and β/b (these ratios control its height,their absolute values its width). Parameter a controls the offset of the cycle and is set suchthat it is located in the [0 , × [0 ,
1] domain. Parameter d mainly controls how a cell goesfrom a stable branch of the cycle to the other: the higher the value of d , the higher thechange in Fas quantity.In equation (2.2), the function γ explicitly describes the influence of Fas-ligand. Thisfunction has been taken as γ = k γ F L , where k γ is the sensitivity of progenitors to Fas-ligand.Finally, α and F cr also depend on extracellular molecules. However, in what followswe will first consider the case when feedback loops are not included, so α and F cr will beconstants first.When Erk and Fas levels are supposed to be zero when the cell starts its cycle, the speedsof reactions are very important: some cells have to remain trapped near the origin at theend of their cycle in order to observe differentiation. The values of b and c are thus veryconstrained. The others simply define a large hysteresis cycle which aims at separating thedifferent cell subpopulations as far as possible from each other. See Table 2 for values.When Erk and Fas values are initially on the hysteresis cycle, then b and c are allowed tovary more largely. It has been decided to put quicker dynamics than in the first case in orderto get rapid variations along the cycle. The width of the hysteresis cycle was also reduced,so values for α and β are higher. See Table 3 for values.Finally, an initial size for the island has to be set and the number of islands that willbe simulated must be chosen. An initial size of 144 cells with central macrophage or 162cells without macrophage was chosen, and parameters were adjusted to keep the size of theisland roughly constant. Therefore, we only simulated one island at the beginning, but oncereasonable parameters are found, it is possible to put several islands side-by-side and see howthey interact. Preliminary experiments have been performed in this direction (not shownhere), see the end of Section 3.1.1. References [1] Ackleh AS, Deng K, Ito K, Thibodeaux J (2006) A structured erythropoiesis model withnonlinear cell maturation velocity and hormone decay rate. Math Bios 204: 21–48.[2] Adimy M, Crauste F, Ruan S (2006) Modelling hematopoiesis mediated by growthfactors with applications to periodic hematological diseases. Bull Math Biol 68 (8):2321–2351.[3] Aispuru GR, Aguirre MV, Aquino-Esperanza JA, Lettieri CN, Juaristi JA, Brandan NC(2008) Erythroid expansion and survival in response to acute anemia stress: the role ofEPO receptor, GATA-1, Bcl-xL and caspase-3. Cell Biol Int 32(8): 966–78.[4] Alarcon T (2009) Modelling tumour-induced angiogenesis: A review of individual-basedmodels and multiscale approaches. In: M.A. Herrero, F. Giraldez, Editors. The Math-ematics of Cancer and Developmental Biology. Contemporary Mathematics 492. pp.45–74. 195] Apostu R, Mackey MC (2008) Understanding cyclical thrombocytopenia: A mathemat-ical modeling approach. J Theor Biol 251: 297–316.[6] Banks HT, Cole CE, Schlosser PM, Hien T (2004) Modelling and Optimal Regulationof Erythropoiesis Subject to Benzene Intoxication. Math Biosci Eng 1(1): 15–48.[7] Bauer A, Tronche F, Wessely O, Kellendonk C, Reichardt HM, et al. (1999) The gluco-corticoid receptor is required for stress erythropoiesis. Genes Dev 13(22): 2996-3002.[8] B´elair J, Mackey MC, Mahaffy JM (1995) Age-structured and two-delay models forerythropoiesis. Math Biosci 128: 317–346.[9] Bernard S, B´elair J, Mackey MC (2003) Oscillations in cyclical neutropenia: New evi-dence based on mathematical modeling. J Theor Biol 223: 283–298.[10] Bernard S, B´elair J, Mackey MC (2004) Bifurcations in a white-blood-cell productionmodel. C R Biologies 327: 201–210.[11] Bessonov N, Crauste F, Fischer S, Kurbatova P, Volpert V (2011) Application of HybridModels to Blood Cell Production in the Bone Marrow. Math Model Nat Phenom 6 (7):in press.[12] Bessonov N, Kurbatova P, Volpert V (2010) Particle dynamics modelling of cell popu-lations. Math Model Nat Phenom 5 (7): 42–47.[13] Chappell, D, Tilbrook PA, Bittorf T, Colley SM, Meyer GT, Klinken SP (1997) Preven-tion of apoptosis in J2E erythroid cells by erythropoietin: involvement of JAK2 but notMAP kinases. Cell Death Differ. 4: 105–113.[14] Chapel SH, Veng-Pedersen P, Schmidt RL, Widness JA (2000) A pharmacodynamicanalysis of erythropoietin-stimulated reticulocyte response in phlebotomized sheep. TheJournal of Pharmacology and Experimental Therapeutics 296: 346–351.[15] Chasis JA, Mohandas N (2008) Erythroblastic islands: niches for erythropoiesis. Blood112 (3): 470–478.[16] Colijn C, Mackey MC (2005) A mathematical model of hematopoiesis – I. Periodicchronic myelogenous leukemia. J Theor Biol 237: 117–132.[17] Colijn C, Mackey MC (2005) A mathematical model of hematopoiesis – II. Cyclicalneutropenia. J Theor Biol 237: 133–146.[18] Crauste F, Pujo-Menjouet L, G´enieys S, Molina C, Gandrillon O (2008) Adding Self-Renewal in Committed Erythroid Progenitors Improves the Biological Relevance of aMathematical Model of Erythropoiesis. J Theor Biol 250: 322–338.[19] Crauste F, Demin I, Gandrillon O, Volpert V (2010) Mathematical study of feedbackcontrol roles and relevance in stress erythropoiesis. J Theor Biol 263: 303–316.[20] Dazy S, Damiola F, Parisey N, Beug H, Gandrillon O (2003) The MEK-1/ERKs signalingpathway is differentially involved in the self-renewal of early and late avian erythroidprogenitor cells. Oncogene, 22: 9205–9216.[21] Demin I, Crauste F, Gandrillon O, Volpert V (2010) A multi-scale model of erythro-poiesis. Journal of Biological Dynamics, 4: 59–70.[22] De Maria R, Testa U, Luchetti L, Zeuner A, Stassi G, et al. (1999) Apoptotic Role ofFas/Fas Ligand System in the Regulation of Erythropoiesis. Blood 93: 796–803.2023] Eller J, Gyori I, Zollei M, Krizsa F (1987) Modelling Thrombopoiesis Regulation - IModel description and simulation results. Comput Math Appli 14: 841–848.[24] Foo J, Drummond MW, Clarkson B, Holyoake T, Michor F (2009) Eradication of chronicmyeloid leukemia stem cells: a novel mathematical model predicts no therapeutic benefitof adding G-CSF to imatinib. PLoS Comput Biol, Sep;5(9):e1000503.[25] Freise KJ, Widness JA, Schmidt RL, Veng-Pedersen P (2008) Modeling time variantdistributions of cellular lifespans: increases in circulating reticulocyte lifespans followingdouble phlebotomies in sheep. J Pharmacokinet Pharmacodyn 35: 285–323.[26] Gandrillon, O (2002) The v-erbA oncogene. Assessing its differentiation-blocking abilityusing normal chicken erythrocytic progenitor cells. Methods Mol Biol 202: 91–107.[27] Gandrillon O, Samarut J (1998) Role of the different RAR isoforms in controlling theerythrocytic differentiation sequence. Interference with the v-erbA and p135gag-myb-etsnuclear oncogenes. Oncogene 16: 563–574.[28] Gandrillon O, Schmidt U, Beug H, Samarut J (1999) TGF-beta cooperates with TGF-alpha to induce the self-renewal of normal erythrocytic progenitors: evidence for anautocrine mechanism. Embo J 18: 2764–2781.[29] Golubev A (2010) Random discrete competing events vs. dynamic bistable switches incell proliferation in differentiation. J Theor Biol 267(3): 341–354.[30] Haurie C, Dale DC, Mackey MC (1998) Cyclical neutropenia and other periodic hema-tological diseases: A review of mechanisms and mathematical models. Blood 92: 2629–2640.[31] Hoehme S, Drasdo D (2010) A cell-based simulation software for multi-cellular systems.Bioinformatics 26(20): 2641–2642.[32] Jeon J, Quaranta V, Cummings PT (2010) An Off-Lattice Hybrid Discrete-ContinuumModel of Tumor Growth and Invasion. Biophys J 98(1): 37–47.[33] Koury MJ, Bondurant MC (1990) Erythropoietin retards DNA breakdown and preventsprogrammed death in erythroid progenitor cells. Science 248: 378–381.[34] Mackey MC (1978) Unified hypothesis of the origin of aplastic anaemia and periodichematopo¨ıesis, Blood 51: 941–956.[35] Mackey MC (1979) Dynamic hematological disorders of stem cell origin. In: G. Vassileva-Popova and E. V. Jensen, Editors. Biophysical and Biochemical Information Transfer inRecognition, Plenum Press, New York. pp. 373–409.[36] Mackey MC, Rudnicki R (1994) Global stability in a delayed partial differential equationdescribing cellular replication. J Math Biol 33: 89–109.[37] Mackey MC, Rudnicki R (1999) A new criterion for the global stability of simultaneouscell replication and maturation processes. J Math Biol 38: 195–219.[38] Mahaffy JM, Belair J, Mackey MC (1998) Hematopoietic model with moving boundarycondition and state dependent delay: applications in erythropoiesis. J Theor Biol 190:135–146.[39] Michor F, Hughes TP, Iwasa Y, Branford S, Shah NP, et al. (2005) Dynamics of chronicmyeloid leukaemia. Nature 435(7046): 1267–1270.2140] Michor F (2007) Quantitative approaches to analyzing imatinib-treated chronic myeloidleukemia. Trends Pharmacol Sci., May;28(5): 197–199.[41] Michor F (2008) Mathematical models of cancer stem cells. J Clin Oncol, Jun 10;26(17):2854–1861.[42] Nagata Y, Takahashi N, Davis RJ, Todokoro K (1998) Activation of p38 MAP kinaseand JNK but not ERK is required for erythropoietin-induced erythroid differentiation.Blood 92: 1859–1869. -[43] Osborne JM, Walter A, Kershaw SK, Mirams GR, Fletcher AG, et al. (2010) A hybridapproach to multi-scale modelling of cancer. Phil Trans R Soc A 368: 5013–5028.[44] Pain B, Woods CM, Saez J, Flickinger T, Raines M, Kung HJ, et al. (1991) EGF-R asa hemopoietic growth factor receptor: The c-erbB product is present in normal chickenerythrocytic progenitor cells and controls their self-renewal. Cell 65: 37–46.[45] Patel AA, Gawlinsky ET, Lemieux SK, Gatenby RA (2001) A Cellular AutomatonModel of Early Tumor Growth and Invasion: The Effects of Native Tissue Vascularityand Increased Anaerobic Tumor Metabolism. J Theor Biol 213: 315–331.[46] Ramis-Conde I, Chaplain MA, Anderson AR, Drasdo D (2009) Multi-scale modelling ofcancer cell intravasation: the role of cadherins in metastasis. Phys Biol 6(1): 016008.[47] Rew DA, Wilson GD, Taylor I, Weaver PC (1991) Proliferation characteristics of humancolorectal carcinomas measured in vivo. Br J Surg 78: 60–66.[48] Rhodes MM, Kopsombut P, Bondurant MC, Price JO, Koury MJ (2008) Adherence tomacrophages in erythroblastic islands enhances erythroblast proliferation and increaseserythrocyte production by a different mechanism than erythropoietin. Blood 111: 1700–1708.[49] Roeder I (2006) Quantitative stem cell biology: computational studies in the hematopoi-etic system. Curr Opin Hematol 13: 222–228.[50] Rubiolo C, Piazzolla D, Meissl K, Beug H, Huber JC, et al. (2006) A balance betweenRaf-1 and Fas expression sets the pace of erythroid differentiation. Blood 108: 152–159.[51] Salazar-Ciudad I, Jernvall J (2010) A computational model of teeth and the develop-mental origins of morphological variation. Nature 464(7288): 583–6.[52] Santillan M, Mahaffy JM, Belair J, Mackey MC (2000) Regulation of platelet production:The normal response to perturbation and cyclical platelet disease. J Theor Biol 206:585–603.[53] Savill NJ, Chadwick W, Reece SE (2009) Quantitative Analysis of Mechanisms ThatGovern Red Blood Cell Age Structure and Dynamics during Anaemia. PLoS ComputBiol 5(6): e1000416.[54] Sawyer ST, Jacobs-Helber SM (2000) Unraveling distinct intracellular signals that pro-mote survival and proliferation: study of erythropoietin, stem cell factor, and constitu-tive signaling in leukemic cells. J. Hematother. Stem Cell Res 9: 21–29.[55] Spencer SL, Gerety RA, Pienta KJ, Forrest S (2006) Modeling somatic evolution intumorigenesis. PLoS Comput Biol 2(8): e108.[56] Spivak JL, Pham T, Isaacs M, Hankins WD (1991) Erythropoietin is both a mitogenand a survival factor. Blood 77: 1228–1233.2257] Scholz M, Engel C, Loeffler M (2005) Modelling human granulopoiesis under poly-chemotherapy with G-CSF support. J Math Biol 50(4): 397–439.[58] Socolovsky M (2007) Molecular insights into stress erythropoiesis. Current opinion inhematology 14: 215–224.[59] Sui X, Krantz SB, Zhao ZJ (2000) Stem cell factor and erythropoietin inhibit apoptosisof human erythroid progenitor cells through different signalling pathways. Br J Haematol110: 63–70.[60] Tsiftsoglou AS, Vizirianakis IS, Strouboulis J (2009) Erythropoiesis: model systems,molecular regulators, and developmental programs. IUBMB Life., Aug;61(8): 800–830.[61] Veng-Pedersen P, Chapel S, Schmidt PRL, Al-Huniti NH, Cook RT, et al. (2002) Anintegrated pharmacodynamic analysis of erythropoietin, reticulocyte, and hemoglobinresponses in acute anemia. Pharm Res 19: 1630–1635.[62] Wichmann HE, Gerhardts MD, SPechtmeyer H, Gross R (1979) A mathematical modelof thrombopoiesis in rats. Cell Tissue Kinet 12: 551–567.[63] Wichman HE, Loeffler M (1985) Mathematical Modeling of Cell Proliferation. BocaRaton, FL, CRC.[64] Wichmann HE, Loeffler M, Schmitz S (1985) A concept of hemopoietic regulation andits biomathematical realisation. Blood Cells 14: 411–429.[65] Wichmann HE, Loeffler M, Pantel K, Wulff H (1989) A mathematical model of ery-thropoiesis in mice and rats. Part 2. Stimulated erythropoiesis. Cell Tissue Kinet 22:31–49.[66] Woo S, Krzyzanski W, Jusko WJ (2006) Pharmacokinetic and pharmacodynamic mod-elling of recombinant human erythropoietin after intravenous and subcutaneous admin-istration in rats. J Pharmacol Exp Ther 319: 1297–1306.[67] Woo S, Krzyzanski W, Jusko WJ (2008) Pharmacodynamic model for chemotherapy-induced anemia in rats. Cancer Chemother Pharmacol 62: 123–133.[68] Wulff H, Wichmann HE, Loeffler M, Pantel K (1989) A mathematical model of ery-thropoiesis in mice and rats. Part 3. Suppressed erythropoiesis. Cell Tissue Kinet 22:51–61. 23 α a d γ Figure 1:
Zero-lines f and g defined in (2.3) and (2.4). Left panel: the f zero-line(straight red line) and the impact of parameters of equation (2.1). Right panel: the g zero-line(straight green line) and the impact of parameters of equation (2.2).24 BCFigure 2: Three possible and characteristic configurations of the system (2.3)–(2.4).
The red line represents the f zero-line, the green line the g zero-line. Arrows indicatedirections of the solutions of system (2.1)–(2.2) in the (Erk,Fas) phase plane. A: Low valueof α , low value of γ . B: Low value of α , medium value of γ (the black curve is the separatrix).C: High value of α compared to β . 25 cr F cr E cr F cr A B E cr F cr CFigure 3:
Illustration of the ‘out of equilibrium’ hypothesis.
All daughter cells areinitially placed at the origin of the (Erk,Fas) plane and, depending on the external Fas-ligandvalue, try to escape the default differentiation zone (bottom left area) into apoptosis (abovethe dashed F cr line) or self-renewal zone (bottom right area) before the end of their cycle.Trajectories from the origin are illustrated by the straight blue curve. A: Low exposure toFas-ligand leads to self-renewal. B: Medium exposure to Fas-ligand leads to differentiation.C: High exposure to Fas-ligand leads to apoptosis.26 poptosisself-renewaldifferentiation Figure 4:
The ‘hysteresis cycle’ scenario.
Daughter cells initially inherit half of theirmother’s protein level, ‘choose’ one of the stable branches depending on their expositionto Fas-ligand, and then move on these branches which form an hysteresis cycle, makingthe decision hardly reversible. Shown with yellow circles are the two possible positions fornewborn cells exposed to exactly the same value of Fas-ligand on the hysteresis cycle. Onecan see how the initial choice is important: cells whose mother had low Fas values will belocated on the lower branch, those whose mother had sufficiently high Fas values will switchto the upper branch. 27
Out of equilibrium' scenario `Hysteresis cycle' scenario
X 00 X X/2X/2YZ ZY
Figure 5:
Schematic representation of the two scenarii considered for distributionof Erk and Fas quantities at cell division.
The left panel illustrates the ‘out of equilib-rium’ case, the right panel the ‘hysteresis cycle’ case. In the ‘out of equilibrium’ scenario acell with a level X = ( E, F ) of Erk and Fas at mitosis gives birth to two daughter cells withno Erk and no Fas ((
E, F ) = (0 , Y and Z .In the ‘hysteresis cycle’ scenario, a cell with a level X = ( E, F ) of Erk and Fas produces twodaughter cells with levels X/ E/ , F/
2) which will reach a priori different levels of Erkand Fas at the end of their own cell cycle, due to stochasticity and response to signaling.28igure 6:
Screenshots of the software program with and without macrophage.
The green cell is a macrophage, yellow cells are immature cells (erythroblasts) and blue cellsare differentiated cells (reticulocytes). The green substance is secreted by the macrophage(growth factors) and the red substance by reticulocytes (Fas-ligand). In each cell, the blackarea represents the incompressible part of the cell.29 N u m b e r o f c e ll s t (cycles)
14 28 N u m be r o f c e ll s t (cycles) Figure 7:
Illustration of the fate of the most stable erythroblastic islands in theabsence of macrophage.
Results were obtained out of equilibrium (left) and using thehysteresis cycle (right). Green lines represent the number of erythroblasts, blue lines reticu-locytes, red lines estimation of circulating red blood cells produced by the island (thick linesare medians, thin lines quartiles over 40 simulations). Observed saturation is only due tospace limitation. 30 .00 0.05 0.10 0.15 0.20 0.25 0.30 0.35 P r o p o r t i o n o f p r o g e n i t o r f a t e F cr .00 P e r c e n t a g e o f p r o g e n i t o r f a t e α Figure 8:
Effect of feedback controls on progenitor subpopulations in the absenceof macrophage in the island.
Each curve represents the percentage (means ± standarddeviation) of self-renewing (green), differentiating (blue) or apoptotic (red) erythroblastsduring the first five cell cycles, as function of F cr (left panel) and α (right panel).31
14 28 N u m be r o f c e ll s t ( cycles ) Figure 9:
Number of erythroid cells in the presence of a macrophage in the centerof the island.
Results were obtained using the hysteresis cycle. Green lines representthe number of erythroblasts, blue lines reticulocytes, red lines estimation of circulating redblood cells produced by the island (thick lines are medians, thin lines quartiles over 40simulations). After addition of a macrophage, islands automatically stabilize for previouslychosen parameters. 32 .1 0.9 N u m be r o f c e ll s Fcr0.2 0.3 0.4 0.5 0.6 0.7 0.8 .0 0 .3 0.6 0.9 1.2 1.5 1.8020040060080010001200 α ₀ N u m be r o f c e ll s Figure 10:
Impact of feedback controls on steady-state values when the erythrob-lastic island contains a macrophage.
Left panel: Influence of F cr variations; Right panel:Influence of α variations. Green lines represent erythroblast counts, blue lines reticulocytecounts and red lines RBC counts. 33arameter Value UnitCell cycle length 24 h Cell cycle variations 4/3 hr Lm Mµ . h − K . M.h − D F L . − L .h − k F L . − molecules.cell − . h − σ F L h − D GF . − L .h − W GF . − molecules. L − . h − σ GF h − Table 1:
Extracellular parameters ( M is an arbitrary mass unit, L an arbitrary lengthunit) 34arameter Value Unit α h − β h − .N U − k a h − b h − .N U − k γ h − .N U − .molecule − c h − .N U − d h − E cr N UF cr N U
Table 2:
Internal parameters when cells begin at the origin (‘out of equilibrium’case).
N U is a normalized quantity unit for the intracellular molecules (maximum possiblequantity is 1). 35arameter Value Unit α h − β h − .N U − k a h − b h − .N U − k γ h − .N U − .molecule − c h − .N U − d h − E cr N UF cr N U
Table 3:
Internal parameters when cells are on the hysteresis cycle.