A combined model of human erythropoiesis and granulopoiesis under growth factor and chemotherapy treatment

Background Haematotoxicity of conventional chemotherapies often results in delays of treatment or reduction of chemotherapy dose. To ameliorate these side-effects, patients are routinely treated with blood transfusions or haematopoietic growth factors such as erythropoietin (EPO) or granulocyte colony-stimulating factor (G-CSF). For the latter ones, pharmaceutical derivatives are available, which differ in absorption kinetics, pharmacokinetic and -dynamic properties. Due to the complex interaction of cytotoxic effects of chemotherapy and the stimulating effects of different growth factor derivatives, optimal treatment is a non-trivial task. In the past, we developed mathematical models of thrombopoiesis, granulopoiesis and erythropoiesis under chemotherapy and growth-factor applications which can be used to perform clinically relevant predictions regarding the feasibility of chemotherapy schedules and cytopenia prophylaxis with haematopoietic growth factors. However, interactions of lineages and growth-factors were ignored so far. Results To close this gap, we constructed a hybrid model of human granulopoiesis and erythropoiesis under conventional chemotherapy, G-CSF and EPO applications. This was achieved by combining our single lineage models of human erythropoiesis and granulopoiesis with a common stem cell model. G-CSF effects on erythropoiesis were also implemented. Pharmacodynamic models are based on ordinary differential equations describing proliferation and maturation of haematopoietic cells. The system is regulated by feedback loops partly mediated by endogenous and exogenous EPO and G-CSF. Chemotherapy is modelled by depletion of cells. Unknown model parameters were determined by fitting the model predictions to time series data of blood counts and cytokine profiles. Data were extracted from literature or received from cooperating clinical study groups. Our model explains dynamics of mature blood cells and cytokines after growth-factor applications in healthy volunteers. Moreover, we modelled 15 different chemotherapeutic drugs by estimating their bone marrow toxicity. Taking into account different growth-factor schedules, this adds up to 33 different chemotherapy regimens explained by the model. Conclusions We conclude that we established a comprehensive biomathematical model to explain the dynamics of granulopoiesis and erythropoiesis under combined chemotherapy, G-CSF, and EPO applications. We demonstrate how it can be used to make predictions regarding haematotoxicity of yet untested chemotherapy and growth-factor schedules.

Influx and efflux of cells at one compartment are amplified, so that the product is the over-all amplification (A in X (t) · A out X (t) = A X (t)), for A X = 0, A X = 1 ln 2 for A X = 1 0 for A X = 0 The effect is a delayed reaction of efflux and compartment size to changes in amplification rates. Amplification splitting is modelled for all compartments with amplification. type/calculation C X content of compartment X function of time t C nor X content of compartment X in steady state parameter, in general we set (normal value) C X (0) = C nor X C rel X content of compartments X relative C rel X (t) = CX (t) C nor X to normal value C in X influx in compartment X function of time C in nor X normal influx parameter, see above C out X efflux from compartment X function of time C out nor X normal efflux parameter, see above a X proliferative fraction in cell compartment X function of state, sometimes constant A X amplification in cell compartment X " A in X amplification of influx " A in X amplification of effflux " n X average number of mitoses n X = ldA X in cell compartment X p self-renewal probability of stem cells function of state τ X average duration of cell cycle function of time, sometimes in compartment X constant (not regulated) T X average transit time of active cells in cell T X = n X τ X compartment X T t X total transit time T t X = nX τX aX k transition, degradation or toxicity coefficients functions of time or parameter Y min quantity Y under minimum stimulation parameter to determine the regulatory function of Y Y nor quantity Y in steady state " Y int quantity Y under intensified stimulation " Y max quantify Y under maximum stimulation " b Y sensitivity of Y under stimulation "

A.2 Regulatory Function
Amplification A and transition time T are regulated between a minimum Y min and a maximum Y max by the growth factors EPO or G-CSF according to the following regulatory function: where Y nor is the steady state value of amplification or transition time, C Cyto ∈ {C int EPO , C rel GCSF } and b Y is the sensitivity of Y under stimulation (see figure A.2 and [1], p. 69).

A.3 Modelling of Delays
Several delays are included into the model such as delayed action of growth factors or chemotherapy. According to [2,3], these delays are modelled by a set of concatenated compartments with first order transitions where C X is a quantity to be delayed, D X is the delay parameter and N X is the number of delay compartments. The delayed quantity is now defined by i.e. as the efflux from the last delay compartment.
As explained in [2,3] this modelling of delays is in between a random, age-independent transition (N X = 1) and a strict "first-in-first-out" kinetic (N X → ∞). By our modelling, one can mimic both, an expectation and a variance of individual delay times. This is achieved by setting D X = T X /N X where T X is the desired expectation of the transition time and N X determines its variance (see [2] for details).

A.4 Chemotherapy
The infusion of chemotherapeutic drugs is modelled by Heaviside functions where N cycle is the number of chemotherapy cycles,t i are the time points with application of a specific drug, and t CHEMO inf is the duration of chemotherapy application. The delayed effect of the drugs is modelled by four compartments drug out (t) − k drug Delay · Ψ  The output-function Ψ (4) drug out (t) is multiplied by the toxicity parameters of the single compartments k drug S , k drug CG , k drug PGB , k drug MGB , k drug BE , k drug CE , k drug PEB , k drug MEB , and k drug RET respectively. If muliple cytotoxic drugs are applied, corresponding toxicity functions are added resulting in an overall toxic effect Ψ X which is cell stage and chemotherapy specific. In complete analogy to our former work [3,4], the effect of chemotherapy is introduced to the balance equations of the bone marrow cell compartments by a first-order loss term. Hence, the schematic compartment equation has the form In clinical practice often only leukocytes are available. We calculate the leukocyte count as the sum of lymphocytes and granulocytes. To avoid a full model of lymphopoiesis, we modelled the reduced lymphocyte count under chemotherapy by an exponential function of the corresponding toxicity function.
where c LY = 3000 cells per µl and c GRA = 4000 cells per µl are the normal concentrations of lymphocytes and granulocytes respectively. Ψ LY is the toxicity function for lymphocytes defined in analogy to the toxicity functions of bone marrow cell stages (see above).

A.5 Stem cell compartment S
Cells differentiating into granulopoietic or erythropoietic lineages originate from the same stem cell compartment. We adopted the corresponding stem cell model of Loeffler & Wichmann [1]. The output of the stem cell compartment is splitted into red or white cell lines under the assumption that 15% of the released cells differentiate into red blood cells (α E = 0.15), and 80% into the white blood cell line (α G = 0.8). The stem cell compartment S has self-renewal capability. Under steady state conditions, 50% of the stem cells remain in this compartment, and 50% differentiate into red or white blood cell lineages. Hence, the stem cell compartment equation is where a S is the proliferative fraction, τ S is the average duration of a cell cycle, and p the self-renewal probability of stem cells.
The self renewal probability p is regulated by a competition of the stem cells (C rel S (t)), the granulopoietic (C rel G (t)), and the erythropoietic (C rel E (t)) bone marrow cells: where the parameters ϑ E = −2, ϑ G = −8, and ϑ S (t) are hypothetical weighting factors representing the strength of the influence of the bone marrow cells C rel S (t) = CS (t) C nor S , C rel E (t) = CBE(t)+CCE(t)+CPEB(t)+CMEB(t) C nor BE +C nor CE +C nor PEB +C nor BE , and C rel G (t) = CCG(t)+CPGB(t)+CMGB C nor CG +C nor PGB +C nor MGB . It is assumed that where the steady state value p nor = 1 2 . Thus, for the initial conditions it holds that The proliferative fraction a S can be interpreted as the percentage of cells which are currently in cell cycle.
The proliferative fractions a X of the compartments S, BE and CG are regulated by the haematopoietic bone marrow system C rel S (t) and C rel G (t) or C rel E (t): The parameters ω E = 0.3, ω G = 0.1 and ω S = 1 are weighting factors. They represent the strengths of the influence of stem cells, erythropoietic and granulopoietic cells on the proliferative fraction of BE, CG and S. With where X is one of the compartments S, BE, or CG. It is a monotone function with ranges between a min X and a max X . Low cell numbers in the bone marrow compartments cause a higher demand of proliferating cells and therefore a larger proliferative fraction a X . The value of y defines the actual point on the regulatory curve.
The variable x is a measure of the total bone marrow content. It is calculated as a weighted sum of the logarithms of the relative counts of stem cells, erythropoietic cells and granulopoietic cells. If any cell counts tend to zero, x tends to minus infinity, and with it, a becomes maximal. Parameter values a int corresponds to x = − ln 2 and a nor corresponds to x = 0 (see figure 2).

Compartment BE
The model equations of BE are similar to [4], but now the additional influence of the granulopoietic cell lineage on the proliferative fraction in BE is taken into account. with the initial values where α E = 0.15 is the proportion of cells differentiating into erythropoietic cell lineage [1].

Compartment CE
The initial values are The amplification and transit time in the compartment CE are regulated by the growth factors EPO and G-CSF. This is modelled by regulatory functions acting on proliferation rate and transition time in the 7 compartment CE. The regulatory function of amplification is regulated by the internalised EPO Z ACE C int EPO . The regulatory function of the transition time is also regulated by internalised EPO, Z TCE C int EPO . But, it is additionally multiplied by a regulatory function regulated by G-CSF (F GCSF (t)), mimicking the G-CSF effect on CE. In summary, amplification and transition time in CE are calculated by where Z T Peg (t) and Z T Fil (t) are the regulation functions of the growth factors endogenous G-CSF and Filgrastim on one hand and Pegfilgrastim on the other hand (equation A.4), 0 ≤ ω P (t) ≤ 1 is the weighting factor to model the superimposing effect of concurrent Filgrastim or endogenous G-CSF, and Pegfilgrastim (see [3]), and C int EPO is the internalised EPO (see [4]): with ω min P = 0 and ω max P = 1

Compartment PEB
The initial values are 8

Compartment MEB
The maturation is modelled by splitting MEB into N MEB = 15 subcompartments, without amplification.
The initial values are

Compartment RET
q RET is the proportion of reticulocytes to the total number of red blood cells in steady state. s nor ERY , T ERY rnd , and T ERY age were explained in the next section.

Compartment ERY
The compartment ERY is split into the compartments "RANDOM" and "AGE". In steady state, most erythrocytes die dependent on age. The age dependent reduction is modelled by divisions into subcompart-9 ments.
Under stimulation, the depletion is more randomly (i.e. exponential decay, see [5]). Hence, the influxes into the compartments "RANDOM" and "AGE" are regulated by the factor s ERY , which depends on the bone marrow output of the reticulocytes. T ERY rnd , and T ERY age are transition times in the subcompartments "RANDOM" and "AGE". (See [1,2,5,6].) with initial conditions

Endogenous production of EPO
According to [6,7], the endogenous production of EPO (EPO prod ) is assumed to depend on the oxygen partial pressure in the kidneys and the number of circulating red blood cells The variables are explained in Table A.2. For further explanation and justification, see [4,6,7].

Exogenous EPO Application
In our model, we use simple pulse functions EPO inj for intravenous injections. Regarding subcutaneous injections, the model adapted from [8] includes direct absorption from the subcutaneous tissue into the bloodstream, or indirect through the lymphatic system. In both processes, a time delay is assumed. A loss of EPO is included at the injection site (k F e ) and in the lymphatic system (k L e ). The structure of this injection model is described in detail in [4,8]. The model equations read as follows: The general EPO injection function EPO inj (t) is modelled by a sum of pulse functions is the Heaviside-function,t i are time points at which EPO at dose EPO dose is administered. The injection time EPO t inf is set to five minutes, and EPO dose is the administered dose in IU/kg. The dynamics of the EPO concentration in the subcutaneous tissue C SC EPO (t) is described by: where k F a is the absorption constant for the direct influx into the central compartment, k F e is a loss term at injection site, and k FL is the absorption constant of the lymphatic system. The efflux from the subcutaneous tissue into the peripheral blood is delayed by four delay compartments, i.e.
with the settings EPO F out (t) enters the central EPO compartment. Analogously, the lymphatic absorption is modelled by a delay function with four compartments: where C EPO L out (t) enters the central compartment. Thus, EPO dynamics in the lymphatic compartment C L EPO are given by where k L a is the absorption constant of the transition between the lymphatic and the central compartment and k L e is a loss term. Summarising the effluxes of the direct and the lymphatic way of absorption yields for intravenous injections.

Central EPO Compartment
Dynamics of Erythropoietin in central compartment (circulation) is modelled in the following way (see [1,4,9]): where EPO prod equals one in steady state, so that equation A.25 equals zero.
The peripheral compartment and the EPO-receptor complex are described by The dynamics of the EPO receptors R is determined on the basis of the bone marrow content of the cell kinetic model, namely the compartments C BE , C CE , C PEB , C MEB , and C RET [4].
To account for different receptor densities of erythropoietic bone marrow cells we introduced weighting factors w RET , w MEB , w PEB , w CE , and w BE . CFU-E have the highest weighting factor due to the highest number of EPO receptors observed [10]: w BE ≤ w CE . The receptor density declines with further maturation: where R rel (t) = wRET·CRET+wMEB·CMEB+wPEB·CPEB+wCE·CCE+wBE·CBE wRET·C nor RET +wMEB·C nor MEB +wPEB·C nor PEB +wCE·C nor CE +wBE·C nor BE is the number of EPO receptors relative to steady-state [4].
The initial values are derived from steady-state conditions where C cent EPO (0) = EPO serum · EPO Vc , EPO serum = 15 IU/l is the basic level of endogenous EPO, and EPO Vc denotes the distribution volume of EPO [9]. According to [9], we set R (0) = 64.31 and from equation A.27 it follows that k syn = k deg · R(0) + k int · C rb EPO (0). The relative internalised EPO (C rel EPO int ) is used as argument of the regulatory functions regulated by EPO:

A.7 Granulopoiesis
Here we present model equations of our granulopoiesis model (see also [3]). Initial conditions are derived again from steady-state conditions.

Compartment CG
The model equations of CG are similar to [3]. The influence of red blood cell line on the proliferative fraction is modelled in analogy to [1].
where a CG is the proliferative fraction and α G = 0.8 is the part of cells differentiating into the white blood cell line [1]. The amplification A CG and the transit time T CG are regulated by the concentration of G-CSF in the central compartment.

Compartment MGB
This compartment is divided into three compartments denoted as G4, G5 and G6. The compartments are again divided into N X subcompartments to model the maturation process by a delay. In these subcompartments, the effect of postmitotic apoptosis is also implemented [2,11] by introducing a postmitotic amplification detoned again as A. It holds that A ≤ 1 for all subcompartments.
The total postmitotic amplification A G4 is equally distributed over all subcompartments in which there is postmitotic amplification, e.g. if there is postmitotic amplification in all subcompartments, it holds that For G5 and G6 the equations are completely analogous to (A.35)-(A.41) if one replaces PGB by G4 and G5 respectively. In the present form of our model, postmitotic amplification is restricted to G6.
where Ψ Pred is the characteristic function of Prednisone applications modelled as a step-function which is equal to one for the duration of one day after Prednisone application.

G-CSF
According to [2], the relative G-CSF production P endo GCSF is a regulatory function of the relative content of segmented granulocytes in bone marrow and granulocytes in circulation.
where ω G6 and ω GRA are weighting parameters and P endo nor The G-CSF injection function reads as follows: where Hv(t) is the Heaviside-function (equation A.18), t i ≥ 0 (i = 1, . . . , L) are the time points of G-CSF injections and d GCSF (t i ) are the corresponding doses (in µg). The injection time t inf is set to 5s. The injection function is specific for each G-CSF derivative (details see [3]).
The subcutaneous compartment is divided into two subcompartments sc 1 and sc 2 where the efflux of the first compartment is the influx to the second compartment. G-CSF is applied to the first subcompartment (second term of A.47). In the first subcompartment there is a dose-dependent loss of G-CSF modelled by a

Michaelis-Menten kinetic (third term of A.47). For Filgrastim injections it holds that
with the initial values C sc 1 GCSF (0) = C sc 2 GCSF (0) = 0. For Pegfilgrastim injections the first term of the righthand side of A.47 is substituted by P exo peg sc GCSF . Likewise, the Filgrastim parameters k F sc , v F max and k F m are substituted by corresponding Pegfilgrastim parameters.
For Filgrastim injections it holds that The balance equation A.49 contains the following terms: P ref GCSF P endo GCSF is the endogenous production, GRA is the specific elimination. The corresponding equation for Pegfilgrastim is the same except for the endogenous production, which is zero, the intravenous injection function which is substituted by P exo peg iv GCSF , the parameters and the initial value which is again zero.
For both G-CSF derivatives, we have transitions between central and peripheral compartment: where the parameters k cp , k pc and V D are specific for Filgrastim and Pegfilgrastim respectively. To model the competition of Pegfilgrastim and endogenous G-CSF with respect to receptor binding, the regulatory functions of Pegfilgrastim and Filgrastim were again combined using the weighting factor ω P defined in equation A. 16: where Y is an arbitrary quantity regulated by G-CSF such as transition times or amplifications. For all these quantities we assumed the same regulatory function of the weighting parameter ω P .

A.8 Parameters
Here, we present all parameters of the model, their values, and how they were determined. In table A.3 we present general parameters of the cell kinetic model of granulopoiesis.  A number of parameters depending on EPO derivative but not on mode of application is presented in Chemotherapy parameters are specific for cell stage, drug, drug doses, and sometimes, risk group of patients such as young or elderly patients. An overview of estimated parameters can be found in tables A.8 (toxicity to stem cells and granulopoiesis) and A.9 (toxicity to erythropoiesis and lymphopoiesis).

Stability analysis
Using the above mentioned model parameters results in a stable steady state of the system. However, modifying certain parameters can result in stable oscillations. This especially applies for parameters of the long range stem cell feedback (see figure 3).    Figure 4: Sensitivity of the parameters describing the regulation in compartment CE. Parameters denoted with "*" could not be changed by 2.5 % without violating constraints. Parameters were modified by ±2.5%. Relative change of the fitness function is shown as length of corresponding bars.

A.9 Simulation Results
In this section, we present all simulation results of our hybrid model and compare it with available data. We also present extensive comparisons of our hybrid model with the single lineage models of erythropoiesis and granulopoiesis.               Figure 21: Serum concentration of erythropoietin, HB and HK value, simulation (black line) and data (circle), data: [17], [57], from subjects required to have a baseline HK less than 48%. If the HK rose above 55%, phlebotomy was performed.    Figure 28: simulation (black line) and data (circle), data: [43,44]. The blood cell donation in the scenario of Bensinger et al. [44] is not taken into account.                [34,36,37,58]. Caused by Hodgkins lymphoma activity, an increased endogenous G-CSF production is assumed in BEACOPP treated patients, see [2].