An Integrative multi-lineage model of variation in leukopoiesis and acute myelogenous leukemia

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 multi-lineage multi-compartment model of the hematopoietic system that captures patient-to-patient 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 self-renewal of progenitor cells are inversely related. Within the data-consistent 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 k-means 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 multi-lineage model of hematopoiesis, we have identified several key features that determine leukocyte homeostasis, including self-renewal 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 self-renewal 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 self-renewal of cancer stem cells; most current therapies target mitosis rate, but based on our results, we propose that additional therapeutic targeting of self-renewal of cancer stem cells will lead to even higher survival rates. Electronic supplementary material The online version of this article (doi:10.1186/s12918-017-0469-2) contains supplementary material, which is available to authorized users.


Supplementary Information 2: Model Equations
: Parameter Set Categories. Seven parameters were determined from the literature. The remainder of the undetermined parameters were constrained using a step-wise approach described in the paper text and Figure 2 of the text. The parameters in the model equations below are highlighted with colors according to the key in this  Figure S1: Negative Feedback Signals of Model. The figure notations are the same as in Figure 1 of the text with the addition of inhibition depicted by red bar-headed arrows. The feedback is described in the methods section of the paper.The numbered labels of the feedback are correspond to the numbered constants described in Table S4 Figure S2: Positive Feedback Signals of Model. The figure notations are the same as in Figure  1 of the text with the addition of positive feedback depicted by green dashed arrows. The feedback is described in the methods section of the paper.The numbered labels of the feedback are correspond to the numbered constants described in Table S4  Table S4: Description of All Feedback Processes. Each of the feedback processes in our model are described in this table with citations for justification of the feedback process. Signals (cytokines, chemokines, transcription factors, etc.) that approximate these processes are shown in this table. Other factors not listed in this table may also regulate these processes. M-M := Michaelis-Menten; TF := transcription factor; PU.1 and C/EBP α := transcription factors that regulate at the multi-potent progenitor level; EPO:= erythropoietin; GM-CSF : = granulocyte macrophage colony stimulating factor; IL := interleukin; LFA-1 := lymphocyte function-associated antigen 1; SDF-1 := stromal cell-derived factor-1; CXC := chemokine motif "C-X-C"; CC := chemokine motif "C-C" ; LPS := lipopolysaccharide; TGF-β := transforming growth factor; TARC : = T cell-directed CC chemokine thymus and activation-regulated chemokine;

Neutrophils
2. Rate of change for committed neutrophils, N2 3. Rate of change for mature neutrophils, N3 4. Rate of change for mature neutrophils in peripheral blood, N3pb

Monocytes
9. Rate of change for committed monocytes, M2 10. Rate of change for mature monocytes, M3 11. Rate of change for mature monocytes in peripheral blood, M3pb 12. Rate of change for mature monocytes in marginal pool, M3mp

Rate of change of cellular debris and apoptotic cells, A
[6] A. Klonz, K. Wonigeit, R. Pabst, and J. Westermann. The marginal blood pool of the rat contains not only granulocytes, but also lymphocytes, nk-cells and monocytes: a second intravascular compartment, its cellular composition, adhesion molecule expression and interaction with the peripheral blood pool. Scandanavian Journal of Immunology, 44:461-9, 1996.
[7] Alessandro Aiuti, IJ Webb, C Bleul, T Springer, and JC Gutierrez-Ramos. The chemokine sdf-1 is a chemoattractant for human cd34+ hematopoietic progenitor cells and provides a new mechanism to explain the mobilization of cd34+ progenitors to peripheral blood. Journal of Experimental Medicine, 185(1):111-120, 1997. The justification for each of the eight processes included in our model is explained in the following table. Represents the majority of white blood cell count and provides a mechanism of G-CSF control [1]; typically largest white blood cell population [2,3]; easy to validate due to being commonly measured Lymphocytes Second most common white blood cell [3]; potential to expand the model for a mechanism for graft-versus-leukemia or graft-versus-host-disease Differentiates between latent mature leukocytes and activated leukocytes [1,[9][10][11][12]

Marginal Pool
Demonstrates rapid demargination of neutrophils and monocytes into peripheral blood after initial chemotherapy or corticosteroid doses [14] Leukemia Represents the leukemic variability in patients; Validates the effectiveness of the treatment variables on patients 2 Self-Renewal Self-renewal and differentiation rates are modeled using the model proposed by Marciniak-Czochra [15]. They are detailed below. Figure S3: Stem cell self-renewal and differentiation.
A stem cell, SC, has the potential to self-renew or differentiate into a mature cell, N , as shown above in Figure S3. When stem cells make the decision to create two daughter cells, a quasi-equilibrium state is maintained, D. These daughter cells have the option of remaining as stem cells or becoming mature cells at some rate. Thus, when excluding terms other than the self-renewal and division rates, the mitosis term can be simplified as such: Quasi-equilibrium is maintained in the daughter state, so setting Equation 2 to 0, we can solve to get where α is a term that represents the fraction of cells that will self-renew when undergoing mitosis. Thus, we can simplify Equations 1 and 3 to In this formulation, α must remain between 0.5 and 1 in order to ensure that dSC dt and dN dt do not go negative due to the mitosis term (Equation 8). For the remainder of the paper, we refer to k d iv as a mitosis rate, mr.
As HSCs progress in their differentiation, they lose their ability to self-renew. In Marciniak-Czochra et al. [16], this theory is maintained by ensuring that α, the fraction of cells that self-renew gets smaller as cells mature. Thus, progenitor cells will have smaller αs than stem cells (Equation 9). Additionally, the mitosis rates of stem cells are very small, and the mitosis rates of more mature cells should be faster (Equation 10).
Stem cells can differentiate into many types of progenitors. In order to maintain mass balance, we have modeled that a fraction of these cells become lymphocyte progenitors, f L , a fraction of the non lymphocytes become monocyte progenitors, (1−f L )f M , and the remainder of the cells become neutrophil progenitors, (1 − f L )(1 − f M ).These progenitor cells have the capability of producing mature cells of their specific lineages [17,18].

Homeostatic Feedback of Stem Cells: Capacity of Bone Marrow
To ensure that the concentration of the marrow cells do not exceed the capacity of the marrow, we include a homeostatic feedback term to the stem cells regeneration capability. This term is a Hill term with a constant (k 1 ) that modulates the concentration of cells in the bone marrow (H) and the number of stem cells in the bone marrow (SC 1 ), divided by constants (n bm and n sc ) to bring these terms to the same scale. We assume that the maximum concentration of cells in the bone marrow is the concentration of a single cell over its own volume. We assume that the concentration of cells in the bone marrow will not exceed the concentration of one cell per volume of one normal monocyte with a diameter of 20 µm [19], as monocytes are the largest inactivated type of white blood cell in our model.
Thus, we limit n bm to between the concentration of one red blood cell over its own volume, with a diameter of 7 µm [20] and the concentration of one monocyte over its own volume because we assume all of the white blood cells are between the size of a red blood cell and a monocyte. We calculate the approximate maximum concentration of cells in the bone marrow to be, converted to cells per microliter (Equation 12): 1 monocyte volume of 1 monocyte = 1 cell 4.18 × 10 3 µm 3 × 1µm 3 1 × 10 −9 µL = 2.39 × 10 5 cells/µL Thus, we started our search for n bm at 3 × 10 5 , with a maximum of 5.5 × 10 6 , if the cells have the diameter of a red blood cell. For simplicity, we set our maximum concentration limit to 3 × 10 5 cells/µL to characterize patient death.

Chemotherapy
In order to elucidate the variation in dynamics of individual patients in response to input and to attain the range of parameters that lead to these dynamics, we incorporate a general term for cell loss due to chemotherapy. In this mathematical model, chemotherapy causes death of the stem cells and progenitor cells at a value that mimics pharmacodynamic (PD) effects. The overall PD effect is modeled to be a scalar greater than zero. This scalar is multiplied to a general decay term for all cells that are affected by chemotherapy, and these cells are moved into the apoptotic state, which are cleared by macrophages, as described in Supplement 2.
The PD effect of one treatment cycle is modeled as a combination of drugs that are administered daily for seven days. The drug administration is a product of two exponentials that takes into account both the response time, t r , and the duration of action, t d (Equation 13).
where P D i is the pharmacodynamic effect of one drug, i, and t i is the time of drug administration for drug i, w i is a constant to modulate the effect of drug i, and dose is the dose of drug i. The effects of the administration of each drug are summed together for the overall PD effect of the drugs. For the purpose of simulating dynamics in virtual 'patients' with unique parameter sets, we simulate the effects of the combination of drugs, cyclophosphamide for the first two days, and fludarabine for the 5 subsequent days. The t r and t d for cyclophosphamide are 0.1/day and 0.31/day, respectively [21], and the t r and t d for fludarabine are 0.125/day [22] and 0.54/day [23], respectively. Thus, a stronger pharmacodynamic effect occurs for the first two days of chemotherapy compared to the last five days of chemotherapy. We do not fit for chemotherapeutic precision, so we assume that the w i × dose i for each drug is 1. For the purposes of this paper, chemotherapy is applied on day -7 through day -1.

Supplementary Information 4: Constraint and Parameter Tables 1 Model Validation
In order to ensure that our model maintained normal homeostatic cell concentrations, we set one of our healthy dynamic criteria to fit concentrations found from literature (below).

Parameter Tables
In this section, the values of the parameters and a description of their mechanism are shown. First, the values of the parameters set from literature are shown in Table S7. Next, the general range (25th percentile to 75th percentile) of all parameter values that were sampled using Latin Hypercube Sampling are shown, along with the median value (Table S8. Finally, a comparison of the parameters that are analagous to each other for each of the uni-lineages are shown ( Figure S4).   Table S9 shows parameters that dictate the ability of macrophages to clean up apoptotic debris and to recruit other cells to assist in clearing or attacking pathogens. For those parameters not found in literature, these parameters were nominally set to illicit minimal responses from macrophages. However, this aspect of the model can be explored further for studies that characterize macrophage behavior.  Figure S4. The self-renewal probability of stem cells is much higher than progenitor cells, and the mitosis rate of stem cells is lower than that of the other progenitor cells, as expected. Additionally, the recruitment rate of monocytes is much larger than the recruitment rate of the other cells, potentially due to the mechanism of forming monocytes for clearing apoptotic debris. Finally, the Michaelis-Menten constants associated with neutrophils and lymphocytes are much larger than that of monocytes, potentially due to the much larger concentrations that exist for neutrophils and lymphocytes than for monocytes. Figure S4: Comparison of Analagous Parameter Ranges from Uni-lineage Models. Abbreviations for these parameters can be found in Table S2  To identify correlations amongst the various lineages, both in parameter space and initial condition (I.C.) space, we computed correlation coefficients for all hematopoieitic non-cancer states and param-eters in our model of all multi-lineage acceptable simulations ( Figure S5). As a reminder, we found initial condition and parameter values simultaneously. We found that strong correlations exist between initial conditions and parameters of analogous lineages, and weak correlations exist between different lineages. However, stem cell parameters are slightly correlated with the initial condition of stem cells and neutrophils in the bone marrow. This is especially true of k 1 , the Michaelis-Menten constant that drives stem cells to maintain homeostasis of the whole system. Figure S5: Correlation Coefficient Heatmap of Initial Conditions and Parameters. The row and column show initial conditions (I.C.) and parameters (Param.) in the same order. The initial conditions of the states are in the same order of the first 12 equations in Supplement 2. Namely, the stem cell initial condition is followed by neutrophil (Neut), lymphocyte (Lymp), and monocyte (Mono) initial conditions. Within each category of initial conditions, the progenitor, mature cell in the marrow, the mature cell in the peripheral blood, and the mature cell in the marginal pool (if applicable) are shown. After the initial conditions are compared, the next values are parameters in the order presented in Table S7. White boxes group all initial conditions of the same cell types together. Black boxes group all parameters of the same cell types together. Finally, dark purple boxes group initial conditions and parameters of the same cell type together. The color bar scale for the value of the correlation coefficient is displayed on the right.

Supplementary Information 5: Physiological Capabilities of Model 1 Summary of Patient Data
Forty-seven adult patients who underwent peripheral blood stem cell transplants at the IU Simon Cancer Center under the supervision of Dr. Robert P. Nelson, Jr. for the treatment of AML were used to determine patient dynamics in response to chemotherapy. Peripheral blood samples were collected from the patients daily for seven days prior to transplantation, and for thirty days post-transplantation.
A complete blood count of the peripheral blood demonstrated neutrophil, lymphocyte, and monocyte concentrations in the peripheral blood. These patients received a non-myeloablative treatment protocol consisting of cyclophosphamide for two days, and fludarabine for the 5 subsequent days. The donor infusion is administered after the chemotherapy regimen is completed. We do not simulate transplantation within our model. We use the daily blood counts of patients to characterize cell dynamics in the peripheral blood during chemotherapy and qualitative characteristics of patient recovery after chemotherapy. We use transplantation data because more discrete data is available for patients after treatment for a longer time period. The real patients have received transplantation from donors, but we assume that the behaviors of the final cell concentrations would be the same between patients undergoing only chemotherapy and those undergoing transplantation in our simulations. We do not characterize the timed dynamics of patient recovery from transplantation or chemotherapy in the context of the current model.

Overshoot
In dynamic acceptability criteria #6, the overshoot of recovery post-chemotherapy in patients was required to be less than 12 times the overshoot value five days post-overshoot. Some simulations displayed overshoot more than we expected, so we wanted to constrain to normal physiological overshoot. As shown in Figure S5, the only one patient demonstrated overshoot of 12 times in the neutrophil lineage. The majority of patients had less than a 6-factor overshoot. Thus, we restricted our criteria to a factor of 12 change for overshoot.
In observing the final multi-lineage simulations ( Figure 6 in text, left column, blue lines), the majority of simulations do not overshoot. For those simulations that do overshoot, the overshoot is within a factor of 3 change. This indicates that though a overshoot of 12 times its equilibrium value was used to constrain the parameter space, a smaller overshoot criteria would have been acceptable. To discriminate for each lineage qualitatively in the patient data, only a small number of patients demonstrate neutrophil overshoot though neutrophils have the potential to overshoot the most. Our multi-lineage model analogously shows the smallest amount of overshoot amongst the lineages, and in the uni-lineage model, neutrophils have the potential to overshoot the most (Figure 6b). Additionally, a small percentage of patients demonstrate lymphocyte and/or monocyte overshoot, which we also show in our model ( Figure  6c and d).

Dampened Oscillations
Simulations that are manually classified as sufficiently dampened versus not dampened are shown in Figure S6. Here, the simulation that does not have a dampened oscillation ( Figure S6a) has period of about 60 days, which is longer than physiologically expected for a patient recovering from chemotherapy. On the other hand, a simulation that we manually classifed as acceptable ( Figure S6b) has a very high frequency oscillation period after recovering from chemotherapy. After manually classifying 100 Figure S5: Histogram of Patient Overshoot Distribution. The overshoot was calculated for the three lineages of data from 47 patients. Neutrophil overshoot is in blue; lymphocyte overshoot is in green; monocyte overshoot is in yellow.
simulations as physiologically acceptable, we found that dynamic criteria #7 uniquely allocated the manually classified simulations correctly. We analyzed the robustness of our model by testing the difference in decoupling the lineages and perturbing one state. We used the 13 centroids of our clustering solutions to change the initial condition of neutrophil progenitors to 10 times its initial value. The results are shown in Figure S7. In the multi-lineage model, this tempers the overshoot of the neutrophil concentrations in the peripheral blood ( Figure S7a) but causes oscillations in lymphocyte concentrations in the peripheral blood ( Figure S7b). However, in the uni-lineage model, the neutrophils demonstrate much higher overshoots ( Figure S7c), and there are no changes in the lymphocyte concentrations ( Figure S7d). We do not show monocytes because the changes in their concentrations are minimal with respect to their concentration. This confirms that combining all three lineages provides some robustness to the system.

Supplementary Information 7: Determining Cluster Number
K-means clustering requires that the number of clusters be pre-defined before computing the clusters. However, determining the number of clusters to group data is a complex problem. One method is to use the Calinski-Harabasz criteria, or F-statistic, to determine the number of clusters that maximizes the between-cluster variance, but minimizes the within-cluster variance [1]. The Calinski-Harabasz criterion is defined as where k is the number of clusters, i is the cluster index, n i is the number of elements in cluster i, m i is the centroid of cluster i, m is the overall mean of the data, N is the number of samples in the data, and x is a data point in the ith cluster, representated by c i . As the number of clusters increases, the within-cluster variance increases, and in order to ensure that the smallest k is chosen that still produces an acceptable within-cluster variance, a curve can be plotted of the F-statistic against the number of clusters. The optimal cluster number can be found at the elbow of this curve. The elbow of the curve is found by first drawing a line between the two ends of the curve. The point that has the maximum distance between the line and the curve defines k, the optimal cluster number. We wanted to ensure that at least one hundred representatives were in each cluster, so we clustered the data over one cluster to 228 clusters (for the 22,796 representatives; Figure S12). Due to computational time, some of the clusters were not computed. The optimal number of clusters, k, was determined to be 13 ( Figure S12). Figure S12: Determining Optimal Cluster Number.The F-statistic is maximized as the number of clusters increases (blue curve). The elbow of the curve determines the optimal number of clusters, which converged to 8 clusters (red +).

Supplementary Information 8: Neutrophil-Derived AML
Neutrophils and monocytes share a common progenitor. AML can arise when the feedback mechanisms of proliferation and/or differentiation are disrupted in a common myeloid progenitor or a granulocytemonocyte progenitor [1]. Thus to demonstrate the differences between AML derived from neutrophil progenitors versus those derived from monocyte progenitors, we modeled AML growth from neutrophil progenitors that do not respond to any feedback signals. This is analagous to how we modeled AML derived from monocytes. We find that AML derived from neutrophils grows much more rapidly and is more likely to lead to faster death in patients who have neutrophil-derived AML ( Figure S13a-c). This is probably due to the faster mitosis rate of neutrophil progenitors than monocyte progenitors that we determined in Table S7 in Supplement 4. Additionally, due to signaling of homeostastic mechanisms to stem cells and normal progenitors, the cell counts of normal cells remain lower in a neutrophil-derived AML, but the overall white blood cell counts would be higher earlier in these patients' time courses. For one subtype of AML, overall high white blood counts are predictive of survival [2], which is what these results indicate.
However, we confirm that the mitosis rate of the AML is the most correlated parameter to time to death for patients with AML derived from neutrophil progenitors with an overall correlation coefficient of -0.64 ( Figure S13d). Additionally, the probability of self-renewal is the second most correlated parameter to survival time of patients with neutrophil-derived AML with an overall correlation coeffeicient of -0.36 ( Figure S13d). This corresponds to the same results we found for monocyte-derived AML, and confirms the potential target of treating AML as the probability of self-renewal of progenitor cells. Figure S13: Survival in neutrophil-derived leukopoesis. Each of the thirteen clusters is represented with a distinct color in a-c, corresponding to the same colors as in Figure 7 in the body of the paper. (a) Progenitor cancer stem cells derived from neutrophil progenitors grow over one year in each of the 22,796 simulations. The gray lines represents the cell concentration and time in which patients may die due to overgrowth of cancer cells. (b) The percentage of cells in the peripheral blood due to AML in each of the 22,796 simulations. Patients die prior to reaching 100% peripheral blasts, as indicated by gray lines where simulation continued. In (a) and (b), some cluster simulations are not visible because they are underneath the other simulations. (c) The simulated survival curves of each of the thirteen representative clusters of patients are shown assuming the patient dies when either cancer in the bone marrow or the peripheral blood reaches a concentration greater than 3 × 10 5 cells/µL. The cumulative survival is shown as a decimal that represents the fraction of patients that are still alive at the time shown on the x-axis.(d) Significant cancer parameter correlation to number of days until death due to bone marrow or peripheral blood density greater than 3 × 10 5 (p < 0.05). Gray boxes represent insignificant correlations, and the correlation values for each of the thirteen clusters are depicted with the color bar.