RoskildeUniversity Dose-dependent mathematical modeling of interferon--treatment for personalized treatment of myeloproliferative neoplasms

Long-term treatment with interferon-alfa (IFN) can reduce the disease burden of patients diagnosed with myeloproliferative neoplasms (MPNs). Determining individual patient responses to IFN therapy may allow for efficient personalized treatment, reducing both drop-out and disease burden. A mathematical model describing hematopoietic stem cells and the immune system is suggested. Considering the bone marrow and the blood allows for modeling disease dynamics both in the absence and presence of IFN treatment. Through comprehensive modeling of the effects of IFN, the model was related to individualized patient-data consisting of longitudinal hematologic and molecular measurements. Treatment responses were modeled on a population level, allowing for personalized predictions from a single pretreatment data point. Personalized fits were found to agree well with data for individual patients. This allowed for a quantitative description of the treatment response, yielding a mechanistic interpretation of differences from patient to patient. The treatment responses of individual patients were combined and a formulation of treatment responses on the population level was described and simulated. Based on pretreatment data and the

modeling of treatment responses from pretreatment data successfully predicted clinical measures for a 5-year period. We believe that this approach could have direct clinical relevance, offering expert guidance for clinical decisions about IFN treatment of MPN patients.

K E Y W O R D S
interferon, mathematical modeling, myeloproliferative neoplasms, personalized treatment INTRODUCTION The Philadelphia-negative myeloproliferative neoplasms (MPNs) are blood cancers, which include the diseases essential thrombocythemia (ET), polycythemia vera (PV), and primary myelofibrosis (PMF) [2,35]. These MPNs develop in a biological continuum over decades from the early cancer stages (ET,PV) to the advanced myelofibrosis stage and ultimately leukemic transformation (acute myeloid leukemia) in a subset of patients. Thromboembolic complications are common, and MPNs are associated with a huge inflammation-mediated comorbidity burden. The MPNs are driven by the so-called "driver" mutations, which include the JAK2 V617F (JAK2) mutation (present in virtually all patients with PV and half of ET and PMF patients), the CALR mutations and MPL mutations. Additional mutations are determinant for disease progression with myelofibrotic and leukemic transformation [34]. Recent studies show that these MPNs evolve from an early precursor MPN stage-clonal hematopoiesis of indeterminate potential (CHIP)-in which the JAK2 allele burden steadily increases to give rise to overt MPN-disease [17]. In this work, we adopt the notion of a biological continuum closely related to the JAK2 allele burden, and consequently focus only on ET, PV and PMF patients that are positive for the JAK2 mutation. Chronic inflammation is considered the driving force for development of MPNs from the CHIP-stage but also from early to advanced MPN disease [14,16]. Another hallmark of MPNs is a disturbed immune system with ensuing defective tumor immune surveillance, which likely accounts for the increased risk of second cancers, not only after the diagnosis but also years before the diagnosis [13,15]. Improved understanding of the link between disease progression, chronic inflammation and the immune system could improve future treatment by halting the disease before the vicious cycle of chronic inflammation becomes uncontrollable. Therapy with pegylated interferon-alpha (IFN) was investigated thoroughly in a recent clinical trial, showing that IFN causes an exponential decline in JAK2 allele burden during long-term therapy [28]. Despite well-documented effects of IFN with depletion of MPN-propagating stem cells in a murine model [23], and a history of IFN therapy for MPN patients [33], IFN is not yet considered the standard of care for MPN patients.
While validation and testing of predictive power of mathematical modeling is necessary before clinical application [5], mathematical modeling is apt to become an important part of clinical assessment of disease stage and decision-making, for both blood cancers and solid tumors [8].
Mathematical models of hematopoietic stem cells (HSC) have given important insight about the hematopoietic system [4,6,11,20,29,30]. Myeloid hematologic cancers arise from a malignant stem cell clone, leading many authors to explicitly model the malignant clone and the competition of stem cells [3,31,38]. Importantly, modeling suggests that treatment of HSC-derived disease must target the malignant stem cells [1,10,24,25]. Consequently, mathematical models of myeloid hematologic malignancies should consider some notion of HSC if successful therapy is expected to be captured by model behavior.
In a recent study, we mathematically modeled the IFN-induced decline in disease burden observed in MPN patients, predicting whether a state of minimal residual disease would be attained [24]. The model described HSC, blood-cell production, and immune-system feedback and was originally calibrated to data from MPN patients [1].
In this work, we propose a related model in which the description of HSC behavior has been expanded. The proposed model is related to both hematologic and molecular measures, demonstrating that both cell counts and the JAK2 allele burden of MPN patients can be modeled simultaneously. We first describe the model, and how IFN therapy was interpreted mathematically. A fittingprocedure relating the model to data of individual patients is described. Personalized modeling suggests that careful consideration of molecular measures is necessary to determine whether treatment cessation will cause an early relapse or not. Population-level model predictions are made for a subcohort of patients, in which the full period of IFN treatment (up to 5 years) is predicted from a single pretreatment data point. This proof-of-concept result suggests that mathematical modeling on both individual and population level could be a useful clinical tool in the near future, predicting, for example, the responses to treatment cessation or an increase in dosage.

Mechanism-based mathematical model of the hematopoietic system
To model the effect of IFN on the hematologic and molecular variables of patients diagnosed with MPN, we first introduce a mathematical model of the interplay of the healthy hematopoietic stem cells (HSC), leukemic stem cells (LSC), blood cells, and the immune cells. The model is based on two previously presented models: the Cancitis model [1,24] and a model of HSC dynamics within the bone marrow [30].
The model is given by a system of ordinary differential equations describing the amount of wild-type HSC denoted , the amount of LSC denoted , the number of mature blood cells arising from the wild-type HSC denoted by , the number of mature blood cells arising from the LSC denoted by , the amount of cellular debris denoted by , and an abstract measure of the immune system response denoted by . We give a brief overview of the concepts behind the model formulation, while additional details are described in the Supporting Information. Dynamics of wild-type HSC and LSC are based on the model of Pedersen et al. [30]. In the previous model, a production of differentiated stem cells were described but not modeled further. In the model presented here, the differentiated stem cells give rise to a production of blood cells. The mature blood cells undergo apoptosis (programmed cell death), which contributes to the cellular debris ( ). As in the model presented in the work of [1], the debris upregulates the immune system ( ), which in turn increases stem cell production, to maintain cell counts. This process is thought to ensure homeostasis within the hematopoietic system and results in a model with both diseasefree hematopoiesis, early-disease stages, and treatment. A schematic diagram of the model is shown in Figure 1.
The model is given as the following system of ordinary differential equations, describing the gain and loss of each of the six variables: F I G U R E 1 Schematic compartment diagram of the mathematical model wild-type HSC and the produced blood cells, and , respectively, are shown as blue circles on the left, while LSC and the LSC-derived blood cells, and , respectively, are shown as red circles on the right. The cellular debris arising from apoptotic cell death of blood cells, , is shown as a gray circle, while the immune system, , is represented by a gray box in the middle. Black arrows represent flows between compartments, and red arrows signify upregulation by the immune system. The circles with × represent a multiplication in numbers to the proliferation of progenitor cells. The self-renewal of HSC and LSC is not depicted in the figure due to its nonlinear form; however, the gray box and the line connecting HSC and LSC represent this interaction of stem cells within the bone-marrow niches (shown as a dark-red box at the bottom)̇= − , (1c) where = (2 − 2 + 2 parameters are nonnegative. Additionally, ≤ 1 and ≥ 1. Default values given in Table 1, along with a brief description of what biological processes the parameters represent. It can be shown that all variables remain nonnegative for nonnegative initial conditions. A numerical solution, representing a typical disease-progression scenario is depicted in Figure 2, together with simulated treatment, as described in Section 2.2. Numerical solutions like these suggest that the model is qualitatively similar to those of the Cancitis model on which it is based. Furthermore, model dynamics are similar to those TA B L E 1 Default parameters used in simulations. Parameters , , , , , , and are in units of [days −1 ], while the other parameters are without unit. Parameters , , , and were determined from the healthy counterparts and were modified to agree with the disease progression as described in [28]. The remaining default values were determined in the work described in [24] and [29]  Clearance rate of cellular debris 0.0003 Debris-dependent immune-system activation 2 Immune-system inactivation 7 External inflammation F I G U R E 2 Illustrative model simulation with treatment affecting , , and . Starting in the healthy steady state at time = 0, one HSC is removed from and added to . Due to the growth advantage of the leukemic clone for default parameters, the disease progresses in the absence of treatment. This progression is shown in black in both panels. The relative frequency of malignant mature cells, observed for similar models of the hematopoietic system seen in the literature; examples of which were mentioned in the introduction. The model considers positive feedback of the inflammatory stimulus on the release of stem cells from the bone-marrow niche. This is based on the idea of inflammation-induced activation of quiescent stem cells (see Chapter A.2 of the Supporting Information). The parameters and relate to the rate of this activation for wild-type HSC and LSC, respectively, allowing for a difference in how strongly inflammation affects either population. For the default parameters given in Table 1, > , consistent with increased inflammationinduced activation of LSC compared to HSC. From the local parameter sensitivity analysis given in the Supporting Information, increasing leads to increased counts of wild-type mature cells , while increased accelerates disease progression and increases . We define the disease level in the model as the relative frequency of blood cells derived from LSC out of all blood cells, For parameters close to the default values given in Table 1, the blood cell count in the model is of the order of 10 12 . The blood cell count represents the entire population of blood cells of a person, regardless of cell type. We emphasize that a value of 10 12 is arbitrary, and that the total blood cell count could vary from patient to patient. For this work, the relative difference from healthy cell counts or the cell counts of particular cell types is more relevant in regard to the health of a patient. To consider any particular cell type, a fraction of the total blood cell count in the model is assumed to be of the given cell type. We define the thrombocyte count as ℎ = ℎ ( + ) and the leukocyte count as = ( + ), where ℎ and are the cell-specific fractions.
Mathematical analysis of the model was described in [27]. We omit a thorough mathematical analysis here; however some brief comments are warranted. For most choices of model parameters, three steady states exists: A trivial steady state with no cells, a healthy steady state with = = 0, and a full-blown myeloproliferative steady state with = = 0. In the full-blown myeloproliferative steady state, only mature cells that have arisen from LSC are present, and the disease level is exactly 1. As such, the steady state represents the most extreme progression of disease, at which the patient is at high risk of life-threatening adverse events. Numerical investigations of local stability reveal that typically only one steady state is stable, while the two others are unstable. This suggests that for a given choice of parameters the system will asymptotically approach the stable steady state. Within the range of model parameters that we here investigate numerically, model solutions (and hence patient trajectories) move from either the healthy steady state toward the full-blown myeloproliferative steady state, or in the opposite direction. Stability of the myeloproliferative steady state is associated with increased cell counts in the model. Hence, solutions with disease progression exhibit increasing ℎ and toward a steady-state count and the disease level increasing from 0 to 1. The opposite is the case when the healthy steady state is stable: cell counts are decreasing toward a steady-state count and the disease level goes toward 0. To investigate the transient behavior during treatment, the model can be solved numerically, simulating particular treatment effects on parameters.

Modeling treatment effects
Treatment with IFN has been found to deplete the population of LSC that give rise to MPN [23]. To model this effect, we consider an IFN-induced decrease of , the parameter related to LSC differentiation. Preliminary model investigations revealed that decreasing resulted in depletion of LSC and a decrease in the disease level. Visual inspection of data suggests a transient effect of IFN on the total cell count, reducing both healthy and leukemic mature cells before the disease level decreases. A change in the proliferation of both healthy and leukemic progenitor cells, as modeled by parameters and , respectively, could explain such a decrease in total cell count. Hence, when modeling IFN treatment, we perturb parameters , , and . To reduce the degrees of freedom, we consider only equal perturbation of and , such that the relation = 2.5 is maintained throughout this work for both default parameter values and perturbed parameters. Perturbation of affects primarily the relative frequency of cells, while perturbation of and affects the total amount of mature blood cells. A numerical example of perturbing these parameters is shown in Figure 2. Within the first year of treatment, a transient effect of perturbing and causes a pronounced decrease in blood cell counts. A corresponding transient increase is seen when treatment is ceased. This transient effect occurs because and affects both healthy and leukemic blood cells, and as a consequence, without perturbation of and , the transient decrease of total blood cell counts does not occur (not shown).
For default parameters, solutions of the model approach the full-blown myeloproliferative steady state for any positive number of malignant stem cells, and, hence, relapse always occurs. However, the degrees of parameter perturbation considered in the simulated treatment affects the time until relapse occurs. To estimate the time between treatment cessation and relapse, we defined an arbitrary relapse threshold corresponding to a 50% increase in blood When total blood cell counts exceed this threshold after treatment cessation, relapse is assumed to have occurred. Note that this arbitrary threshold is a simplification compared to diagnostic criteria for relapse. For the example shown in Figure 2, relapse occurs 6 years after treatment cessation, around year 28. Simulating a range of treatment-specific values of , , and , we investigated the resulting time to relapse. The results are shown in Figure 3. While changes to and are important for a transient reduction in blood cell counts, the effect on the time to relapse is minor. Conversely, reduction of lengthens the time to relapse.

Data
We consider data from the prospective randomizedcontrolled open-label phase III clinical trial "DALIAH" (EudraCT number: 2011-001919-31) [24,28]. In the trial, a cohort of MPN patients received IFN monotherapy (either interferon alfa-2a "Pegasys R " or interferon alfa-2b "PegIntron R "). Data consisted of IFN dosage and timing, as well as longitudinal hematologic and molecular measurements. In particular, the thrombocyte and leukocyte counts as well as the JAK2 allele burden were measured. We consider the JAK2 allele burden as a measure of the disease progression. Here, a total of 63 patients [24] are included in the analysis: 17 were diagnosed with ET, 35 with PV, and 11 with PMF.

Pharmaco-kinetic modeling of IFN-dose
IFN was given once every week, once every second week, or once every third week in a range of dosages. For simplicity, we here report a daily average dose. A pharmacokinetic (PK) model of IFN was considered, in agreement with our previous work [24]. We assume an equal rate of uptake and clearance of IFN, = 1 7 g day , inspired by previous work on PK modeling of IFN [32]. Letting ( ) denotes the stepwise constant daily dose in units of g IFN, we choose the simplest possible PK description of the blood concentration: where the dot denotes the time derivative. The IFN blood concentration ( ) is in units of g IFN, assuming a constant blood volume and ignoring patient-specific variations in volume. For constant ( ) = 0 , administered with onset at time = 0 , this differential equation may be solved explicitly:

Pharmacodynamics of IFN
The pharmaco-dynamics (PD) of IFN is modeled such that parameter perturbation depends directly on the blood concentration of IFN. For the perturbation of and , we definêand̂as functions of the blood concentration :̂( where the parameter describes how significant the parameter is perturbed. For < 0, the parameter is reduced, while it is increased for > 0. For all values of ,̂≥ 0 is fulfilled. A different perturbation for is considered, since ≤ 1 must be maintained for all doses. We definêaŝ where determines the degree of response to treatment. For the patients considered, the time-dependent blood concentration, ( ), is described in Section 2.4. Hence, botĥ,̂, and̂are time dependent.

Procedure for obtaining individualized patient fits
Combining data for IFN dose with the parameter perturbations described by the PK and PD modeling described in Sections 2.4 and 2.5, the response of an individual patient to IFN treatment could be represented by two parameters, and . Only these parameters were considered in the fitting procedure described below, with all other parameters being fixed at the default values given in Table 1.
Three measures were considered for determining how well the model agreed with data; The disease-level error, , defined as the sum of squared errors (SSE) between model disease level and the JAK2 allele burden, the thrombocyte error, ℎ , defined as the SSE of the thrombocyte counts and the (scaled) blood cells counts and finally the leukocyte error, , defined as the SSE of leukocyte counts and the (scaled) blood cell counts.
For the 63 patients considered, personalized model fits were determined. A model simulation without treatment was shifted in time such that at = 0 the model disease agreed with the baseline measurement of the JAK2 allele burden of the patient. This could also be used to estimate the time of initial mutation and disease onset (see [28]).
Subsequently, we used an iterative three-step datafitting procedure. In all steps, the MATLAB function fminsearch was used, with options TolFun:10 −3 and TolX: 0.1. Initially, a value of was found that minimized the disease-level error, . Second, blood cell scaling factors ℎ and and treatment parameter were fitted such that the thrombocyte error, ℎ , and leukocyte error, , were minimized. As a final step, an additional optimization of with the disease-level error, , was carried out, and scaling factors ℎ and were reoptimized to minimize thrombocyte error, ℎ , and leukocyte error, . In part D of the Supporting Information, individual patient fits are shown. In all figures, the fitted values of and are shown, along with the adjusted R squared measure for each of the three error measures considered. The adjusted R squared measure is shown to allow for comparison of patient fits.
We emphasize that the procedure favors agreement between the disease burden, as given by the JAK2 allele burden in data and model disease level given by + , rather than agreement between scaled model counts of blood cells and the corresponding data.

Patient fits
Personalized values of , , , and were determined for each patient, describing the patient-specific response to IFN of the given patients. This was done for all 63 patients from the DALIAH trial. The resulting numerical solutions of the model are shown in the Supporting Information, and two examples are shown in Figure 4. In all figures, patient-IDs are used that correspond to those used in [24]. The model reproduces both the dynamics of hematologic measures (blood cell counts) and molecular measures (the JAK2 allele burden) and was found to agree with patient data irrespective of treatment response.
Good responders are characterized by a significant decrease in JAK2 allele burden and a normalization of cell counts. As changes to primarily affects the JAK2 allele burden while and affect the total cell count, the fitted values of and identify good responders. In particular, good responders are found to have fitted values lower than those of bad responders. Hence, the patient fits provide quantitative measures for how well a patient responds to treatment.
Patient-specific fits allow us to simulate hypothetical treatment-scenarios. In particular, we can consider the effect of halting treatment. In Figure 5, two simulated scenarios are shown for a particular patient, halting treatment after 0.5 and 5 years. When treatment was ceased early, cell counts rapidly relapse to elevated levels, with thrombocytes above the reference interval within half a year. After 5 years of treatment, we find that approximately 9 years without treatment is possible before the thrombocyte count exceeds the threshold. Note that cell counts at the time points chosen were approximately equal (just below 300 × 10 3 cells per L), both in the model and in data. Hence, the stage of disease at year 0.5 and year 5 would be indistinguishable, if based solely on cell counts. Our findings suggest that monitoring the JAK2 allele burden is important to determine the time to relapse and avoid premature treatment cessation.

Cell counts at steady state
For all patients, scaling factors were found, relating the modeled sum of mature cells to the blood cell counts observed in data. Scaling the steady state value of in the healthy steady state with the scaling-factor gives a patient-specific estimate of pre-disease cell counts. For most patients, these cell counts are found to be within the F I G U R E 4 Examples of patient-specific fits of the model. The left panels display the blood cell counts as gray circles. The black curve displays the sum of mature cells, + , scaled by the appropriate scaling factors. A simulation without treatment is shown in dotted black for comparison. Approximate healthy intervals of cell counts are shown in dashed gray, defined as between 145 × 10 3 and 390 × 10 3 thrombocytes per L and between 4 × 10 3 and 11 × 10 3 leukocytes per L. On the right-hand panel, the model disease level is shown together with the JAK2 allele burden data. A treatment-free simulation is shown in dotted black. The bottom-right panel depicts the pharmacokinetically modeled estimate of the IFN concentration. For all possible doses between 0 and 20 g daily IFN, model stability for the given patient-specific parameters was determined numerically. Doses for which the healthy steady state was locally stable are shown as a green background, while doses where the full-blown myeloproliferative steady state was found to be locally stable are shown in red. Panel (a) patient "P082" is depicted, with fit parameters = −0.0091, = −0.0821, = 4.1 × 10 −9 , and = 1.5 × 10 −7 . Panel (b) depicts patient "P198" which had fit-parameters = −0.004, = −0.0455, = 6.5 × 10 −9 , and = 3.1 × 10 −7 healthy interval, validating the model in the absence of disease. Scaling the steady-state value of in the fullblown myeloproliferative steady state provides an estimate for the cell count that would be approached if no treatment had been initiated. Significantly increased blood cell counts were predicted at the full-blown myeloproliferative steady state. As elevated blood cell counts are included in the diagnostic criteria, this further validates the model in the absence of treatment. Additional details along with histograms of the steady-state values are given in the Supporting Information.

Population modeling of patient response
While the model was found to agree well with most patients, the model reproduced data of a subcohort of 20 patients particularly well. These were not necessarily patients that respond well (or poorly) to treatment, but rather the subcohort most accurately reproduced by the model, in terms of the adjusted 2 of model error. Details are given in Section C of the Supporting Information. Based on this subcohort, a two-dimensional F I G U R E 5 Simulations of halted treatment can estimate the time to relapse. Based on the model fit to data for patient "P198," as shown in Figure 4A, we simulated two additional scenarios: One where treatment was halted after 0.5 year, shown in dash-dotted red, and a scenario where treatment was halted at the end of the study, year 5, shown in dash-dotted green. (For a full figure legend, see Figure 4). The colored background of the bottom-right panel shown in Figure 4 was here omitted for visual clarity F I G U R E 6 Virtual patient responses based on baseline measurement and IFN dosing shows good agreement with real patient response in some cases. Patient data for patient "P082" are shown as black circles ○. 1000 virtual patients were simulated, and the sum of mature cells were scaled to agree with the baseline data point for either leukocytes or thrombocytes. The blue curve shows the median response curve. The shaded gray areas display the distribution, with the darkest gray showing the interval from 25% to 75% of values at the given time points, the next darkest interval shows from 10% to 90% while the final interval from 5% to 95% of virtual patient responses is shown in light gray. The bottom right panel displays the modeled IFN blood concentration used for both the real patient and the virtual patients log-normal distribution of and was determined. This distribution describes the treatment-response parameters on a population level. Details of the method and the validation are given in Section C of the Supporting Information.
From the distribution, 1000 sets of dose-dependent parameter perturbations were chosen, each describing the response of a virtual patient. For each of these virtual patients, we simulated treatment schedules identical to each of the 20 patients in the subcohort. For these simulations, initial conditions corresponding to those of the baseline measurement of the real patient were used. To avoid overfitting, the fit of the given patient was not used to determine the distribution used for the virtual patients. In Figure 6

CONCLUSION
A stem cell extension of the mechanism-based Cancitis model was presented, describing simultaneously the blood production in the human body and the behavior of HSCs.
Through pharmacokinetic and pharmacodynamic modeling of the effect of IFN, the model was further extended. Considering a cohort of MPN patients and using patientspecific data for IFN doses, the model reproduced both hematologic data for blood cell counts and the JAK2 allele burden on the level of individual patients. We find that normalization of cell counts during IFN treatment could be explained by an induced decrease in self-renewal of leukemic stem cells and decreased proliferation of all progenitors cells. This finding agrees with and expands our previous work [24].
Apart from population-level parameters shared between patients, model-fits consisted of four fitted parameters: Two scaling parameters of the blood cell counts and two parameters, and . The parameter described the number of blood cells produced per differentiated stem cell and was found to relate to a transient change in blood cell counts.
determined the proliferation of malignant stem cells, and primarily affected long-term disease burden. Furthermore, we found that the time from treatment cessation until blood cell counts return to a heightened value also depended on . Patient-specific fits were made to data for the JAK2 allele burden and the thrombocyte and leukocyte counts. Such a simultaneous agreement with both hematologic and molecular markers is novel, and our work represents initial efforts of modeling different types of data of individual patients. The agreement between model and data demonstrates that mechanism-based mathematical modeling can accurately capture patient behavior, while allowing for a biological interpretation of the treatment response. The treatment response was quantified on a personalized level by the separate effect of and , with the former parameter describing the decrease in the JAK2 allele burden and a lengthening of the time to relapse, while the latter describes a transient decrease in blood cell counts. In addition, our results highlight the importance of monitoring the JAK2 allele burden, by demonstrating that relapse occurs earlier if treatment is ceased before the JAK2 allele burden is sufficiently reduced (see Figure 5).
Quantifying the twofold effect of IFN revealed a clinical challenge. A significant hematological response, as quantified by , suggests that high IFN doses could lead to an excessive decrease in blood cell counts, which could constitute a risk for the health of the patient. However, our results simultaneously suggested that for long-term molecular response and lengthening of the time to relapse, as quantified by , the IFN dose must be increased. While combination therapy or novel dose scheduling could solve this challenge, the challenge of maximizing long-term benefit while minimizing short-term risk remains. Personalized estimates of and quantifies this challenge on a personal level, identifying patients suitable for increased doses and patients that are not.
The estimates of and were based on individual patients. By relating the patient-specific estimates across the cohort of patients considered, treatment responses could be described on a population level. From this, proofof-concept population modeling was presented. In this population modeling, the effect of IFN treatment for virtual patients could be considered and simulated. Comparing the distribution of 1000 virtual patients with data from a subcohort showed that patient data could be predicted by the model using only pretreatment data and information about future IFN dose and scheduling. While further validation of our approach is necessary, predicting patient responses before treatment initiation could be an important clinical tool. Furthermore, updating and refining predictions during subsequent clinical follow-up, could give the clinician an estimate of the patient trajectory, guiding future clinical decisions.
In conclusion, we find that modeling the quantitative dynamics of blood cell counts and the JAK2 allele burden on a patient-specific level during IFN treatment is possible. In order to make accurate predictions for treatment response based on mathematical modeling, both hematologic and molecular data are needed, highlighting the importance of routinely obtaining and analyzing both clinical and molecular data. Further model validation and evaluation of the predictive power of the model are required before the model may be used as a predictive clinical tool. However, we believe that personalized mathematical modeling as described here could potentially improve the outcome of treatment with IFN in patients with MPNs. Furthermore, we believe that future work based on the modeling described here could be relevant to other diseases of the blood and thereby improve treatment outcomes of an even wider range of patients.

A C K N O W L E D G M E N T S
The authors would like to thank colleagues, clinical research units and patients who participated in the DALIAH trial.