Figure 1l.
We present the underlying modelling concepts in brief in the following.
Semi-mechanistic haematopoiesis model
The semi-mechanistic ODE model of Henrich at al (Figure 1A) is proposed as improvement of the Friberg model. Compared to this model, it adds a slowly replicating pluripotent stem cell compartment (Stem ) at the beginning of the compartment chain intended to represent general bone marrow proliferation and maturation . This chain includes a proliferation compartment (Prol ), three equal transit compartments mirroring maturation of blasts (Transit ) and one circulating compartment (Circ ). Proliferation in Stem andProl are controlled by proliferation rate constantskstem and kprol , respectively. For the sake of parsimony, the same parameterktr describes transition rates between the five sequential compartments starting from Prol , the rate of cell division under steady state condition at the proliferation compartments and a degradation rate of circulating cells.
A feedback loop from relative counts of circulating cells is imposed on both, Stem and Prol compartment. More precisely, the term\(\left(\frac{\text{Circ}_{0}}{\text{Circ}}\right)^{\gamma}\)serves as a feedback factor of the proliferation, whereCirc0 is the steady-state value of Circand γ is a sensitivity parameter. As a consequence, more circulating cells imply reduced proliferation. Chemotherapy toxicity is modelled as a factor reducing proliferation in Stem and Prol (see Supporting information S.2 for details).
We borrowed PK models of etoposide, cyclophosphamide and doxorubicin from other authors and attached it to that model. One pharmacodynamical (PD) effect per cytotoxic drug is assumed for the proliferation compartments (Stem and Prol ) only. Relations between PD effects of different drugs are derived from other studies. We assumed inter-individual variability (IIV) for four parameters of this semi-mechanistic model based on parsimony analyses (see Table S.4.4 and explanations provided in the Supporting information S.4).
Comprehensive mechanistic thrombocytopoiesis model
Our mechanistic model is an updated and individualized version of our previous work . Model structure as well as all biological assumptions, mathematical equations and parameter estimation can be found elsewhere . In brief, the model describes the dynamics of active and dormant stem cells, colony-forming units of megakaryocytes, megakaryocytes, of different ploidies and platelets in both spleen and circulation. The model contains the following three feedback loops (see also Figure 1B and Supporting information S.1):
  1. Autoregulation feedback loop: Self-renewal of stem cells is negatively regulated by the relative count of active stem cells.
  2. Intermediate range feedback loop: There is a negative feedback loop of higher megakaryocyte numbers increasing the transition of active stem and CM cells to the respective dormant states effectively blocking further stem cell differentiations into thrombocytopoiesis lineage in case of sufficiently high numbers of megakaryocytes.
  3. Feedback loop of blood to bone marrow. This feedback is mediated by the growth factor thrombopoietin (TPO).
  4. TPO is produced in the liver and the kidneys at constant rate.
  5. TPO activates dormant stem and progenitor cells in a delayed way, increases the ploidy of megakaryocytes, activates dormant megakaryocytes and suppresses platelet formation.
  6. TPO is actively consumed by megakaryocytes and platelets in blood. It is also cleared via the kidneys.
An important novelty of our model compared to other mechanistic models of thrombopoiesis is the incorporation of a model of bone-remodeling developed by Komarova et al. . This model describes the interaction of bone-supporting osteoblasts and bone-destructing osteoclasts. The model is linked to our model of thrombocytopoiesis by assuming that osteoblasts are required to maintain dormant cells of stem cell, CM and megakaryocyte compartments. In more detail, we assume that osteoblasts define the capacity of bone marrow, i.e. higher osteoblast count implies a higher capacity for dormant cells of all types. This is supported by studies showing an involvement of megakaryocytes and osteoblasts in supporting stem cells niches .
The effect of cytotoxic drugs on bone marrow is modelled by a depletion of proliferating stem cells, progenitor cells, megakaryocytes, osteoblasts and osteoclasts. Elimination of osteoblasts through multi-cyclic chemotherapy results in a long-term reduction of the bone marrow capacity, and with it, reduction of dormant stem cells and megakaryocytes. This assumption enables us for the first time to mechanistically describe cumulative long-term toxicity as observed in multi-cycle chemotherapy. We used the same PK models of etoposide, cyclophosphamide and doxorubicin and the same relations between PD effects of different drugs as described for the semi-mechanistic model.
To calibrated our mechanistic model, we developed and applied a novel approach of integration of several clinical data sets, the so-called principle of virtual participation . In brief, we assumed that single patients also did participate in other studies and penalize large deviations from these data. Details are described in the Supporting information S.4. As a consequence, the model is required to be consistent with several other data resources which greatly stabilizes individual parameter estimates even if the individual data are relatively sparse. As a result of parsimony assumptions, a total of 12 parameters of the mechanistic model are assumed to show IIV (see Table S.2 of the Supporting information). Without applying the principle of virtual participation, the number of identifiable parameters would drop to four. Goodness of fit would also drop considerably (not shown).
We implemented and simulated the semi-mechanistic model in Matlab, using 15s ODE (ordinary differential equations) solver . Implementation and simulations of our mechanistic model is performed in R and C++, using package Rcpp. ODE equations have been approximated by the Bogacki–Shampine method .

Comparison of the predictive potential of the models regarding next cycle toxicity

Here we describe, how the models are compared regarding their predictive potential of individual next cycle toxicity. The general idea is that the models are fitted on platelet time series of patients observed in a prescribed number of cycles. Based on the resulting individual parameter estimates, we then simulated the next cycle for that patient, i.e. we predicted the thrombotoxicity for the next cycle. Individual treatment adaptations are considered throughout these processes by simulating the therapy actually applied to that patient.
Simulation results of both models are compared with observed platelet counts using a negative log-likelihood function (see section S.4 in the Supporting information). To assess the relative advantages of the models, we also performed comparisons of model predictions based on differing numbers of cycles used for parameter fitting. Quality of agreement of observed and simulated platelet counts was mainly assessed on the basis of the smallest observed cell counts because these values are typically used to define the degree of toxicity even though smaller values are possible when there is no measurement in the true nadir phase of platelets. We are confined to the following WHO degrees of toxicity for that purpose:
To evaluate the prediction quality of a model, we determined the average of the absolute differences between individually observed and predicted degrees of thrombocytopenia (also called difference of degrees (DD) in the following)
\(\text{DD}_{r,k}=\frac{\sum_{i=1}^{L_{k}}\left|\text{SD}_{i,r,k}-\text{OD}_{i,r,k}\right|}{L_{k}}\par \begin{matrix},&r=1,\cdots,N_{i}-1\par \begin{matrix},&k=r+1,\cdots,N_{i}\\ \end{matrix}\\ \end{matrix}\), (1)
where r is the number of cycles used for calibration, k is the cycle for which the predictions are made, andLk denote the number of patients analyzed in cycle k . SDi,r,k andODi,r,k correspond respectively to the simulated (predicted) and observed degrees of thrombocytopenia of patient iin cycle k . The reciprocal of this error term served as an operationalization of the predictive power of the model.
As an alternative measure of prediction quality, we only counted larger differences between observed and predicted degree of thrombocytopenia, namely differences of more than one degree. Again, we averaged these differences over the number of patients for which a prediction is made (Large difference of degrees, LDD):
\(\text{LDD}_{r,k}=\frac{\sum_{i=1,\left|\text{SD}_{i,k}-\text{OD}_{i,k}\right|>1\ }^{L_{k}}\left|\text{SD}_{i,k}-\text{OD}_{i,k}\right|}{L_{k}}\), (2)
To compare the semi-mechanistic and the mechanistic model, we introduced specific notations for the outcome measures of interest, namely {SDDr,k, {SLDDr,k} , {MDDr,k} and {MLDDr,k}. Here prefix S corresponds to the semi-mechanistic model, while prefix M corresponds to the statistics of the mechanistic model.
Due to clinical importance, we also compared receiver operating characteristics for the prediction of severe grade 3-4 thrombocytopenia for both models.

Clinical data

Primary data set
Our model comparison is based on clinical data of the NHL-B trial of elderly patients treated for aggressive non-Hodgkin’s lymphoma . Patients are randomized to one of the four arms 6xCHOP or 6xCHOEP with either 14 or 21 days of cycle duration. Schedules with 14 day cycle duration are supported by G-CSF at cycle days 4-13. Thrombocytopenia was treated with platelet transfusions, postponement of therapy or reduction in chemotherapy dose by 25% or 50%.
We selected a total of 135 patients whose platelets counts are measured during four or more cycles with at least five measurements per cycle to obtain sufficiently detailed individual time series data. Distributions of thrombocytopenia degrees of these patients in each treatment cycle is shown in Figure 2. Average degree of thrombocytopenia increases from 0.74 at cycle 1 to 1.53 at cycle 6. The fraction of severe thrombocytopenia (degrees 3-4) increases from 3% at cycle 1 to 24.3% at cycle 6. For 16 patients, the last one or two cycles could not been applied due to excessive toxicity. Nine of these had thrombocytopenia of degrees 3-4 prior to therapy withdrawal.
Validation set
Since only 8.5% of the ~1600 NHL patients met our inclusion criteria, we estimated individual parameters of additional 143 NHL patients (validation set) fulfilling relaxed inclusion criteria (monitoring during four or more cycles with at least four measurements per cycle) in order to compare their distributions with those obtained for the primary data set. In this data set, platelet transfusions were not excluded. Thus, it is necessary to account for this additional treatment by assuming an efficacy parameter of platelet transfusions, which was estimated on the basis of full time series.