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):
- Autoregulation feedback loop: Self-renewal of stem cells is negatively
regulated by the relative count of active stem cells.
- 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.
- Feedback loop of blood to bone marrow. This feedback is mediated by
the growth factor thrombopoietin (TPO).
- TPO is produced in the liver and the kidneys at constant rate.
- TPO activates dormant stem and progenitor cells in a delayed way,
increases the ploidy of megakaryocytes, activates dormant
megakaryocytes and suppresses platelet formation.
- 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:
- Degree 0, nadir > 150·109 cells/L.
- Degree 1, nadir between 75·109 and
150·109 cells/L.
- Degree 2, nadir between 50·109 and
75·109 cells/L.
- Degree 3, nadir between 25·109 and
50·109 cells/L.
- Degree 4, nadir < 25·109 cells/L.
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.