 Research Article
 Open Access
 Published:
An Integrative multilineage model of variation in leukopoiesis and acute myelogenous leukemia
BMC Systems Biology volume 11, Article number: 78 (2017)
Abstract
Background
Acute myelogenous leukemia (AML) progresses uniquely in each patient. However, patients are typically treated with the same types of chemotherapy, despite biological differences that lead to differential responses to treatment.
Results
Here we present a multilineage multicompartment model of the hematopoietic system that captures patienttopatient variation in both the concentration and rates of change of hematopoietic cell populations. By constraining the model against clinical hematopoietic cell recovery data derived from patients who have received induction chemotherapy, we identified trends for parameters that must be met by the model; for example, the mitosis rates and the probability of selfrenewal of progenitor cells are inversely related. Within the dataconsistent models, we found 22,796 parameter sets that meet chemotherapy response criteria. Simulations of these parameter sets display diverse dynamics in the cell populations. To identify large trends in these model outputs, we clustered the simulated cell population dynamics using kmeans clustering and identified thirteen ‘representative patient’ dynamics. In each of these patient clusters, we simulated AML and found that clusters with the greatest mitotic capacity experience clinical cancer outcomes more likely to lead to shorter survival times. Conversely, other parameters, including lower death rates or mobilization rates, did not correlate with survival times.
Conclusions
Using the multilineage model of hematopoiesis, we have identified several key features that determine leukocyte homeostasis, including selfrenewal probabilities and mitosis rates, but not mobilization rates. Other influential parameters that regulate AML model behavior are responses to cytokines/growth factors produced in peripheral blood that target the probability of selfrenewal of neutrophil progenitors. Finally, our model predicts that the mitosis rate of cancer is the most predictive parameter for survival time, followed closely by parameters that affect the selfrenewal of cancer stem cells; most current therapies target mitosis rate, but based on our results, we propose that additional therapeutic targeting of selfrenewal of cancer stem cells will lead to even higher survival rates.
Background
Acute myelogenous leukemia (AML) is a cancer of the white blood cells stemming from the myeloid lineage that produces cells including neutrophils and monocytes [1]. On average, four in 100,000 individuals will develop AML, with a median age of 67 years at diagnosis. AML has the highest mortality rate of all of the different types of leukemia [2]. The FrenchAmericanBritish (FAB) cooperative group identified eight different subcategories of AML, M0M7, based on the specific cell population from which AML arises [3]. More recent classification systems now exist that include other criteria besides white cell morphology [4, 5]. Regardless, the majority of patients with AML receive essentially identical induction chemotherapy with drugs that target DNA replication [6]. Unfortunately, for the patients over 55 years of age diagnosed with AML, the 5year survival rate is 10% and the relapse rate is 80% [7]. Given the variability of AML progression among patients, understanding determinants of disease progression will lead to therapeutic advances that result in lower rates of relapse and higher rates of survival.
Computational modeling is an effective tool to personalize therapies for improved patient outcomes. Modelbased personalized therapy has improved several medical interventions including: glucose control in type I diabetes [8, 9]; automatic control of anesthesia dosage [10–12]; and pacemaker control of heart rate variability [13]. Presently, AML treatment is not benefiting from modelbased personalized approaches, in part due to a lack of multilineage AML computational models [14, 15]. Key features of a computational model of hematopoiesis includes: normal formation of mature cells from hematopoietic stem cells (HSCs) [16–24]; abnormal development leading to leukemia [18, 25–27]; treatment using chemotherapy or transplantation [17, 25, 26]; movement of cells between the bone marrow/tissue and the peripheral blood where AML is frequently measured [20, 28–35]; and the feedback signals from cytokines that participate in the regulation of hematopoiesis [20, 25, 26, 28]. Previous models that contain various elements of these key features exist (summarized in Additional file 1: Table S1 in Supplement 1). Here, we focus on the development of a multilineage AML model that captures patient variability and use the model to identify patient subtypes and identify potential novel therapeutic targets.
In this work, we develop a semimechanistic multilineage multicompartment computational model of hematopoiesis and AML. The model describes the interactions of various lineages of hematopoietic stem cells, including neutrophils, lymphocytes, and monocytes. Patient data is used to constrain the model, creating a large number of dataconsistent solutions that differ in their responses to chemotherapy. We apply the constraints from patient data to characterize normal hematopoiesis versus leukemic differentiation and growth in unique parameter sets representing 22,796 simulated patient dynamics. To identify trends in dynamical variability, we clustered the simulated dynamics into thirteen ‘representative patient’ groups. We find that our mathematical model can be used to demonstrate the typical variation in the dynamics of AML progression and patient clinical outcomes. Using sensitivity analysis and correlation analysis, we identified important parameters whose variation most affects the development of leukemia.
Methods
Herein, we develop a semimechanistic multilineage model of leukopoiesis and the formation of AML that includes the following: (1) variability in healthy and malignant white blood cell maturation; (2) interactions between progenitor cells and differentiated cells; (3) interactions amongst macrophages, neutrophils, and lymphocytes; and (4) dynamics of hematopoiesis and AML progression (Fig. 1; model equations in Additional file 1: Supplement 2). This ordinary differential equations (ODE) model of hematopoiesis describes the physiological processes of cell maturation, cell mobility between different blood compartments in the body, interactive responses to cytokines and the cells that produce them, cellular activation and suppression, and cell loss due to natural death or chemotherapy.
Model Configuration
The leukopoiesis process consists of a set of discrete stages of cellular differentiation. There are three different stages of cell development representing the maturation of the stem cells into mature leukocytes (Fig. 1, solid black circles). The first stage contains the hematopoietic stem cells (HSCs) that either asymmetrically selfrenew or differentiate. We assume this HSC compartment consists of many stages of cells that have the capability to differentiate into various lineages, as previous work has shown that the dynamics of mature cells are appropriately represented in a twocompartment system [36–38]. The stem cells proportionally differentiate into distinct progenitors, including common lymphoid progenitors and granulocyte/myeloid progenitors that then clonally amplify their branching lineages [39]. All three lineages are necessary to demonstrate lineage dominance and cellcell interactions amongst the three mostcommonly measured white blood cell types (see Additional file 1: Supplement 3, Table S5 for further justification of these three lineages). Other multilineage models do not consider these three lineages [21–24]. Granulocyte/myeloid progenitors produce mature neutrophils and monocytes, which have relatively short lifespans and are replaced daily. Lymphoid progenitors branch into immature Blymphocytes and Tlymphocytes, yielding more complex and less predictable dynamics than those of myeloid cells. Thus, we focused on studying myeloid cell production, destruction, and deviations associated with myeloid leukemias instead of the less predictable lymphoid leukemias. The process of proliferation and differentiation in our model ensures a physiological mechanism that repopulates the bone marrow after chemotherapy and ensures that the formation of a single mutated cell in the HSC compartment can result in AML.
Mature cells are maintained in three different compartments: the bone marrow or tissue (yellow compartment in Fig. 1), the circulating peripheral blood (red compartment in Fig. 1), and the noncirculating peripheral blood, or marginal pool (blue compartment in Fig. 1). Cells migrate from both the bone marrow and marginal pool to the peripheral blood and are recruited to the bone marrow or marginal pools from the peripheral blood [40–42]. Since the concentration of normal white blood cells in the bone marrow and peripheral blood may vary over several orders of magnitude [40], it is essential to compute the dynamics of the proliferation of cells within the bone marrow and the mobilization of cells. Other multilineage models do not demonstrate movement of cells amongst compartments [21–24], and mobilization is typically limited to neutrophils or monocytes [20, 27–31, 33–35]. Additionally, neutrophils and monocytes also remain in noncirculating peripheral blood pools, or marginal pools, to allow for massive demargination upon a trigger event, such as chemotherapy [40, 43].
Our final model contains 17 states that represent cell populations and 51 parameters. The equations for this model can be found in Additional file 1: Supplement 2, and the model is depicted in Fig. 1. Details of these biophysical processes that regulate cellular dynamics are in Additional file 1: Supplement 3. The parameter values can be found in Additional file 1: Supplement 4. The general form of the model for each cellular state is a sum of all of the production terms, the movement terms, and the loss terms (Eq. 1).
Here, X _{0} is the state of interest; p _{0} is the proliferation term, typically described as asymmetric selfrenewal (circular arrows in Fig. 1; based on [19]; described further in Additional file 1: Supplement 3); s _{−1} is the differentiation (specialization) term from a stem cell or progenitor state, X _{−1}, to the state of interest, X _{0}; s _{0} is the differentiation term from the state of interest, X _{0}, to a more mature state (straight solid arrows in Fig. 1); m _{−1} is the mobilization or margination term from the corresponding bone marrow or peripheral blood state, X _{ b m/p b }, to the state of interest, X _{0} (forward dashed arrows in Fig. 1); r _{0} is the recruitment term of X _{0} to another compartment (backward dashed arrows in Fig. 1); and d _{0} is the loss, or death rate (x’s in Fig. 1). Further details of this model are described in Additional file 1: Supplement 3. A few key physiological details, including feedback and cellcell interactions in the multilineage model, are described below.
Feedback
Cells secrete complex cytokines and chemokines that play a pivotal role in hematopoiesis. Cytokines regulate different processes that occur in the hematopoietic system, such as selfrenewal and movement [44]. Mutations in the production or signaling of cytokines can lead to negative outcomes, such as the formation of AML [45–48]. These cytokines are not specified in our model, but are assumed to be proportionate to the number of cytokine releasing cells (i.e. macrophages). This is incorporated into our model using MichaelisMenten type kinetics to modify the rates of cytokinetargeted processes. The general form of negative feedback is in Eq. 2 and the general form of positive feedback is in Eq. 3.
Here, k is the MichaelisMenten constant; C is the concentration of the cell state that produces cytokines for feedback; r is the rate the feedback modifies; and X is the cell state targeted by cytokines. To model the formation of AML from monocytes, we simulated the multiplehit hypothesis and assumed all feedback signals did not effect a single normal monocyte progenitor, as expected from literature [49].
Inhibitory feedback prevents cells from growing uncontrollably in any state. Thus, our model has several inhibition mechanisms to regulate the concentration of cells. In particular, inhibition feedback modifies every rate except cell death in our model. All movement is inhibited by the concentration of cells in the compartment to which cells are moving to prevent overcrowding using MichaelisMenten type kinetics (Eq. 2). Proliferation and differentiation are hindered by inhibiting the associated selfrenewal probability of stem cells [19]. For progenitor cells, this inhibition occurs from cytokines produced by mature cells of the same lineage in the peripheral blood [50]. However, for stem cell proliferation and differentiation, this inhibition occurs from a scaled combination of the concentration of stem cells to ensure that sufficient stem cells are in the hematopoietic system and the concentration of cells in the bone marrow (yellow compartment in Fig. 1) do not exceed the capacity of the marrow. All negative feedback signals are pictorially described in Additional file 1: Figure S1. The derivation of both the form for asymmetric selfrenewal and the capacity of the bone marrow are in Additional file 1: Supplement 3.
In contrast, several physiological processes in hematopoiesis are triggered by positive feedback. In response to an inflammatory event, mature white blood cells proliferate to accommodate this response. We include two states in our model that demonstrate the effect of debris clearance due to chemotherapy. Monocytes are recruited into the bone marrow to become activated macrophages to help clear excessive cellular debris [41], which we have modeled as the ‘Apoptosis’ state. Macrophages assist in recruiting other white blood cells, such as lymphocytes, neutrophils, and monocytes, during inflammation or other diseases [41], which we have also incorporated into the model. Positive feedback is used in our model to demonstrate several processes, including inducing proliferation of macrophages by apoptotic debris; promoting the recruitment of neutrophils, lymphocytes, and other monocytes into the tissue in response to high levels of macrophages; and increasing the clearing rate of apoptosis due to a high level of macrophages. All positive feedback signals are pictorially described in Additional file 1: Figure S2. Additional file 1: Table S4 in Supplement 2 describes each of the feedback processes in detail.
Results
Parameter selection
The model has 51 parameters and 17 states that describe the cell concentrations of various white blood cells and leukemia (Fig. 1) and reflects seven qualitative and quantitative features found in patient responses to chemotherapy from literature (Table 1). Seven of these parameters were found from literature: the death rates of progenitor cells [51], mature neutrophils [28, 52, 53], mature lymphocytes in the bone marrow [54] and peripheral blood [25, 55], mature monocytes in the bone marrow [28, 56, 57] and peripheral blood [56], and macrophages [58]. We fit the remaining 44 uncertain parameters to reflect the seven dynamic acceptability criteria (Table 1) of our model; these parameter values are in Additional file 1: Supplement 4. Parameter optimization proceeded in sequence by first identifying dataconsistent parameters in the unilineage models and then extending the identified parameters to the integrated multilineage model (Fig. 2).
Qualitative features and properties and quantitative patient data were used to test models and optimize parameters. However, given the amount of patienttopatient variability, this only informs regions of the parameter space consistent with observed clinical data as opposed to optimized parameter values. We explored for parameters in a nested approach, where unilineage models were evaluated independently before being combined with the multilineage model (Fig. 2). Specifically, we were interested in demonstrating a realistic range of patient white blood cell counts and the rates that dictate those counts. Manesso et al. also used a “divide et impera" stepwise approach to parameterizing their multilineage model to demonstrate variability in [23]; however, a key difference of their model to ours is that we searched for both initial conditions and parameter sets that met expected steady state configurations instead of using various parameter set perturbations on one initial condition. In order to find these initial conditions and parameter values, we searched over large parameter ranges for the 44 undetermined parameters using MATLAB’s lsqnonlin search algorithm. Selfrenewal fractions were bound between 0.5 and 1 to maintain stem cell mass without becoming zero. Nominal values for the remainder of the parameters were approximated based on expected dynamics. We searched between 1/5 to 5 times these nominal values to determine the bounds of the parameter space. We simulated 100,000 Latin Hypercube samples, an efficient, semirandom but uniform method of sampling, in log space for each unilineage model (Fig. 3, n _{ LHS }). Almost all of the Latin Hypercube samples produced nonzero equilibrium solutions for further analysis (Fig. 3, n _{ NTEQ }), but not all monocyte simulations produced equilibrium solutions due to the stiffness of the system. In order to determine the parameter ranges that met the dynamical criteria from Table 1 (further detailed in the following section), we simulated chemotherapy for 7 days on each of the nontrivial equilibrium solutions (details of chemotherapy in Additional file 1: Supplement 3). Acceptable solutions met all seven of the dynamical acceptability criteria (n _{ F }=177 for neutrophils, n _{ F }=1796 for lymphocytes, and n _{ F }=1049 for monocytes). Figure 2 depicts the process, and Fig. 3 shows the intersection of all acceptable parameter regions.
Physiological capabilities of model
Dataconsistent simulations not only fit the patient data, but they also meet a number of qualitative and quantitative criteria (summarized in Table 1). These seven capabilities were mostly determined from literature, though we also used data from Dr. Robert Nelson of the Indiana University Simon Cancer Center from patients who received hematopoietic stem cell transplants (HSCT) to determine additional empirical evidence of the patient dynamics (IRB Protocol# 1011009928). Further information on the patient data is included in Additional file 1: Supplement 5. The patient data was used to determine cell concentration bounds, but were not used to validate the time course of patient recovery after transplantation. Thus, a donor infusion was not modeled. The pharmacokinetics and pharmacodynamics of chemotherapy was simulated as described in Additional file 1: Supplement 3, and was implemented as a ‘chemo’ term in the model equations (Additional file 1: Supplement 2). The details of these criteria are below.

1.
Peripheral blood recovers > 80 cells/ μ L: This ensures that the peripheral blood sufficiently recovers after chemotherapy treatment, and the final peripheral blood cell count must be > 80 cells/ μL [58].

2.
Stem cells recovery > 1 cell/ μ L: The final stem cell concentration must be between 1 and 100 cells/ μL. If the stem cell concentration is within the acceptable bounds, the stem cells should also show recovery after chemotherapy occurring [58].

3.
Final dynamic values are within acceptable ranges: The final unilineage peripheral blood cell counts after chemotherapy must be within the upper limits of expected normal peripheral blood cell counts [40, 59]. Neutropenia (low number of neutrophils in the peripheral blood) is very common after undergoing chemotherapy [60], but neutrophilia mostly occurs if another underlying disease exists to cause the high number of neutrophils. We do not want to model neutrophilia during recovery from chemotherapy because that is not a typical response in patients with AML. The values each cell state must be within are depicted in Additional file 1: Supplement 4, Table S6.

4.
Marginal pool is within one order of magnitude of peripheral blood: The final concentration of the cells in the marginal pool must be within one order of magnitude of the final concentration of the peripheral blood because the size of the marginal pool and peripheral blood compartments are approximately the same [43, 61, 62]. Though a marginal pool can exist for lymphocytes, it is much smaller than that of neutrophils and monocytes [43]. Thus, we do not incorporate a lymphocyte marginal pool in our model.

5.
Cell counts reduce with chemotherapy to <20 % of original value: The concentration of cells remaining at the end of chemotherapy should be less than 20% of the concentration of cells before chemotherapy begins [63].

6.
Recovery overshoot < 12 times value five days after overshoot: Patients experience lymphopenia, or low white blood cell concentrations, during chemotherapy. After the chemotherapy regimen is completed, patients’ cell counts almost always increase. In some instances, patient cell counts increase and peak above their final homeostatic cell concentration. Here, we define this overshoot level as the maximum concentration of cells after chemotherapy administration is completed. If the overshoot level is more than twelve times the concentration of cells five days after the overshoot (or the concentration of cells at the end of the simulation, whichever comes first), the solution is rejected. This overshoot was the maximum observed in the HSCT patient data from the Simon Cancer Center (Additional file 1: Figure S6 in Supplement 5).

7.
Amplitude of the second peak of oscillations deviates < 18% from the amplitude of the first peak: Several simulations produced oscillations. Though some small oscillations are reasonable and can occur physiologically, large oscillations are unlikely because they will not maintain normal homeostasis in the body and most patients who are treated for AML do not report large oscillations in cell counts. We manually classified 100 simulations for sufficiently dampened oscillations. We calculated that a simulation was sufficiently dampened if the cell concentration of the second peak value postchemotherapy was less than an 18% decrease from the first peak postchemotherapy. Thus, we implemented this cutoff for all of our simulations. An example of a simulation that does not meet this criteria and a simulation that does meet this criteria are shown in Additional file 1: Figure S7 in Supplement 5.
For each unilineage parameter search, two criteria were salient in determining which models met the physiological features required based on responses to chemotherapy. The remainder of the criteria are mostly robust in constraining the parameter space. One of the salient criteria was the lack of oscillations for each unilineage model. This implies that the need to maintain a steadystate concentration of cells for each of the lineages independently is one of the key constraints in normal hematopoiesis and is consistent with parameters that lead to balanced feedback. The second criteria that largely determined the acceptable unilineage models were the final peripheral blood acceptable dynamics of the neutrophil unilineage model, stem cell recovery for the lymphocyte unilineage model, and sufficient response to chemotherapy for the monocyte unilineage model (indicated with asterisks in Fig. 3). This reflects the necessity of maintaining high concentrations of neutrophils in homeostasis and the ability of cells to repopulate cell populations for appropriate recovery. If we relaxed the constraints imposed on the model set by the features of normal physiological dynamics, the acceptable parameter sets would be larger as well.
Parameter constraints with separatrix method
We identified several relationships among pairs of uncertain parameters in the models that met the dynamic capabilities of unilineage models. This demonstrates parameter relationships that maintain the normal physiological responses to stimuli, and reveals certain physiological processes that are directly related in a specific manner to each other. We discovered that negative exponential relationships exist between pairs of parameters, especially between the fraction of cells that selfrenew and the mitosis rates of cells (Fig. 4). This implies that either the mitosis rate or the probability of cell renewal needs to be high for a progenitor cell to sustain its population. If both of these rates are high, however, instability occurs. Traditional methods of separating data in multiple dimensions into groups do not define the bounds of the groups, such as in clustering [64], or do not separate multiple linear separations observed in Fig. 4, such as in support vector machines [65]. A computationally efficient method of separating multiple linear constraints is needed. Thus, to constrain the parameter space based on these relationships, we developed a separatrix method that defined the bounds of the relationships between parameters.
The separatrix method distinguishes between parameter regions of acceptable solutions and regions of unacceptable solutions, and in logspace, this manifests as a clear separatrix. We first constrain to the acceptable range of each parameter (dashed black box in Fig. 4). For each unique pair of unilineage parameters, we then divide the parameter range rectangle into 10×10 equally sized bins. If the acceptable parameter set is uniformly distributed in parameter space, every bin has a cutoff of (number of acceptable simulations)/100 acceptable parameter sets. Thus, each of the 100 bins is either acceptable or unacceptable depending on whether more acceptable parameter sets existed in the bin than the cutoff value. We then further constrain this acceptable parameter density by removing corners of the range rectangle that did not maintain a cutoff density of acceptable parameter sets (dashed green lines in Fig. 4). Thus, we found a set of inequalities for each unique pair of parameter sets in the unilineage parameterization in log space. We used this separatrix constraint on the parameter search space for multilineage parameterization.
The separatrices of each pairwise parameter relationship confirm parameter relationships from previous works, and we identified additional required relationships in the acceptable solutions. As Getto et al. [38] and Stiehl et al. [52] show, the mitosis rate of stem cells and the probability of selfrenewal are inversely related. However, we found that this also extends to all progenitor cells (Fig. 4 ac). Specifically, a linear relationship emerges between the log of the selfrenewal probability and the log of the mitosis rate. Solutions that have larger mitosis rate or selfrenewal probabilities produce oscillations (red x’s in Fig. 4), which we constrained against in our dynamical acceptability criteria (Criteria 7 in Table 1). We also found a similar inverse linear relationship between the log of the selfrenewal probability of neutrophil progenitors and the log of the homeostatic constant of neutrophils in the peripheral blood; solutions that have lower values of either of these two parameters becomes unacceptable (Fig. 4 d). Overall, this finding means that a specific relationship is maintained between the probability of selfrenewal and the mitosis rate for all cells that are capable of selfrenewing to ensure that the cell population remains steady. Additionally, due to the high concentration of neutrophils, a constrained feedback mechanism exists between the selfrenewal probability of neutrophil progenitors and the cytokines produced by neutrophils in the peripheral blood to maintain appropriate homeostatic concentrations of those cells. Overall, we confirmed the dependence on the mitosis rate and selfrenewal probability of stem cells found by other groups (Fig. 4 e), but we also found additional dependencies on these same parameters of progenitor cells, as well as a specific dependence of the cytokines produced by neutrophils in the peripheral blood and the selfrenewal probability of neutrophil progenitor cells.
Multilineage acceptable solutions
Following analysis of unilineage models, we integrated the unilineage models into a single multilineage model of leukopoiesis and AML to identify the constraints on the rates that regulate normal hematopoietic and leukemic proliferation. To evaluate model performance, we required the complete model to be consistent with all of the constraint data that informed the unilineage models. We used the specific parameter relationships gleaned from using separatrices on the unilineage models to select a set of models that meet the capabilities of the unilineage models. Thus, we simulated 100,000 Latin Hypercube samples of the 44 undetermined parameters in the constrained multilineage parameter space (Fig. 5, n _{ LHS }). 40,368 parameter sets maintained homeostasis (Fig. 5, n _{ NTEQ }); of which, 22,796 simulations with specific parameter sets and initial homeostatic conditions were deemed acceptable simulations using the same seven dynamic criteria from Table 1 (intersection of acceptable parameter spaces in Fig. 5, n _{ F }). The lymphocyte unilineage models caused the most constraints on the full multilineage models, specifically to restrict the oscillations and promote recovery postchemotherapy in lymphocytes (Fig. 5). This is potentially due to the large concentration of lymphocytes and their rapid selection in the thymus [54] in contrast to their long halflives in the peripheral blood [25, 55]. Overall, in the final multilineage model, the acceptable parameter sets represent 22,796 virtual healthy patients. All of the mean and ranges of the parameters that were varied are in Additional file 1: Supplement 4, Table S8. Correlation coefficients of the initial conditions and parameter values of the multilineage model are shown in Additional file 1: Figure S5 in Supplement 4.
A comparison of the simulations of the multilineage solutions to the simulations of the unilineage solutions is shown in Fig. 6 ad. In the stem cell and all peripheral blood cell compartments, the unilineage models allow for a greater spread of cell dynamics than in the multilineage models. Specifically, the overshoot of recovery of cells postchemotherapy is dampened in the multilineage solutions. This means that the combined feedback from cytokines from different lineages has a stabilizing effect. Additionally, in observing the distributions of the concentration of cells of the multilineage solutions in comparison to the unilineage solutions at specific time points (Fig. 6 eh), the multilineage solutions (blue lines in Fig. 6 eh) have lower cell concentrations at the early time points (dotted and dashed lines in Fig. 6 eh) but higher cell concentrations at later time points near homeostasis (solid lines in Fig. 6 eh). This implies that in order to maintain acceptable cell dynamics in a multilineage system, the concentration of cells needs to be at a higher homeostatic concentration than if only considering the unilineage model. The overshoot of the lineages in Fig. 6 ad are all well within patient variation (Additional file 1: Figure S6 in Supplement 5); though, our simulation dynamics may be too constrained in comparison to patient data. Forty out of the 47 patients (85%) in Additional file 1: Figure S6 do not demonstrate neutrophil overshoot, but none of our multilineage neutrophil simulations demonstrate this overshoot. We find that constraining the unilineage models to multilineage data may have constrained the multilineage model more than patient overshoot data reflects. We additionally test the robustness of the multiple lineages in the model by perturbing progenitor neutrophils and testing the dynamics before reaching homeostasis on the different lineages when coupled together compared to when decoupled. We found that the multilineage model tempers the effects of this perturbation in the neutrophil lineage but causes lymphocyte dynamics to change (Additional file 1: Supplement 6, Figure S8).
Sensitivity analysis reveals importance of selfrenewal probability and cytokines
We carried out a sensitivity analysis on the multilineage solutions to determine the most important parameters whose variation most greatly affected output. Using a partial rank correlation coefficient (PRCC) method [66–68], we determined the parameters that were most impactful in determining concentrations of stem cells, neutrophils in the peripheral blood, lymphocytes in the peripheral blood, and monocytes in the peripheral blood. We found that, in general, the probability of selfrenewal and the MichaelisMenten constant that modifies both the rate of selfrenewal probability and the movement of mature cells into the peripheral blood were the most important parameters for their respective cell states. This indicates that the probability of a cell to selfrenew is very sensitive to changes in the system and can cause very drastic outcomes in the final cell concentration if this rate is affected, which corroborates with findings in MarciniakCzochra et al. [19]. Additionally, for stem cells, the feedback cytokines from the concentration of cells in the bone marrow and the concentration of stem cells were very important parameters. These constants also modify the probability of selfrenewal in the stem cells. Finally, many parameters that are associated with neutrophil homeostasis (selfrenewal probability of neutrophil progenitors, mitosis rate, MichaelisMenten constant, and mobilization rate) all were very important in the stem cell concentration. This is probably due to the large concentration of neutrophils affecting the homeostasis values of stem cells. The PRCC results are in Additional file 1: Supplement 6.
Representative patient clusters
In order to simulate cancer in a feasible number of patients, we clustered the 22,796 solutions in the normalized dynamic space. This means that we normalized each state in each parameter set simulation by the maximum value simulated. We clustered these normalized simulations across all states and 162 time points (daily from day 7 to day 150 and once every subsequent 30 days until day 300) using kmeans, where chemotherapy was applied from day 7 to day 1. We found that 13 clusters maximized the dynamical similarities within a cluster over all states and minimized the overlap in the dynamics across all states between clusters. These 13 clusters represent 13 subgroups of patients, and the centroid of each of these clusters is the most representative patient for each cluster. For more information on determining the number of clusters, see Additional file 1: Supplement 7. In Fig. 7 a and b, we show how the dynamics are clustered across some of the most varied dynamical regions from Fig. 6. Stem cells group into two distinct behaviors; stem cell concentrations reach homeostasis within two months postchemotherapy for all but three clusters (Fig. 7 c). Normalized neutrophil dynamics are fairly constrained, and the majority of its variation occurs during the rapid recovery of neutrophils postchemotherapy (Fig. 7 d). However, the overall dynamics of the normalized clustered representatives demonstrate the variability in both lymphocytic (Fig. 7 e) and monocytic dynamics (Fig. 7 f). In general, the various clusters differ from each other primarily by the normalized monocyte concentration in the peripheral blood postchemotherapy, overshoot in the concentration of lymphocytes in the peripheral blood, recovery time of neutrophils in the peripheral blood, and the homeostatic concentration of stem cells relative to their initial concentration postchemotherapy. A correlation analysis did not reflect any significant differences in the independent parameter values amongst the thirteen clusters. However, simulations within individual clusters may have parameter dependencies that are different between clusters. Real patients can potentially be assigned to clusters based on their dynamics following one round of chemotherapy, though outlier patients need to be sufficiently explored to encompass true patient dynamics.
Modeling leukemia
AML is a derivative of common granulocyte macrophage progenitor cells [3, 69]. Thus, leukemia is modeled parallel to the manner in which monocytes are modeled. We assume that the cancer stem cells selfrenew and proliferate at the same rate as monocytic progenitor cells. These then differentiate into mature cancer within the bone marrow, which can mobilize into the peripheral blood. The only difference in the model between leukemic stem cells and monocytes is that the leukemic cells do not maintain homeostasis, so all feedback that regulates selfrenewal [18, 70], differentiation, and movement were removed and the cancer cells cannot be recruited back to the tissue. Since we have modeled feedback from various signals, the mutated cancer cell that grows indefinitely is the product of several mutations that cause all feedback to be inhibited, as currently hypothesized in literature [49].
We simulated AML with one cancer stem cell, and allowed the leukemia to clonally accumulate in each of the 22,796 simulated patients (Fig. 8 a), as expected from literature [71]. The variation shows that some patients may survive only a few months after the cancer starts (i.e., pink cluster lines in Fig. 8 a; gray lines indicate death) and some patients may live close to one year without any other intervention (i.e., green cluster lines in Fig. 8 a). In order to determine how long patients in each cluster would live without any treatment, we assumed patients would survive until the concentration of cancer cells in either the peripheral blood or bone marrow was greater than 3×10^{5} cells/ μL (Fig. 8 ad). This is the theoretical maximal concentration of cells in the bone marrow calculated in Additional file 1: Supplement 3, given the diameter of monocytes and assuming the volume in the bone marrow and peripheral blood are the same. Figure 8 b shows the AML growth as a percentage of the peripheral blood. AML is diagnosed when over 20% of the bone marrow or peripheral blood contains cancer blasts, or immature cancer cells [72], and our model indicates that patients can remain undiagnosed with AML for weeks to months, as previously suggested [73]. The growth of AML is simulated to stop when the theoretical maximum concentration of bone marrow cells is reached. In Fig. 8 b, patients with slower growing cancers die with AML comprising a higher percentage of their peripheral blood. Figs. 8 cd indicate that patients within some clusters have a much better prognosis than others, and generally correlates with a slower recovery of monocytes postchemotherapy (corresponding colors in Fig. 7 e).
For a simple validation to ensure that cancer is not growing too slowly in our model, we can compare overall survival times of the simulations without treatments to trials in which patients received lowdose treatments. The majority of clinical patients who are unable to receive traditional chemotherapy and received hydroxyurea as a treatment postdiagnosis with AML died within one year of diagnosis even with favorable stratification [74]. We find the survival time in all of our simulations to be less than four months in after diagnosis of AML >20% blasts (Fig. 8 b). Thus, cancer growth fits within an expected bound of patient survival.
To test which parameters correlate with the survival time of patients with AML, we correlated survival time to the four parameters that mediate cancer in our model: (1) the death rate of all cancer cells (dMc), (2) the fraction of cells that selfrenew (aMc2), (3) the mitosis rate of cancer progenitor cells (mrMc2), and (4) the mobilization of cancer cells into the peripheral blood (mbMc). For the 22,796 cancer simulations, these parameter values were assumed to be the same as those of monocytes in healthy hematopoiesis. In order to determine which of the parameters that determine survival time the most, we sampled 625 Latin Hypercube Samples for each of the 13 representative patients. For this sampling, we varied the parameters, dMc, aMc2, mrMc2, and mbMc within the ranges shown in Table 2. These bounds were chosen to ensure the death and mobilization rates were slower than those of normal monocytes, but mitosis rates are higher. We computed a correlation of all of the diagnosed cancers in these simulations (bone marrow or peripheral blood blast concentration >20% and cancer stem cell concentration >1 cell/ μL) to the number of days of survival. We found that across all clusters, the mitosis rate of cancer was correlated with survival without chemotherapy. The overall negative correlation of 0.7431 (Fig. 8 e) indicates that the higher the mitosis rate, the more likely the patient will die earlier without treatment. We constrained this same search by fixing the mitosis rate. and found that the next most correlated parameter to survival time was the probability of selfrenewal of cancer, with an overall correlation of 0.8774 (Fig. 8 f). This result aligns with the inverse relationship found between mitosis rates and selfrenewal probability (Fig. 4). Stiehl et al. have similarly found that high proliferation rates and high selfrenewal rates can also lead to earlier death [52]. We simulated AML derived from neutrophil progenitors analogously to simulating monocytic AML and confirmed these results (Additional file 1: Supplement 8).
Discussion
When a patient is diagnosed with AML, physicians assess the white blood cell counts of the patient to guide the course of action for treating the leukemia. Many patients are successfully treated with a stem cell or bone marrow transplant from a matched donor. However, for those patients who are unable to receive a transplant, chemotherapy is administered periodically to lower the cancer load on the patient. Generally, a high white blood cell count makes the patient a candidate for chemotherapy. However, it is likely that factors other than the absolute cell count are important for determining a patient’s prognosis. Using a multilineage model of the formation of white blood cells, we find that (1) certain physiological rate relationships are necessary to prevent unstable cell population growth and (2) the rate of growth of the cancer is an important prognostic factor in determining the survival time of patients.
We determined that there are specific physiological rates and cell concentrations that are crucial in maintaining homeostasis. As expected, progenitor cells are very important in homeostatic mechanisms. A sensitivity analysis demonstrated that the probability of selfrenewal of all progenitor cells and parameters that modify this probability are the most important parameters in hematopoietic dynamics. Other groups have found similar results [19] in the context cancer [75, 76] and specifically AML [77]. Additionally, we found an inverse relationship between the selfrenewal probability of stem cells and the mitosis rate of stem cells (Fig. 4 e), corroborating findings of other groups [38, 52]. This relationship between mitosis rate and selfrenewal probability also exists for progenitor cells (Fig. 4 ac), which was not found previously. In general, solutions that have high mitosis rates or selfrenewal probabilities lead to oscillations in the cell concentrations. This oscillation is highly undesirable and can occur in diseases such as cyclic neutropenia. Previous work has shown that “reentry" into the stem cell compartment was found to be one of the factors that control oscillation in cyclic neutropenia [22]. Our work suggests that the analogous selfrenewal probability and mitosis rate of the same cell, whether it is the stem cell or another progenitor cell, are both crucial to control to prevent diseases such as cyclic neutropenia. Chemotherapy that targets only one of these rates may not be sufficient in controlling physiological oscillations. Thus, we recommend that the feedback mechanisms that govern the selfrenewal probabilities are explored as potential pharmaceutical targets in AML treatment.
We found that peripheral blood concentration levels of neutrophils and lymphocytes are important factors in maintaining homeostasis in our model (Fig. 6), specifically in the multilineage context. This was evidenced in three ways. First, when we identified the capabilities of our multilineage model, we found that lymphocyte dynamics constrained the range of acceptable solutions for all three lineages (Fig. 5). This is likely due to the high selection rate of lymphocytes in the thymus (dL3, which is modeled as the tissue/bone marrow compartment) [54] and low death rate of lymphocytes in the peripheral blood (dL3pb) [25]; this can lead to fast dynamic changes if not constrained appropriately. Second, the inverse relationship between the selfrenewal probability of neutrophil progenitors and the homeostatic term for neutrophil progenitors is a novel relationship in hematopoietic modeling that tightly regulates healthy neutrophil behavior (Fig. 4 d). This means that in order to maintain healthy levels of neutrophils at homeostasis, either the fraction of neutrophil progenitor cells that selfrenew has to be low or the feedback mechanism that maintains neutrophil homeostasis has to have a low threshold for turning selfrenewal off. Thus, selfrenewal is very tightly regulated to ensure that cells do not grow indefinitely. To reduce the cancer load for potential treatment, physicians could target the feedback mechanisms associated with self renewal probability, in particular focusing on cytokine signaling. Third, the large concentration of neutrophils cause the parameters associated with neutrophils to be important parameters in determining the concentration of stem cells, as determined by global sensitivity analysis. We confirmed the importance of neutrophils in relation to lymphocytes in the multilineage model by testing the robustness of the multilineage model in Additional file 1: Supplement 6. We also found that the neutrophils in our multilineage model might be overconstrained by comparing to clinical overshoot data in comparison to the unilineage models. Specifically, the neutrophil overshoot of the multilineage model only reflected the neutrophil overshoot in 85% of clinical data (Additional file 1: Figure S6 in Supplement 5). Thus, our model encompasses most of the dynamics of the patient population. Furthermore, this reflects the importance of the sensing mechanism of each of these cells to maintain appropriate homeostatic levels. When we modeled cancer, removing these homeostatic terms from our model allowed the cancer to grow indefinitely, as expected. Thus, treatment that targets the feedback receptors on cancer cells can drastically help reduce the cancer load.
Individual patients are likely to develop unique phenotypes of AML. The dynamics and relative ratios of the absolute concentrations of neutrophils, lymphocytes, and monocytes vary widely across the patient population. Furthermore, phenotypic variability of leukemia may depend on the specific source cell within the population of common granulocyte progenitors that first becomes cancerous. Here, we develop a mathematical model that is able to describe patient variability in AML. We find that though cancer will form when homeostatic mechanisms are altered, additional mutations that increase the mitosis rate of cancer will reduce the survival of the patient without intervention (Fig. 8 e). Current chemotherapy regimens inhibit DNA replication, which corresponds to inhibiting the mitosis rate of the cancer [78]. However, we also find that if the mechanism that determines the probability of selfrenewal of cancer cells is mutated, then this can be an additional target for pharmaceutical treatments (Fig. 8 f), similar to what has been found in other work [75, 77].
The multilineage model developed in this work can be modified to explore mechanisms governing hematopoiesis and leukopoiesis. For example, the model could be adapted to discriminate amongst chemotherapy regimens on simulated patients to identify characteristics of patients that would benefit from certain regimens. Then, this could lead to identification of treatment schedules optimized for individual patients. The patient clusters could be used to stratify real patients by using the dynamics of patients’ responses to chemotherapy to match to a specific simulation and its outcome. Additionally, since we found that the probability of selfrenewal is a potential secondary target for chemotherapy, various specific feedback mechanisms could be incorporated to identify the exact cytokine signal that would be the best target for therapy. This could be further used to identify the feedback signals that are most sensitive to cause AML growth. Furthermore, the linear relationship in loglog space between mitosis rates and selfrenewal probabilities that was found in our model could be tested experimentally to ensure these relationships exist, though mediating cytokines properly is very difficult. Though we did not find specific relationships amongst the different lineages, the signaling mechanisms that control bone marrow size could also be explored experimentally and compared with the model to regulate overproduction of a particular cell lineage with respect to other lineages. Specifically, GMCSF is often used to regulate the stem cell production of neutrophils and macrophages [79]. This could be manipulated to identify effects in lymphocyte counts and the production of AML. Overall, the multilineage model we present here can be extended to characterize many aspects of hematopoiesis.
The multilineage model of hematopoiesis and leukopoiesis developed in this work can be readily adapted and expanded to incorporate many other immunological effects. Various groups have already modeled bone marrow transplantation and its potential complications [31, 33, 34]. Using our model, bone marrow or stem cell transplantation and transplant rejection can be integrated to predict graft rejection, and different lymphocyte subtypes can be added to the model, such as natural killer cells, to model an alternative outcome of transplantation: graft versus host disease. Lymphocytic leukemias can be explored by altering the homeostatic mechanisms of lymphocytes; analogously, myelodysplastic syndromes and minimal residual disease could be further characterized by determining parameter changes that lead to appropriate model behavior. Additionally, a wide array of other immunological diseases or the wide array of patient responses (including a large Tcell repertoire) can be adapted into this multilineage model to characterize the complexity of normal function and other immunological diseases. More specifically, future work could characterize the mechanisms that cause spontaneous remission of AML in the presence of bacterial infection [80].
Conclusions
The multilineage mathematical model we have created of hematopoiesis and leukopoiesis can help identify how individuals differ in their white blood cell and leukemia production. This model is useful for several reasons. We have determined crucial parameters and parameter relationships that can be used as potential drug targets for both AML and other potential immunological disorders. For many chemotherapeutic drugs, DNA replication is targeted [78], which aligns with targeting the mitosis rate. However, a combination therapy that also addresses the probability of a cancer cell to selfrenew could be potentially helpful for patients whose cancer becomes resistant to the initial therapy. There is no current way to experimentally determine how these different lineages interact and limit each other. Thus, this multilineage model is a very powerful tool that can aid in understanding how blood forms normally and misforms into AML in individual patients. In addition, the model captures multiple dynamics that represent specific patient subgroups. By clustering the wide population of individual differences, we find that one advantage of this multilineage model is that it can be readily extended to investigate personalization of treatment schedules of individual patients to prolong overall survival.
References
 1
Koeffler HP, Golde DW. Acute myelogenous leukemia: a human cell line responsive to colonystimulating activity. Science. 1978; 200:1153–4.
 2
Epidemiology S, Program ERS. Research Data (19732012), National cancer institute, DCCPS, Surveillance research program, Surveillance systems branch. 2015. www.seer.cancer.gov. Accessed 10 Jan 2017.
 3
Bennett JM, Catovsky D, Daniel MT, Flandrin G, Galton DA, Gralnick HR, Sultan C. Proposals for the classification of the acute leukaemias. frenchamericanbritish (fab) cooperative group. Br J Hematol. 1976; 33:451–8.
 4
Grimwade D, Walker H, Harrison G, Oliver F, Chatters S, Harrison CJ, Wheatley K, Burnett AK, Goldstone AH. The predictive value of hierarchical cytogenetic classification in older adults with acute myeloid leukemia (aml): analysis of 1065 patients entered into the united kingdom medical research council aml11 trial. Blood. 2001; 98:1312–9.
 5
Vardiman JW, Thiele J, Arber DA, Brunning RD, Borowitz MJ, Porwit A, Harris NL, Beau MML, HellströmLindberg E, Tefferi A, Bloomfield CD. The 2008 revision of the world health organization (who) classification of myeloid neoplasms and acute leukemia: rationale and important changes. Blood. 2009; 114:937–51.
 6
Kimby E, Nygren P, Glimelius B. A systematic overview of chemotherapy effects in acute myeloid leukaemia. Acta Oncol. 2001; 40:231–52.
 7
Baron F, Storb R. Hematopoietic cell transplantation after reducedintensity conditioning for older adults with acute myeloid leukemia in complete remission. Curr Opinions Hematol. 2007; 14:145–51.
 8
Parker RS, Doyle FJ, Peppas NA. A modelbased algorithm for blood glucose control in type i diabetic patients. IEEE Trans Biomed Eng. 1999; 46:148–57.
 9
Wong XW, SinghLevett I, Hollingsworth LJ, Shaw GM, Hann CE, Lotz T, Lin J, Wong OS, Chase JG. A novel, modelbased insulin and nutrition delivery controller for glycemic regulation in critically ill patients. Diabetes Technol Ther. 2006; 8:174–90.
 10
Dua P, Dua V, Pistikopoulos EN. Model based parametric control in anesthesia. Eur Symp Comput Aided Process Eng. 2005; 20:1015–20.
 11
Mendez JA, Torres S, Reboso JA, Jagannivas HR, Hepsiba D. Modelbased controller for anesthesia automation. IEEE Int Conf Autom Sci Eng. 2009;379–84.
 12
Jagannivas N, Hepsiba D. Control of anaesthesia concentration using model based controller. Int J Innov Res Technol. 2014; 1:237–44.
 13
Bogdan P, Jain S, Marculescu R. Pacemaker control of heart rate variability: A cyber physical system perspective. ACM Trans Embed Comput Syst (TECS); 12(1s):1–22.
 14
Chakrabarty A, Pearce SM, R P Nelson J, Rundell AE. Treating acute myeloid leukemia via hsc transplantation: A preliminary study of multiobjective personalization strategies. Am Control Conf (ACC). 2013:3790–5.
 15
Noble SL, Sherer E, Hannemann RE, Ramkrishna D, Vik T, Rundell AE. Using adaptive model predictive control to customize maintenance therapy chemotherapeutic dosing for childhood acute lymphoblastic leukemia. J Theor Biol. 2010; 264:990–1002.
 16
Peng CA, Koller MR, Palsson BO. Unilineage model of hematopoiesis predicts selfrenewal of stem and progenitor cells based on ex vivo growth data. Biotechnol Bioeng. 1996; 52:24–33.
 17
Scholz M, Engel C, Loeffler M. Modelling human granulopoiesis under polychemotherapy with gcsf support. J Math Biol. 2005; 50:397–439.
 18
Michor F, Hughes TP, Iwasa Y, Branford S, Shah NP, Sawyers CL, Nowak MA. Dynamics of chronic myeloid leukaemia. Nature. 2005; 435:1267–70.
 19
MarciniakCzochra A, Stiehl T, Ho AD, Jager W, Wagner W. Modeling of asymmetric cell division in hematopoietic stem cells–regulation of selfrenewal is essential for efficient repopulation. Stem Cells Dev. 2009; 18:377–85.
 20
Ho T, Clermont G, Parker RS. A model of neutrophil dynamics in response to inflammatory and cancer chemotherapy challenges. Comput Chem Eng. 2013; 51:187–96.
 21
Colijn C, Mackey MC. A mathematical model of hematopoiesis–i. periodic chronic myelogenous leukemia. J Theor Biol. 2005; 237:117–32.
 22
Colijn C, Mackey MC. A mathematical model of hematopoiesis–ii. cyclical neutropenia. J Theor Biol. 2005; 237:133–46.
 23
Manesso E, Teles J, Bryder D, Peterson C. Dynamical modelling of haematopoiesis: an integrated view over the system in homeostasis and under perturbation. J R Soc Interface. 2012;10(80).
 24
Székely T, Burrage K, Mangel M, Bonsall MB. Stochastic dynamics of interacting haematopoietic stem cell niche lineages. PLOS Comput Biol. 2014; 10(9):e1003794.
 25
Moore H, Li NK. A mathematical model for chronic myelogenous leukemia (cml) and t cell interaction. J Theor Biol. 2004; 227:513–23.
 26
DeConde R, Kim PS, Levy D, Lee PP. Posttransplantation dynamics of the immune response to chronic myelogenous leukemia. J Theor Biol. 2005; 236:39–59.
 27
Stiehl T, Baran N, Ho AD, MarciniakCzochra A. Clonal selection and therapy resistance in acute leukemias: mathematical modelling explains different proliferation patterns at diagnosis and relapse. J R Soc Interface. 2014; 11.
 28
Takumi K, Garssen J, de Jonge R, de Jong W, Havelaar A. Release kinetics and cell trafficking in relation to bacterial growth explain the time course of blood neutrophils and monocytes during primary salmonella infection. Int Immunol. 2005; 17:85–93.
 29
LahozBeneytez J, Elemans M, Zhang Y, Ahmed R, Salam A, Block M, Niederalt C, Asquith B, Macallan D. Human neutrophil kinetics: modeling of stable isotope labeling data supports short blood neutrophil halflives. Blood. 2016; 127:3431–8.
 30
Craig M, Humphries AR, Mackey MC. A mathematical model of granulopoiesis incorporating the negative feedback dynamics and kinetics of gcsf/neutrophil binding and internalization. Bull Math Biol. 2016; 78(12):2304–57.
 31
Stiehl T, Ho A, MarciniakCzochra A. The impact of cd34+ cell dose on engraftment after scts: personalized estimates based on mathematical modeling. Bone Marrow Transplant. 2014; 49(1):30–7.
 32
Stiehl T, Ho AD, MarciniakCzochra A. Assessing hematopoietic (stem) cell behavior during regenerative pressure. In: A Systems Biology Approach to Blood. New York: Springer: 2014. p. 347–67.
 33
Ostby I, Rusten LS, Kvalheim G, Grottum P. A mathematical model for reconstitution of granulopoiesis after high dose chemotherapy with autologous stem cell transplantation. J Math Biol. 2003; 47:101–36.
 34
Østby I, Kvalheim G, Rusten LS, Grottum P. Mathematical modeling of granulocyte reconstitution after highdose chemotherapy with stem cell support: effect of posttransplant gcsf treatment. J Theor Biol. 2004; 231:69–83.
 35
Engel C, Scholz M, Loeffler M. A computational model of human granulopoiesis to simulate the hematotoxic effects of multicycle polychemotherapy. Blood. 2004; 104:2323–31.
 36
Stiehl T, MarciniakCzochra A. Characterization of stem cells using mathematical models of multistage cell lineages. ath Comput Model. 2011; 53(7):1505–17.
 37
Nakata Y, Getto P, MarciniakCzochra A, Alarcón T. Stability analysis of multicompartment models for cell production sytems. J Biol Dyn. 2012; 6:2–18.
 38
Getto P, MarciniakCzochra A, Nakata Y, dM. Vivanco M. Global dynamics of twocompartment models for cell production systems with regulatory mechanisms. Math Biosci. 2013; 245:258–68.
 39
Metcalf D. Clonal analysis of proliferation and differentiation of paired daughter cells: action of granulocytemacrophage colonystimulating factor on granulocytemacrophage precursors. Proc Natl Acad Sci. 1980; 77(9):5327–330.
 40
Smith CW. Production, distribution, and fate of neutrophils In: Kaushansky K, Lichtman MA, Prchal JT, Levi MM, Press OW, Burns LJ, Caligiuri M, editors. Williams Hematology. New York: McGrawHill: 2015. Chap. 61.
 41
Douglas SD, Douglas AG. Production, distribution, and activation of monocytes and macrophages In: Kaushansky K, Lichtman MA, Prchal JT, Levi MM, Press OW, Burns LJ, Caligiuri M, editors. Williams Hematology. New York: McGrawHill: 2015. Chap. 68.
 42
Seet CS, Crooks GM. Lymphopoiesis In: Kaushansky K, Lichtman MA, Prchal JT, Levi MM, Press OW, Burns LJ, Caligiuri M, editors. Williams Hematology. New York: McGrawHill: 2015. Chap. 74.
 43
Klonz A, Wonigeit K, Pabst R, Westermann J. The marginal blood pool of the rat contains not only granulocytes, but also lymphocytes, nkcells and monocytes: a second intravascular compartment, its cellular composition, adhesion molecule expression and interaction with the peripheral blood pool. Scandanavian J Immunol. 1996; 44:461–9.
 44
Zhang CC, Lodish HF. Cytokines regulating hematopoietic stem cell function. Curr Opin Hematol. 2008; 15:307–11.
 45
Gonda TJ, D’Andrea RJ. Activating mutations in cytokine receptors: Implications for receptor function and role in disease. Blood. 1997; 89(2):355–69.
 46
Mizuki M, Schw able J, Steura C, Choudhary C, Agrawal S, Sargin B, Steffen B, Matsumura I, Kanakura Y, B ohmer FD, M ullerTidow C, Berdel WE, Serve H. Suppression of myeloid transcription factors and induction of stat response genes by amlspecific flt3 mutations. Blood. 2003; 101:3164–73.
 47
Kohl TM, Ellwart SSJW, Hiddemann W, Spiekermann K. Kit exon 8 mutations associated with corebinding factor (cbf)–acute myeloid leukemia (aml) cause hyperactivation of the receptor in response to stem cell factor. Blood. 2005; 105:3319–21.
 48
Levine RL, Huntly MLBJP, Loh ML, Beran M, Stoffregen E, Berger R, Clark JJ, Willis SG, Nguyen KT, Flores NJ, Estey E, Gattermann N, Armstrong S, Look AT, Griffin JD, Bernard OA, Heinrich MC, Gilliland DG, Druker B, Deininger MWN. The jak2v617f activating mutation occurs in chronic myelomonocytic leukemia and acute myeloid leukemia, but not in acute lymphoblastic leukemia or chronic lymphocytic leukemia. Blood. 2005; 106:3377–9.
 49
Knudson AG. Mutation and cancer: statistical study of retinoblastoma. Proc Natl Acad Sci. 1971; 68(4):820–3.
 50
Rubinow SI, Lebowitz JL. A mathematical model of the acute myeloblastic leukemic state in man. Biophys J. 1976; 16:897–910.
 51
van den Akker E, Satchwell TJ, Pellegrin S, Daniels G, Toye AM. The majority of the in vitro erythroid expansion potential resides in cd34– cells, outweighing the contribution of cd34+ cells and significantly increasing the erythroblast yield from peripheral blood samples. Haematologica. 2010; 95.
 52
Stiehl T, Baran N, Ho AD, MarciniakCzochra A. Cell division patterns in acute myeloid leukemia stemlike cells determineclinical course:a model to predict patient survival. Cancer Res. 2015; 75:940–9.
 53
Shochat E, RomKedar V, Segel LA. Gcsf control of neutrophils dynamics in the blood. Bull Math Biol. 2007; 69:2299–338.
 54
Holländer GA, Widmer B, Burakoff SJ. loss of normal thymic repertoire selection and persistence of autoreactive t cells in graft vs host disease. J Immunol. 1994; 152:1609–17.
 55
Kim PS, Lee PP, Levy D. Minitransplants for chronic myelogenous leukemia: A modeling perspective. Biol Control Theory: Curr Challenges. 2007; 357:3–20.
 56
van Furth R, Dulk MMCDD, Mattie H. Quantitative study on the production and kinetics of mononuclear phagocytes during an acute inflammatory reaction. J Exp Med. 1973; 138(6):1314–30.
 57
Whitelaw DM, Batho HF. The disribution of monocytes in the rat. Cell Tissue Kinet. 1972; 5:215–25.
 58
Schmaier AH, Lazarus HM, (eds).Concise Guide to Hematology. West Sussex, UK: WileyBlackwell; 2012.
 59
Le T, Bhushan V, (eds).First Aid for the USMLE Step 1 2012. New York: McGrawHill Education; 2012.
 60
Morstyn G, Souza L, Keech J, Sheridan W, Campbell L, Alton N, Green M, Metcalf D, Fox R. Effect of granulocyte colony stimulating factor on neutropenia induced by cytotoxic chemotherapy. Lancet. 1988; 331(8587):667–72.
 61
Harlan JM. Leukocyteendothelial interactions. Blood. 1985; 65:513–25.
 62
Jedrzejczak WW. Mobilization of the marginal pool of neutrophils with epinephrine. results in healthy persons, patients with neutropenias, patients with neutrophilias, and patients with changes in neutrophil count induced by cancer chemotherapy. Haematologica. 1979; 64:586–96.
 63
Wäsch R, Reisser S, Hahn J, Bertz H, Engelhardt M, Kunzmann R, Veelken H, Holler E, Finke J. Rapid achievement of complete donor chimerism and low regimenrelated toxicity after reduced conditioning with fludarabine, carmustine, melphalan and allogeneic transplantation. Bone Marrow Transplant. 2000; 26:243–50.
 64
Jain AK, Murty MN, Flynn PJ. Data clustering: a review. ACM Comput Surv (CSUR). 1999; 13:264–323.
 65
Amari S, Wu S. Improving support vector machine classifiers by modifying kernal functions. Neural Netw. 1999; 12:783–9.
 66
Zheng Y, Rundell A. Comparative study of parameter sensitivity analyses of the tcractivated erkampk signaling pathway. IEEE Proc Syst Biol. 2006; 153:201–11.
 67
KinzerUrsem TL, Linderman JJ. Both ligand and cellspecific parameters control ligand agonism in a kinetic model of g protein–coupled receptor signaling. PLOS Comput Biol. 2007; 3(1):e6.
 68
Marino S, Hogue IB, Ray CJ, Kirschner DE. A methodology for performing global uncertainty and sensitivity analysis in systems biology. Theor Biol. 2008; 254:178–96.
 69
Krivtsov AV, Twomey D, Feng Z, Stubbs MC, Wang Y, Faber J, Levine JE, Wang J, Hahn WC, Gilliland DG, et al. Transformation from committed progenitor to leukaemia stem cell initiated by mll–af9. Nature. 2006; 442(7104):818–22.
 70
Haeno H, Levine RL, Gilliland DG, Michor F. A progenitor cell origin of myeloid malignancies. Proc Natl Acad Sci U S A. 2009; 106:16616–21.
 71
Zeijlemaker W, Gratama JW, Schuurhuis G. Tumor heterogeneity makes aml a “moving target” for detection of residual disease. Cytom Part B: Clin Cytom. 2014; 86(1):3–14.
 72
Döhner H, Estey EH, Amadori S, Appelbaum FR, Büchner T, Burnett AK, Dombret H, Fenaux P, Grimwade D, Larson RA, LoCoco F, Naoe T, Niederwieser D, Ossenkoppele GJ, Sanz MA, Sierra J, Tallman MS, Löwenberg B, Bloomfield CD. Diagnosis and management of acute myeloid leukemia in adults: recommendations from an international expert panel, on behalf of the european leukemianet. Blood. 2010; 115:453–74.
 73
Schiffer CA, Gurbuxani S, Larson RA, Rosmarin AG. Clinical manifestations, pathologic features, and diagnosis of acute myeloid leukemia. Waltham. http://www.uptodate.com/contents/clinicalmanifestationspathologicfeaturesanddiagnosisofacutemyeloidleukemia. Accessed 10 Jan 2017.
 74
Burnett AK, Milligan D, Prentice AG, Goldstone AH, McMullin MF, Hills RK, Wheatley K. A comparison of lowdose cytarabine and hydroxyurea with or without alltrans retinoic acid for acute myeloid leukemia and highrisk myelodysplastic syndrome in patients not considered fit for intensive treatment. Cancer. 2007; 109(6):1114–24.
 75
Ashkenazi R, Gentry SN, Jackson TL. Pathways to tumorigenesis—modeling mutation acquisition in stem cells and their progeny. Neoplasia. 2008; 10(11):1170–111826.
 76
Gentry S, Ashkenazi R, Jackson T. A maturitystructured mathematical model of mutation, acquisition in the absence of homeostatic regulation. Math Model Nat Phenom. 2009; 4(3):156–82.
 77
Stiehl T, MarciniakCzochra A. Mathematical modeling of leukemogenesis and cancer stem cell dynamics. Math Model Nat Phenom. 2012; 7:166–202.
 78
Cytarabine: Drug Information. Waltham. http://www.uptodate.com/contents/cytarabinepatientdruginformation. Accessed 10 Jan 2017.
 79
Metcalf D. Hematopoietic cytokines. Blood. 2008; 111(2):485–91.
 80
Maywald O, Buchheidt D, Bergmann J, Schoch C, Ludwig WD, Reiter A, Hastka J, Lengfelder E, Hehlmann R. Spontaneous remission in adult acute myeloid leukemia in association with systemic bacterial infection—case report and review of the literature. Ann Hematol. 2004; 83(3):189–94.
Acknowledgements
Not applicable.
Funding
Support for this research comes from the National Institutes of Health under Grant R01HD073156 and the Indiana CTSI TL1 Award, funded in part by Grant UL1TR001108.
Availability of data and materials
The clinical data that support the findings of this study are available from Dr. Robert P. Nelson, Jr., but restrictions apply to the availability of these data, which were used under IRB approval for the current study, and so are not publicly available. The remaining data are however available from the authors upon reasonable request.
Author information
Affiliations
Contributions
JMS and SR created the multilineage model. JMS completed the model analysis and drafted the written work. RPN contributed to the conception, interpretation, and revision of the work. AER, DMU, and TKU contributed to the study design, interpretation, and revision of the written work. All authors read and approved the final manuscript.
Corresponding author
Ethics declarations
Ethics approval and consent to participate
Deidentified clinical data was used under Purdue University’s IRB Protocol# 1011009928 from the Indiana University Simon Cancer Center from patients who received hematopoietic stem cell transplants during a study approved by the IRB of Indiana UniversityPurdue University Indianapolis.
Consent for publication
Not applicable.
Competing interests
The authors declare that they have no competing interests.
Publisher’s Note
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Additional file
Additional file 1
Supplementary information is provided as pdf files. Supplementary Information 1. Existing Mathematical Models. Current models that describe a subset of mechanisms described in our model are compared in this section. Supplementary Information 2. Model Equations. All ODE model equations are included, as well as diagrams of the positive and negative feedback included in the model. Supplementary Information 3. Detailed Model Description. The model descriptions explain the mechanisms included in the model. A derivation of the form of asymmetrical selfrenewal, the capacity of the bone marrow, and the pharmacodynamics of chemotherapy are also detailed in this section. Supplementary Information 4. Constraint and Parameter Tables. Constraint and parameter value ranges, descriptions, and citations are described in this supplement. Supplementary Information 5. Physiological Capabilities of Model. Qualitative patient data are shown to demonstrate validation of dynamical acceptability criteria. Supplementary Information 6. Sensitivity Analysis. Results from completing PRCC are shown. Supplementary Information 7. Determining Cluster Number. A derivation of the optimal number of clusters to describe the patient variation is explained in this supplement. Supplementary Information 8. NeutrophilDerived AML. Simulations of AML derived from neutrophil progenitors. (PDF 1502 kb)
Rights and permissions
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License(http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.
About this article
Cite this article
Sarker, J., Pearce, S., Nelson, R. et al. An Integrative multilineage model of variation in leukopoiesis and acute myelogenous leukemia. BMC Syst Biol 11, 78 (2017). https://0doiorg.brum.beds.ac.uk/10.1186/s1291801704692
Received:
Accepted:
Published:
DOI: https://0doiorg.brum.beds.ac.uk/10.1186/s1291801704692
Keywords
 Acute myelogenous leukemia
 Mathematical model
 Personalized medicine
 Hematopoiesis
 Leukopoiesis