Quantifying Drug-Induced Bone Marrow Toxicity Using a Novel Haematopoiesis Systems Pharmacology Model.

Haematological toxicity associated with cancer therapeutics is monitored by changes in blood cell count, and their primary effect is on proliferative progenitors in the bone marrow. Using observations in rat bone marrow and blood, we characterize a mathematical model that comprises cell proliferation and differentiation of the full haematopoietic phylogeny, with interacting feedback loops between lineages in homeostasis as well as following carboplatin exposure. We accurately predicted the temporal dynamics of several mature cell types related to carboplatin-induced bone marrow toxicity and identified novel insights into haematopoiesis. Our model confirms a significant degree of plasticity within bone marrow cells, with the number and type of both early progenitors and circulating cells affecting cell balance, via feedback mechanisms, through fate decisions of the multipotent progenitors. We also demonstrated cross-species translation of our predictions to patients, applying the same core model structure and considering differences in drug-dependent and physiology-dependent parameters.

WHAT IS THE CURRENT KNOWLEDGE ON THE TOPIC?

☑ Cancer therapeutics impact on proliferative progenitors in the bone marrow, disrupting the process of haematopoiesis and leading to observed changes in blood cell counts. Theoretical models have been used to characterize these variations in peripheral blood cells using assumptions based on nonquantified bone marrow toxicity.

WHAT QUESTION DID THIS STUDY ADDRESS?

☑ This study aims to mechanistically characterize the overall carboplatin‐induced myelosuppression profile, linking carboplatin pharmacokinetics with variations in bone marrow progenitors and in circulating blood cells.

WHAT DOES THIS STUDY ADD TO OUR KNOWLEDGE?

☑ Our model describes, quantitatively and mechanistically, carboplatin‐induced DNA damage on bone marrow progenitors and how toxicity propagates across cell lineages over time, modifying multiple clinical endpoints in peripheral blood.

HOW MIGHT THIS CHANGE DRUG DISCOVERY, DEVELOPMENT, AND/OR THERAPEUTICS?

☑ Our model captures the complexity and dynamics of haematopoiesis and reveals new understanding essential for better prediction of drug‐induced effects on specific lineages in rat and in man. This approach may be adopted for other drugs where myelosuppression is of concern.

Drug‐induced myelosuppression is a common dose‐limiting adverse event for cancer therapeutics,1 causing delays and discontinuation of treatment, which compromise the therapy's success. The cytotoxic activity of antineoplastic agents affects highly proliferative bone marrow progenitors, reduces bone marrow cellularity, and disrupts blood cell production.2 Carboplatin, which is the standard of care treatment for several solid tumors,3 induces DNA damage in proliferating bone marrow cells, causing cell death and consequent myelosuppression.1 The dose‐limiting toxicities in patients receiving carboplatin are thrombocytopenia and leukopenia.4, 5

Haematopoiesis is an extremely dynamic process that comprises highly complex cell differentiation cascades. Around 109 red blood cells and 108 white blood cells are generated every hour from a common set of haematopoietic stem cells (HSCs) in the bone marrow of a healthy adult.6 To enhance our understanding of this process, as well as of the mechanisms behind drug‐induced myelosuppression, it is essential to quantify cell dynamics and interactions across multiple branches of the haematopoietic phylogeny rather than modeling the differentiation of cell lineages in isolation.7 This task is not trivial because the complexity of the system's dynamics arises from the interaction of multiple cellular prosurvival mechanisms primed by suppressive drug effects.

Clinically, drug‐induced myelosuppression is monitored by tracking the variation in blood counts,8 which are relatively easy to obtain. However, haematological toxicity originates in the bone marrow, where proliferative progenitors reside. Theoretical models have been used widely to explore hypotheses on nonquantified drug‐induced toxicity in the bone marrow,7, 9 and empirical adjustments are sometimes required to capture drug effects over multiple cycles.7, 10, 11, 12 In a recent review,7 we discussed the possibility of enhancing clinical data with complementary information from preclinical studies, where measurements in the bone marrow are more feasible.

Here we present a novel quantitative systems pharmacology model to describe the link between carboplatin pharmacokinetic (PK) profiles and variations in bone marrow progenitors and in multiple lineages of mature blood cells (Figure  ). Combining this theoretical framework with experimental results in rats, we investigate hypotheses on progenitor cell fate, feedback mechanisms, and the disruption of cell dynamics in the haematopoiesis cascade disrupted by carboplatin. Furthermore, we translate our preclinical findings and predictions to patients by scaling the essential features of the haematopoietic process and adjusting the drug potency for interspecies variations7, 13 (Figure  ).

Figure 1

Carboplatin toxicity in the bone marrow is mechanistically linked with effects in the blood. The diagram provides a graphical representation of our quantitative systems pharmacology model defined in Eqs. (1), (2), (3), (4), (5), (6), (7), (8), (9). The formation of every mature blood cell is pictured as a path through the haematopoiesis differentiation tree, with haematopoietic stem cells at its top. Haematopoietic stem cells mature into MPPs, which then commit to specific lineages giving rise to fully differentiated blood cells. Red triangles highlight populations affected by carboplatin. Rat and human symbols mark model variables for which observations in rat and human were available. Cell surface markers used to quantify bone marrow progenitors in rats are also reported in the diagram. Proliferative cells are affected by carboplatin. Proliferation and maturation in the bone marrow are regulated by the interplay of feedback mechanisms from different circulating cells (dotted arrows). CD90+/Lin‐, CD90 positive and lineage negative cells; CMP, common‐myeloid progenitor; MEP, megakaryocyte–erythrocyte progenitor; Mono, monocytes; Neut, neutrophils; MPP, multipotent progenitor; Plt, platelets; RBCs, red blood cells; Ret, reticulocytes; T1, T2, T3, transit‐compartments.

Figure 2

Graphical representation of our theoretical approach developed to investigate drug‐induced haematological toxicities. The first step consisted of defining a novel global and quantitative system pharmacology model able to integrate data from different sources (i.e., rat carboplatin PK, bone marrow effects, and peripheral blood counts) and describe carboplatin‐induced myelosuppression profiles in rats. Then, we considered the cross‐species differences between rat and human to update model parameter values. Last, we generated clinical predictions. In addition, when clinical data are available, back‐translation can also be performed (dotted gray arrow). BM, bone marrow; PK, pharmacokinetic.

Methods

In vivo rat studies

All animal care and experimental procedures used at the AstraZeneca facility in the United Kingdom were performed under the authority of a valid Home Office Project License and they were conformed to the UK Animals (scientific procedures) Act, 1986. Animal studies are reported in compliance with the ARRIVE (Animal Research: Reporting of In Vivo Experiments) guidelines.14, 15

A total of 178 animals were used in these experiments. Details of the study design, haematology analysis, flow cytometry analysis, and bioanalysis are reported in Table and Supplementary Text , respectively.

Theoretical description of carboplatin‐induced myelosuppression

Details on carboplatin pharmacokinetic model are reported in Supplementary Text .

We described haematopoiesis as a set of compartments connected in branches (one for each lineage), which originate from common sources of HSCs and multipotent progenitors (MPPs; Figure  ). Branches represent distinct cell lineages (neutrophils, monocytes, platelets, erythrocytes), and form our model subunits. Note that we did not consider the lymphoid branch because the proliferation, activation, and regulation of lymphoid cells occur mostly in lymphoid organs.16 Within lineages, compartments correspond to specific haematopoietic cell stages (progenitors, transit cells, mature cells), and model variables describe the cell count in each of these compartments.17 Our model describes a dynamic system of cells continuously proliferating and differentiating, able to maintain and recover homoeostasis after perturbations by balanced cell flows across the connected compartments (predecessor–successor populations). The value of these flows is determined from the lifespan of cells reported in the literature and observed size of compartments.

We described haematopoietic cell dynamics with a system of ordinary differential equations, which are here grouped together by the stage of the haematopoietic cell population under analysis (namely, progenitors, transamplifying cells, and mature cells).

Progenitor equations

We assumed that all haematopoietic cells come from a common pool of MPPs (CD90+/Lineage– cells):

A source of HSCs6 (k stem), whose flow is mediated by a feedback function (fdbkstem), feeds into the MPP compartment, and MPPs also increase their number via proliferation (). MPPs differentiate into the more mature common myeloid progenitors (CMPs; CD45+ cells) and megakaryocyte–erythrocyte progenitors (MEPs; CD71+ cells), with the rates , respectively. MPPs are affected by carboplatin, whose effect is described by the parameter SlopeMPP. We assumed that the drug effects were driven by the cell damage (Damage) and that drug‐induced damage affected the proliferative fraction of the cell compartment, which is denoted by the product of the proliferation rate (k prol) and the compartment size (MPP):

We quantified cell damage as the result of (i) the process of adducts formation, or DNA Damage, with rate k dam proportional to the concentration of free platinum in plasma, and (ii) the DNA repair process, with first‐order kinetics (k rep), as follows:

The CMP and MEP compartments follow similar rules to those previously described for MPP:

They proliferate with rates , respectively, and differentiate into the successors of their lineage, i.e., transit‐neutrophils (with rate ) and transit monocytes (with rate ) for CMP and transit reticulocytes (with rate ) and transit platelets (with rate ) for MEP. Differentiation rates (k tr) were derived resolving the steady‐state equations reported in Supplementary Text . Similar to MPP, CMP and MEP are affected by the cytotoxic activity of carboplatin (SlopeCMP, SlopeMEP).

Transit compartment equations

We used a chain of three non proliferative transit compartments to describe maturation of both neutrophils (T1Neut, T2Neut, T3Neut) and monocytes (T1Mono, T2Mono, T3Mono):

As previously reported,18 we introduced a new parameter (a) to describe maturation rate within each lineage‐specific transit compartment. We defined this rate as the ratio between the number of transit compartments and the maturation time:18

For reticulocytes and platelets, we used a chain of two proliferative transit compartments (T1Ret, T2Ret, T1Plt, T2Plt) or late progenitor cell populations, thus susceptible to carboplatin activity, followed by non proliferative transit compartments (T3Ret, T3Plt):

As previously described for neutrophil and monocyte maturation, the MTTs of reticulocytes and platelets were defined as:

Lifespans of transit compartments (1/a) and progenitor cells (1/k tr) were obtained from the steady‐state equations (Supplementary Text ) and by fitting the model to the data (Supplementary Text ).

The cells in the proliferative transit compartment are affected by the cytotoxic activity of carboplatin with rates (δRetSlopeMEP, δPltSlopeMEP), which are proportional to the damage induced by carboplatin (Damage) and the proliferative fraction of the compartment.

Circulating cell equations

After maturation, fully differentiated cells are released into the blood:

Circulating rates (k circ) were derived from cell lifespan values (1/k circ) in the literature.

All modeling and simulations were performed in Matlab (Mathworks, Natick, MA) toolbox Simbiology (constrained nonlinear least‐square optimization) and Phoenix non linear mixed effect (NLME) 7.0 (Certara, Princeton, NJ; Supplementary Text and ).

Identification of values for the model parameters

The lifespans of mature blood cells were derived from rat values found in the literature.5, 13, 17, 19 The size of all compartments, lifespan and proliferation of transit compartments, and progenitor cells were calculated by solving the steady‐state equations of the untreated bone marrow (Supplementary Text ) and by fitting the model to our experimental data sets (Supplementary Text ). The parameters describing feedback strength and drug effects were calculated by fitting the model to our experimental data sets (Supplementary Text ). Stability analysis was also performed with Oscill8 Dynamical Systems Toolset (Copyright (C) 2005, Emery Conrad, http://oscill8.sourceforge.net/), Supplementary Text , to identify parameter regions leading to oscillatory solutions, which do not describe a physiological response of haematopoiesis to carboplatin treatment.

Model translation to the clinic

The parameters of our model were scaled to reflect the differences between rat and human (i) haematopoiesis physiology and (ii) carboplatin bone marrow toxicity.

where correcting for the cross‐species differences in the fraction of unbound drug was not needed for carboplatin.

The lifespan and the counts of peripheral blood cells were obtained from the literature, whereas for MPP counts we applied allometric scaling from rat to human20 given the specific marker‐dependency of flow cytometry data. Then we derived CMP and MEP counts using the proportions between these populations and MPP, which are experimentally reported for humans21 (Supplementary Text ). Mean transit times in humans were calculated from the literature5 according to the number of non proliferative transit compartments in our model.

The drug‐dependent parameters, denoted Slope, were scaled using species‐specific in vitro drug potency experimental results (Table and Supplementary Text ), which express the drug effect on several proliferative cell types. To scale these parameters, we used a previously reported expression13 with the following form:

Last, the human carboplatin PK was simulated using a previously reported model22 (Supplementary Text ).

Results

A new quantitative systems pharmacology model for haematopoiesis

Our novel approach describes cell dynamics in the haematopoietic system based on the proliferation, differentiation, and interaction of feedback processes across some of the major branches of the haematopoietic phylogeny.23 This theoretical framework was able to quantify the temporal dynamics of carboplatin‐induced DNA damage on progenitor cells in the bone marrow, which is the drug site of action, together with the downstream impact of this bone marrow toxicity on mature blood cells.

The final model structure contains the formation of distinct blood cell lineages as branches of transit compartments originated from a common source of early progenitors (Figure  ). Each branch incorporates the dynamics of cell proliferation, lineage differentiation, and negative feedback processes in response to stimuli.7 These feedback loops connect cells across lineages and regulate the proliferation and differentiation of common progenitors to rescue the depleted populations.24 This results in a haematopoietic dynamic system in homeostatic equilibrium capable to recover25, 26 after perturbations (Figures  , , ). Parameter values are reported in Tables  , , , and .

Figure 3

Overall myelosuppression profile induced by repeated cycles of carboplatin 40 mg/kg. Rats were dosed with carboplatin 40 mg/kg on day 1 of repeated 14‐day cycles (Table ). Plots show the time course of carboplatin effects across the main bone marrow progenitor populations (multipotent progenitors, common‐myeloid progenitors, megakaryocyte–erythrocyte progenitors) and mature blood cells (neutrophils, monocytes, platelets, reticulocytes, red blood cells). Multipotent progenitor, common‐myeloid progenitor, megakaryocyte–erythrocyte progenitor, neutrophil, monocyte, and red blood cell total counts go down over multiple cycles, whereas reticulocyte and platelet recovery to baseline (after rebound) is still possible. Points show the observations, red lines show average model predictions, and blue shadow areas show percentiles in the simulated data. Visual predictive checks were generated from 1,000 simulations using the log‐additive residual errors reported in Table . d, days.

Table 1

Physiology‐dependent system parameters

Description	Parameter	Unit	Rat			Human
			Value	CV (%)	Source	Value	Source
Circulating rates
Neutrophil circulating rate	k_circ_Neut	1/hour	0.17	–	45	0.1	22
Monocyte circulating rate	k_circ_Mono	1/hour	0.06	–	45	0.04	41
Platelet circulating rate	k_circ_Plt	1/hour	0.01	–	46	0.0052	5
Red blood cell circulating rate	k_circ_RBC	1/hour	0.0007	–	45	0.00037	13
Baseline values
Neutrophil baseline value	Neut0	109 cells/L	1.32	2.6	Estimated	2–7	22
Monocyte baseline value	Mono0	109 cells/L	0.16	2.6	Estimated	0.2–10	41
Platelet baseline value	Plt0	109 cells/L	967	2.5	Estimated	345	5
Reticulocyte baseline value	Ret0	109 cells/L	235	2.7	Estimated	40–115	1
Red blood cell baseline value	RBC0	109 cells/L	8,161	2.2	Estimated	4,100–5,900	32
Multipotent progenitor baseline value	MPP0	109 cells/L	69	4.6	Estimated	1.3	Allometric scaling ref. 42
Common myeloid progenitor baseline value	CMP0	109 cells/L	214	4.1	Estimated	20.9	Proportions from ref. 43
Megakaryocyte–erythrocyte progenitor baseline value	MEP0	109 cells/L	156	5.2	Estimated	15	Proportions from43
Transient times
Neutrophil mean transient time	MTTNeut	hour	61.5	56	Estimated	210	Scaled from ref. 5
Monocyte mean transient time	MTTMono	hour	62.2	55	Estimated	121.5	Scaled from47
Platelet mean transient time	MTTPlt	h	74	5.6	Estimated	168	Scaled from5
Reticulocyte mean transient time	MTTRet	hour	75	4.8	Estimated	66	Scaled from internal AZ study (data not shown)
Steady‐state equation parameters
Reticulocyte 1st transit compartment baseline value	λ1	Dimensionless	2	–	Fixed	2	As in the rat
Platelet 1st transit compartment baseline value	λ2	Dimensionless	2	–	Fixed	2	As in the rat
Proliferation rate CMPs	λ3	Dimensionless	1.8	–	Fixed	1.8	As in the rat
Proliferation rate MEPs	λ4	Dimensionless	1.8	–	Fixed	1.8	As in the rat
Proliferation rate MPPs	λ5	Dimensionless	1.8	–	Fixed	1.8	As in the rat
Feedback powers controlling
MPP zero order proliferation rate	γstem	Dimensionless	0.07	–	Fixed	0.07	As in the rat
MPP maturation rate toward CMPs	γ_Mat_CMP	Dimensionless	0.6	–	Fixed	0.6	As in the rat
MPP maturation rate toward MEPs	γ_Mat_MEP	Dimensionless	0.3	–	Fixed	0.3	As in the rat
First‐order proliferation rates of T1Ret, T2Ret, T1Plt, T1Plt	γ_prol_Trans	Dimensionless	0.7	–	Fixed	0.7	As in the rat
Damage
Damage creation	k dam	1/hour	0.017	–	Fixed	0.017	As in the rat
Damage repair	k rep	1/hour	0.017	–	Fixed	0.017	As in the rat

Values for rat were obtained from the literature, when available; estimated by fitting the model to the data; or derived/fixed from experiments. Human values were derived from literature or scaled as explained in the main text and in the Supplementary Material. Coefficient of variation (CV%) values show the confidence in our parameter estimates. AZ, AstraZeneca; T1, T2, proliferative transit‐compartments.

Table 2

Drug‐dependent system parameters

Description	Parameter	Rat			Human
		Unit	Value	CV (%)	IC50 scaled value
Carboplatin effects in multipotent progenitors	SlopeMPP	1/μM	2.05	11.8	0.79
Carboplatin effects in common myeloid progenitors	SlopeCMP	1/μM	1.47	8.6	0.57
Carboplatin effects in megakaryocyte–erythrocyte progenitors	SlopeMEP	1/μM	2.19	10.3	0.66
Carboplatin effects in reticulocyte proliferating transit‐compartments	δRet	Dimensionless	2.8	21	2.8 (as in rat)
Carboplatin effects in platelet proliferating transit compartments	δPlt	Dimensionless	0.54	30.4	0.54 (as in rat)

Carboplatin effects in rat were estimated by fitting the model to the data. Human values were calculated by adjusting the rat values for species‐specific drug sensitivity , Eq. (10). The IC50 values used to scale drug effects are reported in Table . Coefficient of variation (CV%) values show the confidence in our parameter estimates. IC50, half maximal inhibitory concentration.

Carboplatin‐induced DNA damage drives proliferative cell loss

We observed that nucleated cells in the bone marrow started to decrease from the first administration of carboplatin. The nadir was observed 3–4 days post dosing, when free platinum, which is the active species of carboplatin and that is cleared relatively rapidly from blood27, 28 (Figure ), was not detected in either blood or bone marrow (data not shown). Seven days from dosing, although the recovery process had started, the bone marrow cells were still compromised (Figures  and ).

Similar to previous reports,29 we accounted for this observed temporal disconnection between free platinum concentration and effects on haematopoietic cells, modeling carboplatin toxicity as a function of the DNA damage (Damage) caused by the formation of platinum DNA adducts, which persist in the system longer than free platinum28, 29, 30 (Eq. (3); Figure ).

Drug effects in Eqs. (1), (4), and (7) are driven by the damage‐induced cell death, which occurs in the proliferative fraction of the cell population compartment (Eq. (2)). Notice that the cell death rate (Slope) is independent of both the time and the structural scale of the system, and it is ready to be used for the extrapolation of predictions across scales and species with the same sensitivity to carboplatin. Moreover, Eq. (2) describes the drug effect as independent of any feedback mechanism activated by perturbations, which improves the computational tractability of our model when compared with previous approaches.17, 31

Under these assumptions, our model explained the lack of synchronicity between free platinum plasma kinetics (Figures  and ) and the depletion of bone marrow populations previously described.

MPPs are rescued during injury

Similar to previously published models,32, 33 we implemented a flow of cells, or source of HSCs (Source), feeding into the MPP compartment (Eq. (1)). In our model, the value of this flow is regulated by negative feedback signals (fdbkstem in Eq. (1)) coming from the MPP descendants CMPs and/or MEPs, which is supported by experimental evidence.24 A decrease of the descendants below their homeostatic values (CMP0, MEP0) increases the flow of cells into the MPPs compartment. This was modeled by including the following feedback factor into Eq. (1):which enabled to upregulate MPP levels after carboplatin exposure and hence model the observed oscillations of CMP and MEP in response to carboplatin (Figures  , ). Because of the rather small value of γ stem when compared with the other feedback values, the effect of fdbkstem function on MPP dynamics were quite mild (i.e., k stem rate increased 27% when the feedback was at its maximum). However, this was enough to counterbalance MPP cell loss as a result of their differentiation into CMP and MEP,24 which are also regulated by feedback functions (Eq. (12)).

CMPs fail to balance carboplatin‐induced cell loss

We observed a mild depletion of CD45+ cells (CMP); however, the recovery of this population was relatively slow, and homeostasis was not restored before the start of the next 14‐day cycle, leading to a progressive reduction of CMP population from the first carboplatin dose (Figures  and ). This was explained assuming that CMP proliferation rate together with feedback mechanisms controlling the in‐flow of cells into CMP compartment were not able to rebalance the cell loss induced by carboplatin. Consequently, white blood cell counts decreased progressively during the administration schedule (Figures  and ), and we estimated a recovery period of ~50 days to reinstate at least 80% of the original CMP population after the last carboplatin cycle.

Circulating cells regulate MPP cell fate

Analysis of the progenitors in the bone marrow revealed that the most affected populations were CD90+/Lin‐ cells and CD71+ cells, i.e., MPPs and MEPs, respectively. In the first 3 dosing cycles, both MPP and MEP experienced great cell loss followed by a recovery process that, because of the well‐established negative feedback loops,34, 35 led to values first exceeding and then converging to baseline before the following cycle started (Figure  ). In our experiments, where bone marrow status was not assessed again until the eighth cycle, full recovery was not achieved (Figure ). We assumed that this resulted from the interaction of different cellular prosurvival mechanisms that regulate bone marrow cell dynamics,25 and we defined not only proliferation but also maturation25, 36 as dependent on the level of circulating cells through feedback loops. We modeled the cell flows from MPP into CMP and MEP as dependent on feedback signals originated from circulating cells with the following functionswhich are integrated in Eqs. (1) and (4). These functions define cross‐lineage effects and most likely describe the impact of signals mediated by lineage‐specific cytokines,7, 37 such as erythropoietin for erythrocytes, granulocyte colony stimulating factor for neutrophils, and macrophage colony stimulating factor for monocytes, on MPP cell‐fate decision.38 This is a vital process enabling the rapid recovery following myelosuppression; however, we do not have experimental measurements of such cytokines in this particular study, and as a result we have taken an empirical approach based on deviations of mature circulating cell counts from baseline values. Differentiation of MPP into CMP is favored by reduced white blood cell counts39 or equivalently by increased granulocyte and monocyte stimulating factors, whereas maturation toward MEP is controlled by red blood cells29 or erythropoietin. The effects of these feedback functions can be quantified from the recovery times required by each progenitor cell population after the last carboplatin cycle. We estimated that MEP require 12 days to reinstate at least 80% of their original population size vs. the 30 days required by MPP.

Figure  shows that, when considering feedback from circulating cells on MPP maturation, as described in Eq. (12), the model predicted the observed decreasing numbers of progenitors in successive cycles, whereas this pattern was not predicted when feedback mechanisms were modulating only proliferation (Figure ). This hypothesis is supported by independent experimental findings, which suggest that high levels of erythropoietin, in response to the progressive depletion of red blood cells, introduces bias in the differentiation process, guiding MPP fate in favor of committed erythroid progenitors.40

Reticulocyte and platelet transit cells proliferate

Carboplatin induced severe depletion in reticulocytes and platelets, for which the nadir was more than 95% and 45% reduction from baseline, respectively. However, as previously reported,29, 41 within each one of the eight cycles, we observed a strong rebound of these populations exceeding by 290% and 180% their baseline, respectively, before the recovery of homeostatic values. These large oscillations could not be explained by toxicity and proliferation of their common progenitor, MEP (Figure ). We hypothesized that the first two transit compartments of reticulocytes35 and platelets42 are proliferative and thus also sensitive to carboplatin. Furthermore, we also considered that erythropoietin and thrombopoietin significantly affected proliferation in these erythroid and platelet transit‐amplifying compartments by a negative feedback loop from circulating cells expressed by the following functionswhich are included in Eq. (7).

Under the hypothesis that immature transit cells are proliferative (i.e., late progenitor cells), the large oscillations observed in circulating reticulocytes and platelets were accurately described by the model (Figure  ). Using stability analysis, we found that the power of the feedback function defined in Eq. (13) (γprolTrans) is a crucial parameter that can lead our system into oscillatory solutions and that does not describe a physiological haematopoiesis process (Figure ).

Clinical predictions

We assumed that parameter values, but not model structure (i.e., haematopoiesis biology), were species dependent. Model predictions generated scaling our mathematical model (Eqs. (1), (2), (3), (4), (5), (6), (7), (8), (9)) to reflect the differences between rat and human haematopoiesis physiology (e.g., population size, maturation time, etc.), but using the drug‐effect parameters estimated in rats did not reflect accurately either clinical neutropenia or thrombocytopenia (Figure  ). Therefore, we scaled our mathematical model to also account for the differences in sensitivity to carboplatin between human and rat bone marrow cells. Then, we provided predictions in patients and compared these with simulated clinical data, which were generated using the model by Schmitt et al.5 (Supplementary Text ).

Figure 4

Clinical predictions of carboplatin‐induced neutropenia and thrombocytopenia were generated translating our quantitative systems pharmacology model of carboplatin‐induced myelotoxicity from rat to human. Clinical predictions (solid red lines) generated using the same slope values as those estimated from the rat (a) and predictions generated using slope values adjusted for species‐specific drug sensitivity, Eq. (10), Table  (b), were compared with clinical data of carboplatin‐induced neutropenia and thrombocytopenia described in ref.5 which are represented by solid (means) and dotted (percentiles) blue lines. Blue shadow areas show percentiles in our simulations (1,000 individuals), dotted horizontal black lines are thresholds for grades 3 and 4 (neutropenia and thrombocytopenia). Two cycles of 21 days of carboplatin were simulated. Doses were calculated with the Calvert formula,51 targeting an AUC (0–24) of 5 (mg/mL) minutes (Supplementary Text ). The predicted percentages of patients with neutropenia (top) and thrombocytopenia (bottom) grades with those in ref.5 are compared (c). d, days. AUC(0‐24), area under the curve in the 24 h time interval; IC50, half maximal inhibitory concentration.

Our in vitro experimental results suggested that carboplatin is more potent in rat than human, with MEP being more sensitive than CMP in both species (Table ). After scaling the drug‐effect parameters according to the observed species‐specific drug sensitivity (Eq. (10)), both clinical neutropenia and thrombocytopenia predictions improved (Figure  ). Precisely, the percentages of patients experiencing grade 3 and grade 4 neutropenia and grade 3 thrombocytopenia were predicted within onefold difference (Figure  c). However, grade 4 thrombocytopenia remained underpredicted by fourfold, suggesting that further understanding of the differences between human and rat platelet biology is crucial for increasing confidence in the predictions of severe clinical thrombocytopenia. Other factors that could affect the predictions are the underlying disease itself or impaired bone marrow function as a result of previous treatments, which would also impact on haematopoietic recovery, and these features were not accounted in our translation.

Discussion

In this study, we described a new systems pharmacology model of haematopoiesis that quantitatively and mechanistically links carboplatin‐induced effects in the bone marrow with multiple clinical endpoints, considering various lineages of the haematopoietic phylogeny, cellular turnover, and interaction among lineage‐dependent feedback loops (Figure  ) rather than effects on single lineages in isolation.7 We used the rat as a model system to investigate haematotoxicity, and we accurately predicted the temporal dynamics of bone marrow progenitors and mature blood cells with our theoretical model.

Carboplatin‐induced loss of proliferative cells is driven by DNA damage, which is caused by the formation of platinum DNA adducts,28, 30 and these adducts are stable for days28, 30 after carboplatin is cleared. We modeled the temporal dynamics of the damage induced by platinum DNA adducts (Eqs. (2), (3)) in proliferative progenitors, and this was sufficient to capture the observed lack of synchronicity between carboplatin plasma kinetics and the response of bone marrow populations (Figures  and ).

When homeostasis is perturbed, the number and type of progenitor cells in the bone marrow, as well as of circulating cells, influence the balance between cell in flow and out flow in early progenitor populations via complex interacting feedback signals to facilitate the recovery of the system. In our model, the feedback signals affect not only the proliferation potential of susceptible populations but also impact cell‐fate decision at a higher level by modulating the competitive equilibrium between pro‐proliferation and pro‐differentiation inputs coming from other progenitors and mature cells. For instance, we observed that MPP progenitors lost their capacity to recover baseline values after successive dosing cycles (Figure ). To explain this result, and in agreement with previously reported findings,40 we hypothesized that, after several cycles of carboplatin, erythropoietin levels raised to rescue the progressive (but late) detectable loss of red blood cells, thus inducing bias toward erythroid lineages.40 This is achieved by disrupting the equilibrium controlling cell‐fate decision and favoring MPP differentiation into the erythroid lineage over proliferation (Figure  ). These model assumptions reflect a degree of plasticity of haematopoietic progenitors to rapidly respond to strong perturbations.43 In this respect, a number of studies suggest that, in some pathological conditions, HSCs are even able to shift out of the haematopoietic hierarchy and transdifferentiate into cell types that belong to other organ systems, although these results remain controversial.44

Furthermore, we observed large oscillations in reticulocyte and platelet counts, which MEP proliferation activity coupled with carboplatin cytotoxicity were not able to capture (Figure ). We described these large variations (Figure  ) by assuming that T1Ret, T2Ret, T1Plt, T2Plt are proliferative compartments (i.e., late progenitor populations) and thus are susceptible to carboplatin toxicity, and with a proliferation potential responsive to erythropoietin‐mediated35 and thrombopoietin‐mediated42feedback, respectively.

Beyond this exemplar case with carboplatin, our approach easily enables the incorporation of additional motifs capturing cell type–dependent and drug‐dependent mechanisms that act at the molecular scale and drive the haematopoietic response. For instance, modeling how specific progenitor cell types cope with DNA damage45 will enable the quantification of the haematological toxicity of DNA damage repair inhibitors. Similarly, the effects of cytokines, such as granulocyte‐colony stimulating factor on progenitors’ proliferation and differentiation46 can be measured and then explicitly considered to predict the rescuing effects of adjuvant therapies during chemotherapy.39 For this purpose, hypotheses derived from in vitro systems, such as bone marrow organ‐on‐a‐chip, where a niche of HSCs and progenitors sustains a dynamic haematopoietic system,47 could be integrated into the model to quantify cytokine‐driven feedback dynamics.7 Other aspects are the investigations of bone marrow microenvironments or niches6 and the spatial properties of the haematopoietic system. Our model hypotheses can be integrated in computational models48 and/or partial differential equations that, combined with high‐resolution imaging techniques, will enable the study of the effect of the microenvironment or niche6 on HSCs and early progenitors cell‐fate decision.49

We demonstrated the application of our model to describe carboplatin toxicity and, based on the assumption that the same core structure would apply to human (i.e., assuming a conserved haematopoiesis biology across species), we showed the translation of our theoretical insights and predictions for the observed clinical neutropenia and grade 3 thrombocytopenia induced by carboplatin treatment (Figure  ). We accounted for known cross‐species differences in model parameters, tested hypotheses about exposure–response relationships, and made assumptions where knowledge was not available (e.g., allometric scaling of parameters).50 Although extrapolation beyond the observed data must always be undertaken with caution, our multidisciplinary approach can be easily adapted to model the impact of other drugs on patients’ bone marrow, where the sampling of such cells is rarely feasible. Likewise, bone marrow disease models, such as for acute myeloid leukemia,31 can be straightforwardly built on our theoretical framework to provide prospective clinical predictions for different patient populations.

Altogether, our quantitative systems toxicology approach enables the investigation of the effects of compounds and combinations, with substantial potential to explore new schedules and inform drug‐induced effects in man, prior to the initiation of clinical trials.

Funding

C.F. was a fellow of the AstraZeneca postdoc program.

Conflict of Interest

C.F., L.O.O.C., C.P., J.W.T.Y., J.T.M., and T.A.C. are AstraZeneca employees; J.W.T.Y., T.A.C., and L.O.O.C. are shareholders of AstraZeneca. D.I.J. receives funding for clinical trials and associated laboratory studies from AstraZeneca.

Author Contributions

C.F., C.P., T.A.C., J.W.T.Y., J.T.M., D.I.J., and A.C. wrote the manuscript. C.F. and T.A.C. designed the research. C.F., L.O.O.C., and A.S. performed the research. C.F. and L.O.O.C. analyzed the data.

Figure S1. Overall myelosuppression profile induced by one cycle of carboplatin 30 mg/kg.

Click here for additional data file.

Figure S2. Residuals.

Click here for additional data file.

Figure S3. Free platinum pharmacokinetics in rat and damage response.

Click here for additional data file.

Figure S4. Bone marrow status after several cycles of carboplatin is compromised.

Click here for additional data file.

Figure S5. Alternative scenarios for competitive feedback mechanisms.

Click here for additional data file.

Figure S6. Proliferation in immature transit cells result in relatively large oscillations in reticulocytes and platelets as experimentally observed.

Click here for additional data file.

Figure S7. Bifurcation plots and model oscillations.

Click here for additional data file.

Table S1. Designs of rat in vivo studies.

Click here for additional data file.

Table S2. Species‐specific drug sensitivity quantified with colony forming unit assays for both human and rat progenitor cells.

Click here for additional data file.

Table S3. Carboplatin pharmacokinetic model parameter values.

Click here for additional data file.

Table S4. Residual variability values.

Click here for additional data file.

Supplementary text.

Click here for additional data file.

Model code.

Click here for additional data file.

Supplementary References

Click here for additional data file.