Towards Quantitative Systems Pharmacology Models of Chemotherapy-Induced Neutropenia.

Neutropenia is a serious toxic complication of chemotherapeutic treatment. For years, mathematical models have been developed to better predict hematological outcomes during chemotherapy in both the traditional pharmaceutical sciences and mathematical biology disciplines. An increasing number of quantitative systems pharmacology (QSP) models that combine systems approaches, physiology, and pharmacokinetics/pharmacodynamics have been successfully developed. Here, I detail the shift towards QSP efforts, emphasizing the importance of incorporating systems-level physiological considerations in pharmacometrics.

As first‐line defenders playing direct and indirect roles in both innate and adaptive immunity,1 neutrophils—the most abundant white blood cell in the body2—are essential components for the maintenance of overall health. As with all blood cells, neutrophils are produced by the hematopoietic stem cells (HSCs) in the bone marrow. Although there has recently been some question as to the length of time neutrophils remain in the blood,3 mathematical modeling efforts support the historic view that neutrophils are short‐lived in the circulation.4, 5 The neutrophil lineage is therefore particularly susceptible to the cytotoxic effects of chemotherapeutic drugs,6 as it is in a state of constant production to ensure adequate numbers are reserved in case of emergency.7 Neutropenia (depressed neutrophil counts) and febrile neutropenia increase the risk of hospitalization and complications from morbidities as a result of a heightened susceptibility to infections. Due to the prevalence of neutropenia during chemotherapeutic treatments, absolute neutrophil counts (ANCs) and the hematopoietic response are often surrogate metrics for treatment success.8, 9, 10 Accordingly, as a leading cause of dose adaptation or complete treatment cessation,11 myelosuppression (the decreased production of blood cells brought on by the undesired cell‐killing mechanisms of chemotherapeutic agents) is a primary concern for drug developers, clinicians, and patients.

Consequently, the call for mathematical models that are able to predict neutrophil dynamics during anticancer therapy is high. In addition to chemotherapeutic dosing levels, the dose schedule, the number of repeated administrations, and the period between doses can each have significant effects on the development of neutropenia and patient outcomes.10 Many authors have studied blood cell dynamics and hematopoiesis using a variety of mathematical models (an extensive survey is available in Pujo‐Menjouet12). In this review, I will highlight significant modeling advances with respect to neutropenia during chemotherapeutic treatment. The studies discussed here can broadly be divided as having been developed in one of two fields, with advances in one discipline having historically not been adopted in the other. However, using progressively translational modeling approaches and techniques, researchers have increasingly aimed to bridge the divide to provide a more holistic and physiologically detailed understanding of the mechanisms underlying the neutrophil lineage's response to anticancer drugs to better predict treatment outcomes. It is my hope that this review will stimulate and encourage the further blending of traditional pharmacokinetic/pharmacodynamic (PK/PD) approaches and applied mathematical methodologies to better respond to clinical and patient needs.

This article is divided as follows: I will begin by providing an overview of the mechanisms and behaviors of the neutrophil lineage from HSCs to commitment, release, and subsequent removal through apoptosis or margination from the circulation. Next, I highlight the major neutropenia model used in the field of pharmaceutical sciences. I will then detail the multiple modeling approaches developed by mathematical biologists before outlining quantitative systems pharmacology (QSP) models of chemotherapy‐induced neutropenia. I will conclude with a more general discussion on the role QSP approaches play in the future of hematological and pharmaceutical modeling as a whole.

NEUTROPHIL PRODUCTION

Today's highly developed stem cell research field was born of the discovery of these highly proliferative cells by Till and McCulloch at the University of Toronto after they injected bone marrow into mice spleens and observed clonal nodules in the excised organs.13 The hematopoietic stem cells (HSCs), a particularly well‐studied stem cell subpopulation, regulate all of hematopoiesis (the production of blood cells) and their proliferative capacity is almost stupefying: HSCs produce nearly 10 times a human's body weight in blood cells per lifetime,14 yet comprise just 0.01–0.2% of the total bone marrow mononuclear cells.15 However, much remains unknown about the mechanisms instructing HSC lineage potential,16 and the correct paradigm of lineage determination remains an open question.17

In the traditional understanding of HSC commitment, HSCs begin by ramping up mitotic events as multipotent progenitors (MPPs) before differentiating into one of two progenitor types, the common myeloid progenitors (CMPs), or the common lymphoid progenitors (CLPs). After commitment into the myeloid lineage, HSCs begin the process of becoming one of a number of terminally differentiated cells, including neutrophils, whose development occurs only in the bone marrow, where they proliferate and mature in a process known as granulopoiesis.18 The formation of neutrophils from CMPs involves several transitional steps. First, CMPs differentiate into myeloblasts, which are the common progenitor of all granulocytes (basophils, eosinophils, and neutrophils) and monocytes. Myeloblasts then further differentiate to become promyelocytes, myelocytes, metamyelocytes, and band neutrophils before they proliferate and mature into neutrophils (see figure 1.1 in Ref. 19 and figure 2 in Ref. 20).21 These mature neutrophils are subsequently stored in the bone marrow prior to their release into the circulation.7

As previously mentioned, the historic perspective is that neutrophils have a short half‐life in the circulation of 7–10 h.2 A recent report using heavy water deuterium suggested that neutrophils have longer half‐lives of 3.7 days,3 a finding that would have significant implications on our understanding of neutrophil kinetics.5 A confounding factor in the determination of neutrophil fates in the circulation is the presence of a large marginated pool (marginated neutrophil pool, or MNP) that exists in several tissues and organs in addition to the capillaries and renders the quantification of circulating half‐lives difficult, as neutrophils leave the circulation not just through cell death.4 The existence of the MNP is experimentally verifiable: during neutrophil reinfusion, one‐third of reintroduced cells were found in the liver and the bone marrow, while ∼15% migrated to the spleen.2 It is further known that the lungs harbor a significant amount of the total blood granulocyte pool. The fluctuating relationship between the freely circulating neutrophils and those lodged in the tissues and organs, together with experimental limitations and complications presented by labeling techniques,22 may account for the discrepancy between historical reports and more recent findings. However, two independent mathematical models have reinforced the previously held belief of short half‐removal times.4, 5

After margination, neutrophil clearance is carried out by both the spleen and the liver, although some suggest that neutrophils that have trafficked back into bone marrow may be cleared from there.7, 23, 24, 25 Neutrophil removal is accomplished by macrophages at the various sites of clearance,7 during which the affected neutrophils release interleukin (IL)‐23, setting off a cascade of cytokine secretion inducing an increase in granulocyte colony‐stimulating factor (G‐CSF) concentrations.2

G‐CSF is considered the primary cytokine driving granulopoiesis,7 although a broad set of signaling molecules are involved in regulating neutrophil counts, including IL‐3, granulocyte‐macrophage colony stimulating factor (GM‐CSF), and IL‐6. When IL‐3, GM‐CSF, or IL‐6 receptors are knocked out in mice, however, no defects are observed in the production of granulocytes. In contrast, in G‐CSF(‐/‐) mice, only 20–30% of normal neutrophil counts are present and the mobilization of neutrophils from the marrow to the circulation is impaired. Absolute neutropenia in humans has also been observed when deficiencies in G‐CSF receptors (G‐CSFRs) are present. In one case, a child born with severe chronic neutropenia (no detectable neutrophils) was found to be unresponsive to treatment with G‐CSF as a result of a point mutation on the extracellular domain of their G‐CSFRs, suggesting that G‐CSF is indispensable for neutrophil regulation.26

G‐CSF is a potent regulator of both neutrophil production (by various means) and of the release of neutrophils from the marrow into circulation,7, 27 but it does not seem to affect the clearance or removal of neutrophils.18 The most immediate effect of an increase in G‐CSF concentrations is emergency granulopoiesis, or the rapid mobilization of neutrophils into the blood.7, 28 Unlike homeostasis, neutrophil transit from the marrow seems to also be mediated by IL‐17 during inflammatory events.7

To compensate for the release of cells from the marrow reservoir, G‐CSF acts upstream on proliferating and maturing cells to ensure the replenishment of the reservoir by increasing cell proliferation, decreasing the speed of transit from the marrow, and ensuring the differentiation and maturation of neutrophils.19 To that end, it has been observed that neutrophil precursors are subject to higher rates of apoptosis (programmed cell death) in the absence of G‐CSF.29 An overview of the various actions of G‐CSF is provided in Figure 1.

Figure 1

An overview of granulopoiesis. As with all blood cells, neutrophils begin as hematopoietic stem cells (HSCs, orange circle) in the bone marrow (pale yellow background), where they develop. HSCs are capable of self‐renewal and are subject to cell death (dashed arrows). HSCs may also differentiate into one of the blood cell lines, including the neutrophils (purple circles). After commitment to the neutrophil lineage, cells undergo a period of proliferative expansion at the end of which they no longer divide. Postmitotic neutrophils then mature, growing in size and gaining receptors. At the end of the maturation process, cells are then stored in the bone marrow reservoir from which they egress to reach the circulation (pale red background) before removal (by margination or death). G‐CSF acts to modulate the rate of exit from the marrow reservoir, increase the rates of maturation and proliferation, and to modulate the rate of differentiation into the neutrophil lineage (G‐CSF actions represented by blue vertical arrows).

The initiation of G‐CSF's action on neutrophils is produced through binding to G‐CSFRs on the surface of the cells that subsequently degrade the G‐CSF and remove it from the blood. G‐CSF is also cleared by the kidneys,19, 29 through a linear clearance pathway activated primarily after the neutrophil‐binding removal mechanism becomes saturated. These dual elimination routes are important drivers of the PK/PD of G‐CSF and they have a notable effect on the half‐life of G‐CSF in the blood: G‐CSF has a half‐life of around 4.7 h without neutrophils available for binding, but just 2 h when they are available.29 The cooperation between the two clearance pathways is particularly relevant to the PKs of G‐CSF, and some modeling work may have previously overestimated the proportion of renal removal vs. that attributable to the neutrophils.30

CHEMOTHERAPY‐INDUCED NEUTROPENIA

Anticancer drugs disrupt the growth of cells by disturbing DNA synthesis and/or inhibiting a cell's ability to properly cleave by interfering with microtubule elongation and contraction. Unfortunately, the nature and efficacy of this cell‐killing mechanism leads to cytotoxic secondary effects.31, 35 At homeostasis, neutrophil precursors in the bone marrow must be steadily dividing to maintain healthy blood neutrophil counts. Circulating neutrophil numbers can therefore become significantly reduced during anticancer treatment as they are subjected to the antimitotic effects of chemotherapy drugs, effects that rarely impact slowly dividing cells like the HSCs. The increased death of early progenitor cells is further aggravated by the repeated exposure to cytotoxic drugs inherent to current periodic combination chemotherapy regimens.10 Recognizing the likelihood of cytotoxicity during chemotherapy, hematopoietic rescue drugs like exogenous G‐CSF can be coadministered (prophylactically or as adjuvants). Indeed, due to the decreased incidences of febrile neutropenia, a reduced risk of infections has been observed in certain cancer patients treated concomitantly with exogenous G‐CSF.32 Further, given that G‐CSF can rescue falling neutrophils counts, its coadministration during chemotherapeutic regimens have been shown to increase survival times due to the ability to increase the intensity of the chemotherapeutic treatment.33, 34

In 1991, the US Food and Drug Administration (FDA) approved the use of G‐CSF mimetics during chemotherapy to address the prevalence of neutropenia.11 To appropriately triage patients for effective care during chemotherapeutic treatments, regulatory agencies have since put forward methods for standardizing clinical decisions on the use of exogenous G‐CSF (typically in biosimilar forms like filgrastim but increasingly in its PEGylated form, pegfilgrastim, which bypasses renal clearance routes to maintain higher circulating drug concentrations for longer) in oncological settings. Differences nonetheless remain between recommendations from different organizations owing to the prevalent practices of their respective regions35 and despite these guidelines, the use of G‐CSF during myelosuppressive chemotherapy has been observed to be nonstandard in practice.36

It is therefore evident that mathematical modeling is well positioned to respond to the need to forecast the neutrophil response during chemotherapeutic administration. Model‐based predictions can help to standardize treatment protocols and take out some clinical guesswork, leading to increasingly evidence‐based practices.10 To that end, conventional PK/PD and pharmacometric approaches continue to hold an undeniably important place along the drug discovery pipeline, from early phase trials to after‐market regimen evaluation and optimization,37 and are heavily relied upon in industry to evaluate candidate drugs and their potential efficacy.38 Examples of PK/PD (including population PK/PD studies) used in drug development are therefore numerous (see Refs. 37 and 38, and references therein, for example). The field of neutropenia modeling is no exception, and the construction of models with consistency across parameters (i.e., incorporating quantifiable biological mechanisms) has been a focus of both the pharmacometrics and mathematical biology communities for decades. One researcher in particular, Lena Friberg of Uppsala University, Sweden, has had perhaps the broadest impact on modeling chemotherapy‐induced neutropenia within pharmacometrics, as discussed in the following section. A summary table of the models discussed in future sections is given in Table 1.

Table 1

Summary of discussed models by discipline and type

	Author(s)	Focus	Citation
	Friberg	Chemotherapy‐induced myelosuppression	39
	Venkatakrishnan et al.	Optimization of oncological therapeutics	44
	Bulitta et al.	Paclitaxel PopPD and neutropenia	45
	Soto et al.	Case study: semi‐mechanistic model + novel chemo drug	46
Semimechanistic models	Fetterly et al.	Paclitaxel PD model	47
	Kloft et al.	Identification of patient subgroups across chemo drugs	48
	Léger et al.	Topotecan‐induced neutropenia	49
	Jayachandran et al.	Optimal chemo regimens for leukemia	50
	Quartino et al.	Endogenous G‐CSF and relationship to myelosuppression	41
	Krzyzanski et al.	PopPKPD filgrastim (no chemo)	55
TMDD	Wang et al.	PopPKPD filgrastim (no chemo)	56
	Pastor et al.	Model of G‐CSF effects during carboplatin treatment	57
	Mangas‐Sanjuan et al.	Semimechanistic cell‐cycle model for diflomotecan schedules	52
Extensions	Parker	PBPK modeling + chemotherapy design	58
	Ho, Clermont, & Parker	Neutrophil response during inflammation + cancer	59
Stochastic	Krinner et al.	ODE granulopoiesis model + stochasticity	60
	Steinbach et al.	Canine granulopoiesis	65
	Loeffler & Wichmann	Stem cell proliferation model	66
	Wichmann, Loeffler, & Schmitz	Hematopoietic regulation model	70
ODEs	Engel, Scholz, & Loeffler	Myelosuppression during multicycle combination chemo	72
	Scholz, Engel, & Loeffler	Design of chemo regimens for leucopenia	73
	Hu et al.	Hemodose (unintented radiation model)	74
	Graessle & Fliedner	Model for severity of radiation exposure on cell renewal	75
	Rubinow & Lebowitz	Neutrophil production/control	83
	Schwegler & Mackey	Effects of noise on cell‐cycle model after chemo	84
DDEs	Zhuge, Lei, & Mackey	Neutrophil dynamics after periodic chemo	85
	Lei & Mackey	Review of age‐structured models + treatment of cytopenia	86
	Brooks et al.	Neutrophil response during chemo + G‐CSF	87
	Mouser et al.	Hematopoietic response during chemo + stimulant	89
	Scholz, Engel, & Loeffler	Model of combination chemo + G‐CSF support	71
	Vainstein et al.	Physiological granulopoiesis model and G‐CSF effects	91
	Vainas et al.	Personalizing docetaxel with G‐CSF regimens	92
	Craig et al.	G‐CSF regimen optimization	88
	Craig et al.	G‐CSF/chemo IIV effects on physiological model	93
	Craig, Humphries, & Mackey	Determining G‐CSF PK from physiological model	30

Pharmaceutical sciences: light gray background; mathematical biology: white background; quantitative systems pharmacology: dark gray background. Articles in bold are of particular interest. TMDD, target‐mediated drug disposition; ODEs, ordinary differential equations; DDEs, delay differential equations.

NEUTROPENIA MODELS: DIFFERENT APPROACHES TO THE SAME QUESTION

The Friberg semimechanistic model

From its introduction, the semimechanistic model of Lena E. Friberg39 has undoubtedly become the most frequent and gold‐standard approach to modeling the myelosuppressive effects of chemotherapy in the pharmaceutical sciences literature. As summarized in Figure 2 , the model describes the phases of neutrophil production through the incorporation of key mechanisms along the developmental pathway from the early neutrophil progenitors to circulating cells. To account for the continuous mitosis of the proliferative phase, neutrophils remain in the proliferative compartment through a self‐renewal process. The path from proliferation to the circulation occurs after cells transit through several artificial compartments that mimic the maturation delay naturally present during neutrophil development. Circulating cells are then lost from the blood through linear cell death, and the size of the circulatory compartment exerts negative feedback on the rate of self‐renewal of the proliferative cells (e.g., if circulating cell concentrations are high, there is less self‐renewal, and vice versa). The myelosuppressive actions of chemotherapeutic drugs affect the rate of self‐renewal of the proliferative cells and thus induce a loss from the proliferative compartment. Friberg originally presented a myelosuppressive model with three transit compartments accounting for the delay between proliferation and circulation. However, since the publication of the Friberg semimechanistic model,39 several authors have extended the model structure beyond the three original transit compartments to n such transitional phases.40, 41, 42 As such, a more general model of n transit compartments is given by the systems of equations in (1): where P is the proliferative cells, k is the rate of proliferation, E is the PD model of the myelosuppressive action of the drug, N 0 is the baseline neutrophil count, N is the current number of circulating cells, γ modulates the effect of the feedback from the circulating compartment to the proliferative compartment, k is the rate of transit, T (i = 1,…,n) is the i transit compartment, and k is the rate of exit from the circulating compartment.

Figure 2

Schematic representation of the general semimechanistic model of neutrophil development developed by Friberg.39 Proliferative cells self‐renew at rate k or begin the transition to the circulation by exiting with rate k (the rate of transit). The delay between the time cells leave proliferation to when they enter the circulation is called the mean transit time and is equally divided between n transit compartment, each connected by the rate of transit. The original model as presented in Friberg39 and used most frequently is restricted to three transit compartments. More recent extensions have more recently appeared.40, 41 Once cells enter the circulation, they are removed with rate k. The number of circulating cells has a negative feedback on the proliferation rate of the proliferative cells. The myelosuppressive action of the drug is assumed to also affect k.

Data from drug and clinical trials are most often limited to plasma drug concentrations and ANCs, constraining the number of parameters that can be estimated directly from clinical results using available statistical techniques. To account for this limitation, predicting the ANC response following chemotherapy using the model in Friberg39 can be carried out using just the model (at an individual or population level43), the dosing history, a baseline ANC measurement, and another neutrophil datapoint around the nadir.44 The Friberg model is therefore straightforward to implement and provides good fits to clinical data; thus, it has since been adapted by many researchers seeking to predict the cytotoxic effects of chemotherapy (for a decidedly nonexhaustive sample, see, for example, Refs. 45–50).

Owing to the lack of clinical data characterizing the upstream processes of granulopoiesis, some simplifications, like a constant transit time between maturation compartments, for example, must be made within the framework presented in Friberg.39 Further, the reliance on data to dictate the model's structure can lead to disparate estimates for otherwise well‐characterized physiological mechanisms (such as the mean maturation time from proliferation to the circulation)30 and restrict the predictive ability of such empirically derived methods.51

To address these limitations, some authors have extended the application of the model in Eq. 1 to problems beyond myelosuppression by combining the Friberg model with increasingly complex mechanisms and relationships. In a study of the neutropenic effects of diflometcan, the cell‐cycle dynamics of the stem cells (quiescent and proliferative) were adapted within the Friberg model framework.52 The estimation of the length of the cell‐cycle was shown to agree with experimental data on the length of HSC division time, and the cell‐cycle neutropenia model accounted for differences between different dosing schedules. Quartino et al. accounted for the PKs of endogenous G‐CSF using stimulatory models for the drug effects in lieu of the myelosuppressive action of chemotherapy41 (model schematic given in Figure 3). The study of the G‐CSF effects during exogenous treatment is a natural offshoot of neutropenia modeling, and various groups have modeled these PK/PD relationships using approaches similar to the Friberg model.39 When the PKs of G‐CSF are the modeling focus, target‐mediated drug disposition (TMDD) models, pioneered through the work of Ref. 53 (tracing back to Ref. 54), are also frequently used. TMDD models have mainly been applied to model the PKs of exogenous G‐CSF without reference to chemotherapeutic scenarios,55, 56 but have also been incorporated to study the coadministration of G‐CSF during anticancer treatment.57 The group of Robert Parker at the University of Pittsburgh and associated authors have also looked extensively at the neutrophil response in chemotherapeutic applications and beyond through the incorporation of physiologically based pharmacokinetic (PBPK) models58 and more complex immunological relationships.59 In Ref. 59, the actions of a variety of cytokines and cell populations are included, and the model captures both fast (inflammatory) and slow (cytotoxic) hematological dynamics.

Figure 3

The integrated G‐CSF‐myelosuppression model describing the dynamics of endogenous G‐CSF and neutrophils following chemotherapy. For the myelosuppression model the parameters are baseline neutrophil count (ANC0), mean maturation time (MMT = 5/ktr), the half‐life of neutrophils in circulation (t1/2circ = LN(2)/kcirc), the feedback parameters of G‐CSF on neutrophil proliferation (γ) and transit time (β) and the drug‐related effect (Edrug). The estimated parameters for the G‐CSF turnover model are baseline G‐CSF (GCSF0), nonspecific elimination rate constant (ke), and ANC‐dependent elimination rate constant (kANC) and cortisol‐induced G‐CSF release (DOSEcort) and the half‐life of cortisol‐induced G‐CSF release (t1/2cort = LN(2)/kcort). Figure reproduced from “Characterization of endogenous G‐CSF and the inverse correlation to chemotherapy‐induced neutropenia in patients with breast cancer using population modeling,” Pharm. Res. 31, 2014, p. 3396, Quartino, A.L. et al.,41 with the permission of Springer.

Modeling perspectives from the mathematical biology community

The models highlighted in the previous section were constructed to respond to very defined clinical and drug development needs, chiefly to explain observed neutrophil responses during chemotherapy treatment. Generally speaking, mathematicians tend to seek out models that are capable of explaining many situations at once,60 and are further interested in investigating dynamic presentations of hematopoietic diseases, including the sustained oscillations in neutrophil counts present in patients with cyclical neutropenia. Cyclical hematopoietic diseases were addressed as early as the late 1960s by many.61 The group of Alec Morley, in particular, was interested by cyclical neutropenia,62 a fascination later picked up by Michael Mackey,63 who would become a pivotal figure in the hematological modeling field.10 Broadly, models arising in mathematical biology use systems approaches64; in contrast to the well‐developed and heavily used methods present in the pharmaceutical sciences, parameter estimation using first principles is undertaken in the applied mathematics community with less emphasis placed on statistical parameter fitting techniques. However, owing to a lack of exposure to the PK/PD modeling field,10 relevant drug disposition and effects models have been largely absent from approaches in mathematical biology until more recently. Nonetheless, the incorporation of PK/PD models is seen throughout the work of the group of Markus Loeffler and associated members, as discussed in greater detail in the following section.

Ordinary differential equation neutropenia models in mathematical biology

As a modeling framework, ordinary differential equations (ODEs) have several significant advantages. The theory of ODEs is well‐developed, allowing for careful mathematical study (including local and global stability, and bifurcation analyses) and the numerical implementation of ODEs is equally well‐characterized with many computing languages and software packages, both free and proprietary, available to perform simulations. Further, most experimentalists and clinicians have been exposed to ODEs in some capacity and are therefore familiar with the assumptions underlying their use (a fact which is particularly relevant when one is interested in translational modeling).

As we have seen, ODE models have been successfully developed to study neutropenia during chemotherapy within the pharmacometrics community.39, 53, 58 In mathematical biology, the 1980 model of Steinbach et al. that included seven cellular compartments and granulocyte releasing factor (what would later become known as G‐CSF) also relied on ODEs to model canine granulopoiesis.65 There, as in the model of Friberg,39 the subdivision of the maturation process into several transit compartments had no direct link to the physiology and, as a rapid equilibrium between the circulation and the pool was assumed, no marginal pool was incorporated. In contrast to Friberg,39 the mean division time of the early neutrophil progenitors, several transit times, and other model parameters were identified from existing experimental work. Concurrent to the development of the model in Ref. 65, the group of Loeffler and Wichmann began modeling various aspects of hematopoiesis in earnest66, 67, 68, 69 (for insight into the mechanisms and considerations of interest to the group at the time, see Ref. 70).

In the early 2000s, Loeffler and collaborators shifted their focus to the neutrophil response to combination chemotherapy71, 72 and the use of G‐CSF during anticancer treatment.71, 73 A detailed overview of the granulopoiesis model central to these studies is given in Ref. 71, where the general equations for compartment size changes are described by:

although the above summary equation certainly does not convey the sophistication underlying the model's construction and terms. Both G‐CSF and GM‐CSF are included in the PK/PD model as regulators of neutrophil production. The effect of chemotherapy on the neutrophil lineage is exerted on the stem cells, the granulopoietic progenitor cells, and the proliferating and maturing granulopoietic precursor cells, and independence between chemotherapeutic agents is assumed. The model successfully recreates clinical data from a large randomized trial of polychemotherapy protocols run by the German Hodgkin's Lymphoma Study Group and the German High‐Grade Non‐Hodgkin's Lymphoma Study Group.72 More recently, Loeffler, Scholz, and coauthors have incorporated stochastic elements into their modeling approach (see the following section).

(As a brief aside, although the subject of interest here is myelosuppression during chemotherapy, exposure to radiation—as an anticancer treatment or unintentionally—can also cause significant damage to circulating blood counts, in particular to ANCs. There has been some focus on the development of computerized diagnostic assessment tools for relating the significance of accidental radiation exposure using blood cell counts that may be of interest to readers. Both Refs. 74 and 75 developed ODE compartment models to represent the production of granulocytes (and also platelets and lymphocytes) to assess the level of exposure to radiation based on hematopoietic responses. The general granulopoiesis modeling structure in each is similar to that of Eq. (2) and the models in both Refs. 74 and 75 are shown to successfully predict patient data from a variety of accident sites.)

Stochastic elements in neutropenia models

The prevailing view of the development of various tumorigenic and blood cancers is that they are derived from cancer stem cells that act similarly to HSCs to maintain cancerous cell populations with increased proliferative characteristics.76 In cancers like chronic myeloid leukemia, mutations in the HSCs and/or the early progenitor cells lead to the uninhibited proliferation of white blood cells.76 The appearance of such aberrations is random in nature, and the path from a somatic mutation to the development of cancer remains elusive.77 On account of this unknown mutational path and owing to open questions remaining about the lineage specification and commitment of HSCs, and therefore of hematopoiesis as a whole,78 the incorporation of stochasticity into models of hematopoiesis is an avenue of active investigation.

Most authors who have studied hematopoietic models with random dynamics have done so with the focus being the development of cancers78, 79, 80, 81 (with detailed review in Ref. 82). However, as previously mentioned, the Loeffler group has published a hybrid deterministic/stochastic model incorporating the effects of G‐CSF, GM‐CSF, and chemotherapy to incorporate the cell cycle‐dependent effects of chemotherapy.60 The deterministic portion was taken from Ref. 71 (and previously briefly described in the section Ordinary Differential Equation Neutropenia Models in Mathematical Biology), and stochastic effects were incorporated into the HSC compartment owing to the open questions surrounding the linear stem cell lineage determination model (see, for example, Refs. 16 and 17). The model is shown to be in good agreement compared to data from the CHOP‐21 and CHOEP‐21 trials, where periodic combination chemotherapy was supported by adjuvant G‐CSF support. As experimental insight into HSC dynamics progresses, the inclusion of stochastic elements to models of hematopoiesis and myelosuppression, as in Ref. 60, may become increasingly important.

Age‐structured and delay differential models

Perhaps the principal approach to modeling hematopoiesis in mathematical biology is through age‐structured partial differentiation equations (PDEs) used to derive delay differential equation (DDE) models. Owing to the delays inherent to the development of various blood cells (time for proliferation, time for maturation, time in the circulation, etc.), a natural mathematical representation of the process of blood cell production is via DDEs, where the current state of the system depends on the present and the past. In its most general form, a DDE is given by: where x(t) is the solution at the current time t, and x(t – τ) (i = 1,…,n) is the solution at some past time T , which may be a fixed constant, a probability density, a time‐dependent (τ (t)), or a state‐dependent (τ (t,x(t))) delay. Hematological DDE models are usually derived from age‐structured PDEs that track the rate of change of cellular populations with respect to time (as in the ODE case) and age. Then, by the method of characteristics, the age‐structured PDEs are transformed into DDEs by the general method described in Ref. 31.

Rubinow and Lebowitz, who incorporated many of the elements central to more recent models including the G0 cell‐cycle, and the bone marrow reservoir,83 published one of the earliest age‐structured model of granulopoiesis. Their model predates the discovery of G‐CSF, so no explicit growth factor effects were included. The modeling approach of Rubinow and Lebowitz was subsequently picked up by others in mathematical biology, including Mackey at McGill University, who was interested in nonlinear dynamics and dynamic hematological diseases like cyclic neutropenia (see the reviews of Refs. 10 and 31). In 1994, Schwegler and Mackey published an age‐structured PDE model with stochastic elements (noise) in both the proliferation and maturation potentials.84 Later the Mackey group began focusing on myelosuppression and the use of exogenous G‐CSF to stimulate neutrophil production during chemotherapy with the aim of optimizing the timing of G‐CSF support after periodic chemotherapy.

Foley and Mackey, who had previously studied the timing of G‐CSF administrations for patients with cyclic neutropenia, showed that delaying the initiation of G‐CSF treatment was beneficial to the neutropenic response during chemotherapy.6 Other authors in the Mackey group opted to study more fundamental elements of the dynamics of neutrophils during anticancer treatment.85, 86 However, in each of Refs. 6, 85, and 86, simple processes were used to mimic the PK models of each drug without reference to common PK parameters. The group began to shift towards incorporating more detailed PK models in Ref. 87, a move that would then see the inclusion of more traditional pharmacometric approaches in Craig et al.88 and beyond (see the following section).

Other authors have also used age‐structured approaches to model myelosuppression and its treatment with exogenous cytokines. Mouser et al.89 used a DDE model of granulopoiesis to predict HSC and neutrophil responses to chemotherapy with hematopoietic induction agents, like G‐CSF and erythropoietin (the cytokine responsible for stimulating the production of red blood cells). The model includes two compartments, the proliferating cells and the nonproliferating cells, and parameter estimation was performed from the existing literature. Perhaps somewhat contrary to some current clinical approaches, their results indicate that the coadministration of chemotherapy with a hematopoietic‐inducing agent can worsen the leukopenic nadir. That said, this neutropenic exacerbation is supported by other modeling efforts, detailed in the next section. Uncovering such previously unrevealed mechanisms highlights the unique avenue offered by mathematical modeling derived from physiological principles, where novel and unexplored relationships can be brought to light less invasively and without requiring a heavy experimental burden.

BRIDGING THE GAP: MECHANISM‐DRIVEN MODELING APPROACHES INCORPORATING DETAILED PK/PD RELATIONSHIPS

For the most part, due to disparities in general research aims, the well entrenched PK/PD modeling strategies of the pharmaceutical sciences and the first principles approaches of mathematical biologists have not coincided in neutropenia modeling efforts. Indeed, the issue of the separation between more traditional PK/PD modeling groups and systems biology/mathematical biology is not specific to groups focused on neutropenia. In 2011, the United States National Institutes of Health released a White Paper following two workshops focused on QSP and its role in shaping the drug development pipeline and translational medicine.90 The participants and authors concluded that industrial/academic partnerships centered on interdisciplinary approaches should be increasingly emphasized. With specific regard to toxicity and undesired reactions, such as myelosuppression during chemotherapy, the QSP workshop group highlighted the necessity of the simultaneous incorporation of physiological mechanisms, drug disposition, and drug effects to respond to drug safety considerations.90

Within the chemotherapy‐induced neutropenia context, researchers emphasizing QSP approaches (without necessarily deeming their work as such) have approached the modeling of the hematopoietic and granulopoietic systems with vastly different models, yet have nonetheless arrived at fairly consistent conclusions. The remarkable similarity of their results highlights the role physiology plays in determining the fate of xenobiotics and drug responses, further supporting the call for the increased integration of systems approaches to pharmaceutical problems.

Groups merging traditional PK modeling with physiological, first‐principles models have largely focused on the model‐based optimization of the use of G‐CSF during periodic chemotherapy. Scholz et al. studied the impact of heterogeneous risk factors on various administration scenarios, including cyclophosphamide dose escalation, additional etoposide administrations, and cycle‐time reductions, within the CHOP/CHOEP‐14 protocols (14‐day periodic combination chemotherapy with G‐CSF support (filgrastim) from days 4–13).73 In Ref. 71, the timing of G‐CSF support during the CHOP‐14 cycle was investigated. The model, based on the one developed in Ref. 72, includes equations for the compartments essential to neutrophil production (including the HSCs, the proliferative neutrophils, the maturing neutrophils, and circulating cells) while also incorporating chemotherapy and G‐CSF PK models, together with their respective pharmacodynamic effects. The conclusions in Ref. 71 are consistent with those of the other groups highlighted later in this section: too many G‐CSF administrations per chemotherapy cycle are detrimental to the neutrophil response (“the rule ‘much helps much’ is not true in the context of G‐CSF application”71). Scholz et al. observed that too much G‐CSF undermines the desired outcome of G‐CSF support during chemotherapy by causing too many neutrophils to release from the reservoir, thereby overstimulating the upstream processes involved in neutrophil production. They further note that more desirable neutrophil outcomes are possible by delaying the G‐CSF support to day 5 postchemotherapy (G‐CSF administered daily from days 5–10), even with a shortening of the chemotherapy period from 14 days to 11 days. These results, along with the others in this section that call for delaying the first day of G‐CSF administration after periodic chemotherapy, do, however, conflict with current clinical practices based on some evidence of an association between delayed G‐CSF treatment and more febrile neutropenic days/longer neutrophil recovery times.36

The group of Zvia Agur also developed a physiological model of granulopoiesis together with PK/PD models to docetaxel/G‐CSF regimen optimization.91, 92 Figure 1 in Ref. 92 provides an overview of the PK and physiological models' components. Akin to the DDE methodology expanded upon in the section Age‐Structured and Delay Differential Models (above), the model in Ref. 91 incorporated gamma‐distributed density transition times between physiological compartments modeled by ODEs to incorporate the accelerated maturation that occurs with high concentrations of G‐CSF. Age‐structured equations are also used to model the cell‐cycles in the proliferative phases. Although some parameters are fit from available data, many are fixed from literature sources or through model‐based calculations. To personalize schedules, interindividual variability (IIV) in the PKs of docetaxel, related to the effects of CYP3A‐mediated docetaxel clearance, was included in the population PK/PD model.92 After validation, where model predictions are shown to agree with clinical data in a variety of docetaxel schedules, Vainas et al. focused on optimizing treatment schedules with regard to three considerations: the first day of G‐CSF treatment postchemotherapy, the number of G‐CSF administrations given per chemotherapy cycle, and the total daily dose of G‐CSF.92 Their results suggest that the early administration of G‐CSF (1 day post‐docetaxel) depletes the neutrophil reservoir in the bone marrow too soon to compensate for the delayed cytotoxic effects of the chemotherapeutic drug. As compared to the standard docetaxel/G‐CSF regimen, where 300 μg of G‐CSF is administered daily beginning 1 day postchemotherapy, the group's model predicts that a more optimal G‐CSF protocol is to administer four daily administrations beginning 7 days post‐docetaxel, again supporting the idea that delaying the beginning of G‐CSF support is beneficial to the overall neutropenic response.

We also undertook the optimization of filgrastim protocols during chemotherapy in Craig et al.,88 where we incorporated a DDE physiological model of neutrophil production with G‐CSF and chemotherapy disposition and effects models for a typical patient. The schematic of the most recent physiological model30 is provided in Figure 4 .

Figure 4

Schematic representation of the production of circulating neutrophils in the bone marrow and the interaction of the system with G‐CSF. Hematopoietic stem cells (HSCs‐Q) enter the neutrophil lineage, the other blood lines, or are removed from the HSC pool. Differentiated HSCs undergo successive divisions during the proliferative phase. Cells then mature before being stored in the marrow reservoir, or dying off during maturation. Neutrophils remain in the reservoir until they are removed randomly or enter the circulation, where they disappear rapidly from the blood. Freely circulating G‐CSF may bind to receptors on the neutrophils. The concentration of bound G‐CSF drives its pharmacodynamic effects. The concentration of G‐CSF bound to mature neutrophils, G 2, determines the rate of release from the marrow reservoir. The concentration of G‐CSF bound to neutrophil precursors, assumed proportional to G 1, the concentration of freely circulating G‐CSF, determines the rate of differentiation from the HSCs, the speed of maturation, and the rate of proliferation. For all four effects, speed and rates increase with increasing G‐CSF concentration. Figure reproduced from “A mathematical model of granulopoiesis incorporating the negative feedback dynamics and kinetics of G‐CSF/neutrophil binding and internalization,” Bull. Math. Biol., 78, 2016, p. 2308, Craig, M., Humphries, A.R., and Mackey, M.C.30 with the permission of Springer.

The work in Craig et al.88 builds upon previous studies carried in the Mackey's group, including Refs. 6 and 87, by updating the physiological model and combining it with previously published PK/PD models constructed using nonlinear mixed effects modeling techniques. All parameter values were fixed via values available in the literature or through equilibrium relationships. Regimen optimization of the CHOP‐14 protocol examined in Refs. 71 and 73 was undertaken with respect to the day of the first G‐CSF administration after chemotherapy, the number of administrations of filgrastim per chemotherapy cycle, and the time between G‐CSF administrations. Our results again suggest that delaying the administration of G‐CSF postchemotherapy, and then administering it daily, is optimal for reducing neutropenia. Instead of administering filgrastim from days 4 through 13 during a 14‐day chemotherapy period, we found that delaying the first G‐CSF administration to day 6 or day 7 required the administration of as few as four or even three subcutaneous injections to obtain improved neutrophil responses.88 As in Refs. 71 and 92, the basis of retarding the beginning of the G‐CSF support is attributable to the significant depletion of the bone marrow reservoir.30 In a later analysis, we included IIV and interoccasion variability in the PK models of the chemotherapeutic and G‐CSF drugs used in Craig et al.88 were shown to not affect the model's predictions,93 suggesting that the first‐principle physiological model inherently incorporates variability by the nature of its construction.

DISCUSSION

It is increasingly recognized that models in the pharmaceutical sciences must begin to incorporate the mechanisms underlying drug disposition and drug responses to better respond to drug development needs and to broader translational goals.94, 95 The traditional PK/PD modeling framework is firmly rooted in the paradigm of descriptive models that characterize data through the quantification of mean behavior and variability. Often, however, these models can be limited in their ability to extrapolate findings to different pathologies and treatment scenarios.95 Conversely, mathematical and systems biologists are largely focused on modeling mechanisms without reference to specific drug models. A better understanding of how and why drugs work in mechanistic and systemic ways contributes to drug discovery targets and to after‐market objectives, like the minimization of toxicity. However, it should be recognized that systems‐level QSP approaches can come at the cost of expedience,96 which can be a significant impediment to their adoption along the drug delivery pipeline. Additionally, there is a potentially higher person‐power cost due to the increased need for trained mathematicians and life scientists capable of constructing, applying, and translating the sometimes quite complex resulting QSP models.96 In response, we have previously argued that the level of model detail should necessarily line up with the goals and expected outcomes of any given study.97 However, as discussed in the following examples, the recognition that QSP models can be broadly applied to the study of questions and pathologies outside of their original scopes by accounting for the physiological determinants of drug disposition and effects should be seen as a strong response to their higher upfront costs and previously mentioned drawbacks.

For example, the model of Ref. 71 is adapted in Ref. 98 to study the treatment of pneumococcal lung infections and their treatment with antibiotics, thereby applying their model of granulopoiesis to the effect neutrophils have on the immune response. We reexamined the model in Ref. 88 in Ref. 30 and showed that previous G‐CSF PK models (exogenous and endogenous) overestimated the contribution of the renal G‐CSF elimination pathway. Accordingly, we introduced a novel PK model for G‐CSF that explicitly accounts for the bound and unbound G‐CSF concentrations. Using exogenous G‐CSF, chemotherapy‐only, and CHOP14 data, we determined that around 70% of the clearance of G‐CSF is performed by binding to neutrophil receptors.

Although my focus here was to highlight models used to characterize and predict myelosuppression during chemotherapy, within the broader context of the future direction of pharmacological modeling, the studies highlighted in the section on Bridging the Gap suggest a path towards a wider‐scale adoption of the QSP vision through the implementation of the general principles underlying the approaches in Refs. 71, 92, and 30: the incorporation of a systems‐level description of the physiology and the inclusion of molecular‐level considerations when necessary; the addition of PK/PD models based on the interactions of xenobiotics with the physiological system; parameter estimation informed by previously published experiments, with fitting undertaken when no identifiable mechanistic data are available; and the inclusion of variability through the a priori examination of where random elements are likely to effect the system. Not limited solely to predicting the myelosuppressive effects of chemotherapy drugs, the models developed by groups using QSP and QSP‐like approaches are all able to tackle larger questions such as the optimization of regimens, the classification of various patient subpopulations (low‐ and high‐risk groups, etc.), and, most crucially, to explain why their findings hold, by identifying the precise physiological mechanisms responsible for an observed result.

Ultimately, improvements to patient care necessitate innovative science and strong partnerships between academia, industry, and clinicians working towards specific goals (addressing a particular toxicity like myelosuppression) while keeping the bigger picture (developing broadly applicable approaches to improve translational medicine techniques, for example) in focus. Motivated by the exorbitantly high financial, time, and patient‐care costs of attrition along the drug development pipeline, quantitative systems pharmacology has emerged from an increasing number of collaborations between pharmaceutical scientists and mathematicians. In our era of big‐data and rapid experimental advances, the paradigm of “one‐gene, one‐receptor, one‐mechanism”90 is limited in its ability to extensively address current and future medical and drug development problems. Recognition of the role QSP modelers have to play, not only early in the drug development pipeline but perhaps more principally as integrative translational researchers, will positively influence the quality of pharmaceutical research and improve the predictive strengths of the resulting modeling work, thereby productively affecting overall patient care.