Mathematical Modeling of MPNs Offers Understanding and Decision Support for Personalized Treatment

(1) Background: myeloproliferative neoplasms (MPNs) are slowly developing hematological cancers characterized by few driver mutations, with JAK2V617F being the most prevalent. (2) Methods: using mechanism-based mathematical modeling (MM) of hematopoietic stem cells, mutated hematopoietic stem cells, differentiated blood cells, and immune response along with longitudinal data from the randomized Danish DALIAH trial, we investigate the effect of the treatment of MPNs with interferon-α2 on disease progression. (3) Results: At the population level, the JAK2V617F allele burden is halved every 25 months. At the individual level, MM describes and predicts the JAK2V617F kinetics and leukocyte- and thrombocyte counts over time. The model estimates the patient-specific treatment duration, relapse time, and threshold dose for achieving a good response to treatment. (4) Conclusions: MM in concert with clinical data is an important supplement to understand and predict the disease progression and impact of interventions at the individual level.

A reduced model is proposed and analysed by Ottesen et al. (2019). The model reduction is based on quasisteady-state approximations ofẋ 1 ,ẏ 1 ,ȧ andṡ. The quasi-steady-state approximations imply the following algebraic expression: where κ = d x0 x 0 + d x1 x 1 + d y0 +d y0 y 0 + d y1 y 1 . The two-dimensional system is then simply: Cell-counts are interpreted from the model as: Dose-dependent perturbation of model parameters Some differences in dosing and timing of IFN treatment occurred in the DALIAH trial. Dosing information was approximated as the average daily dose, but to have gradual changes in IFN-perturbed parameters, we made use of a simple pharmakokinetic model of IFN. Absorption and clearance of IFN was modelled to occur with the same rate, τ , independently of the current dose. The IFN blood-concentration can hence be described by: where D(t) is the dose at time t and B(t) is the blood concentration at the given time. Note that D(t) is a stepsize constant function describing the average daily dose. Considering an interval t ∈ [t 0 , t 1 ] where D(t) is constant, D(t) = D 0 , the differential equation (7) has solution: for t ∈ [t 0 , t 1 ]. This allowed for a very simple estimate of the blood concentration level with smooth transitions from one level of treatment to another. The parameterd y0 was modelled to increase with the current blood-concentration of IFN. This was done by defining a functiond y0 (ρ, t) such that:d where B(t) is the time-dependent blood-level concentration of IFN and ρ is a parameter describing the strength of the patient-specific response to IFN treatment. Hence, a large ρ corresponds to a significant increase tod y0 during treatment, while a low ρ signifies a low response. To ensure positive parameters, only positive values of ρ were considered.

Fitting procedure
The JAK2 V617F allele burden was interpreted in the model as i.e. as the contribution of leukemic mature cells to the sum of mature cells. Hence, solutions of equations (5) together with calculation of x 1 and y 1 from equations (4) can be used to determine a model estimate of the JAK2 V617F allele burden.
The model was simulated numerically in MATLAB using the function ode15s. One healthy stem cell, x 0 was removed and a leukemic stem cell y 0 was added at time 0. Due to the default parameters used, the system approaches a leukemic steady-state asymptotically. The JAK2 V617F allele burden is initially close to 0 and approaches 1 as time increases.
The development of MPN before diagnosis was assumed to be equal for all patients and hence that the untreated disease development to be the same for the population. The initial baseline measurement of the JAK2 V617F allele burden at diagnosis was used to determine where on the common disease development curve the patient were before treatment initiation. Hence data and model were shifted in time such that baseline JAK2 V617F coincides with the model. Simulations of treatment-response used the model variables at the corresponding time as initial conditions.
The model was simulated for each patient over the time-period where data for the given patient was available. By separately modelling the blood-concentration of IFN as described above, and simulating a dose-dependent perturbation ofd y0 as shown in equation (9), the model could be simulated for different values of ρ, yielding different responses and developments of the JAK2 V617F for the simulated time-period.
Model error was determined by the difference between the measurements of the JAK2 V617F allele burden and the model estimate at the same time-points. Note that the baseline our methodology forces the model to be in agreement with baseline measurement. The value of ρ found to (locally) minimize the sum of squared errors was found using the MATLAB function fminsearch.

Cell-count measurements
The sum of mature cells in the model was assumed to correlate with both the cell-counts of thrombocytes and of leukocytes, see equations (6). As fluctuations between cell-lines were expected, and indeed observed in data, a complete correlation with both cell-counts cannot be expected for all measurements.
By assuming a linear scaling, we were able to determine two scaling-parameters, K leukocytes and K thrombocytes , that minimized the sum of squares difference between cell-count data and the sum of x 1 (t) and y 1 (tx) the JAK2 V617F-error minimizing fits described above. These scaling parameters then describe the proportion of mature cells in the model which can be attributed to the given cell-line; thrombocytes or leukocytes.
As some significant differences were observed between thrombocytes and leukocytes (e.g. one patient had increasing thrombocyte-counts but decreasing leukocytes counts over the study period), such linear scaling of cell counts provides a rough estimate of the general level of mature blood-cells, and cannot be expected to be in perfect agreement with data. However, for a subset of patients, these ad-hoc blood-data fits were found to agree well with both JAK2 V617F data and blood-cell data simultaneously even though the model was only optimised for agreement with JAK2 V617F data, see supplementary material C.
In summary, a single parameter,d y0 was modelled to increase linearly for increasing IFN-dose. An errorminimizing fitting procedure of the JAK2 V617F allele burden identified a patient-specific response-parameter describing the strength of the perturbation tod y0 . A subsequent fit of the resulting model estimate of mature blood-cells counts to the thrombocyte and leukocyte counts found excellent agreement between the response to IFN in data and in the model for most patients. Hence a mechanistic interpretation as IFN resulting in an dose-dependent increase of malignant stem-cell clearing (d y0 ) appears as (part of) a plausible explanation for the treatment-responses observed in patient data.
Supplementary material B: Fits to JAK2 V617F data Below we show fits of the Cancitis model to all available data for the JAK2 V617F allele burden. In the fits, the PD-effect of the IFN was modelled to increase the LSC death rate, dependent on the current blood concentration of IFN, modelled by a simple PK model as described in the main text.
The top panel shows patient-data and the model disease burden. The dotted black curves shows disease development, assumed to be the same for all patients. Baseline JAK2 V617F is shifted in time to coincide with this curve. Grey * signify data-points used for the fitting procedure, while black * signify later data-points not used for the procedure, shown only for comparison and evaluation of the model fit. The full black line shows the best fit to the given data, defined as the fit that minimizes the sum of squared distance between the used data-points and the model at the given time-points, see the supplementary describing the mathematical methods. For time after the last data-points used for fitting, the actual future dosing information was used.
The bottom panel shows the PK-modelling of the blood-concentration of IFN as a black curve. The background color signifies the long-term state that the system will move toward: For doses where the background is green, the model implies that long-term the healthy state will be approached, while the red background shows doses that will approach full-blown disease. The unit of the dose is the average daily dose of IFN in µg. No distinction was made between Peginterferon-α-2a (Pegasys) and Peginterferonα-2b (PegIntron) At the top of all figures, a goodness-of-fit-metric of the best fit is shown, defined as the sum of squared errors, divided by the number of points used in the fit. A low number signifies a good fit. Only the error of points used (the grey points) were included in the error estimate.

Patient P199
Supplementary material C: Fits to leukocyte and thrombocyte data Here we show the best fit to JAK2 V617F allele burden, along with the subsequent fits to leukocytes and thrombocytes. The right-hand panels displays the full model fit to all data, see the description in supplementary material B. The two left-hand panels display the sum of mature cells in the corresponding model simulation, scaled to fit with either the leukocytes or the thrombocyte. The resulting scaling-factor hence describes the assumed fraction of mature cells that correspond to the given cell-type. In the figures, the approximate healthy intervals are also shown, defined as between 4 · 10 3 to 11 · 10 3 per microliter of blood for leukocytes and 145 · 10 3 to 390 · 10 3 per microliter of blood thrombocytes. We emphasize that the patient-specific response to IFN in the model was only fit to the data for the JAK2 V617F allele burden, and the scaling of the mature cell counts was done subsequently, independent of the IFN blood-concentration.
At the top of the leukocytes and thrombocyte figures, a goodness-of-fit-metric of the best fit is shown, defined as the sum of squared errors, divided by the number of points used in the fit. A low number signifies a good fit. Note that the error metric depends on the order of magnitude of the data, and hence the metric is typically higher for thrombocytes than for leukocytes.