Simulation model
CERONCO is a predominantly discrete multiscale computational model of cervical cancer response to treatment (external beam radiotherapy with concomitant weekly cisplatin, followed by pulsed dose rate brachytherapy) in the patientindividualized context. It makes use of the available multiscale (e.g. imaging, histological, treatment) longitudinal data of the patient. CERONCO follows a cellular automaton approach.
Core algorithms of the model can be found in our previous publications. These algorithms have been adapted and combined as necessary in order to address the particularities of the considered cervical cancer treatment. For the general features of the model see^{9,10,11} For basic concepts with regards to cisplatin chemotherapy modeling see^{11} and with regards to External Beam Radiation Therapy (EBRT) modeling see^{9,12,13}, For related sensitivity analyses see^{10,11}, A new simulation module has been developed for Pulsed Dose Rate Brachytherapy (PDRBT) modeling. This module is presented below, along with a basic outline of the main simulation model, to facilitate the understanding of the current work.
The region of interest (Gross Tumour Volume, GTV) as derived from the imaging data is represented by a threedimensional discretization mesh. The elementary volume of the mesh is called Geometrical Cell (GC). At initialization, each GC accommodates a number of biological cells (defined based on the typical cell density of 10^{9} cells/cm^{3}, unless more specific information for a particular tumour is available); each cell can belong in any of the following five classes: Stem cell; LIMP cell (LImited Mitotic Potential or committed progenitor cell); Terminally differentiated cell; Apoptotic cell; Necrotic cell. The cell cycle phases (G1, S, G2, M) and the dormant (G0) phase constitute subclasses in which stem or LIMP cells may reside.
Figure 1 depicts CERONCO’s cytokinetic decision calculator, which dictates the transitions between cell states (1hour time step). The cytokinetic model incorporates several cellularlevel phenomena: cycling of proliferating cells, symmetric and asymmetric stem cell division, terminal differentiation of LIMP cells, transition of proliferating cells to dormancy, reentrance of dormant hypoxic cells into the active cell cycle, necrosis of inadequately nourished tumour cells, spontaneous apoptosis, and cell necrosis and apoptosis due to therapy. Cell kill by EBRT is modelled based on the LinearQuadratic (LQ) model^{9,12,13} Additive toxicity of radiotherapy and chemotherapy is considered^{14}. The model offers the possibility of assigning increased radiosensitivity/chemosensitivity to stem cells compared to LIMP cells^{15,16,17}.
A new module simulating PDRBT treatment has been developed, based on the modified LQ model with correction for incomplete repair^{17,18,19,20}. Considering a fraction of PDRBT consisting of N pulses of dose d and an interpulse interval on the order of one hour, sublethal damage may not be completely repaired and the final survival fraction is given by:
$${{rm{S}}{rm{F}}}_{{rm{N}}}({rm{d}})=exp [{textstyle text{}}(alpha {rm{N}}{rm{d}}+{beta {rm{N}}{rm{G}}}_{{rm{N}}}{{rm{d}}}^{2})]$$
(1)
where G_{N} is the LeaCatcheside factor, and the α (alpha) and β (beta) parameters are the linear and the quadratic radiosensitivity coefficients, respectively, of the irradiated cells.
The LQ model is based on the curvilinear nature of doseresponse curves for the log of cell survival^{21}. It considers two cellkilling components: a linear component (alpha component) and a quadratic component (beta component). The parameters alpha and beta determine the initial slope and the degree of downward curvature, respectively, of the survival curve. According to the most common mechanistic interpretation, the yield of lethal lesions is the sum of lethal lesions produced from single radiation tracks (which are linearly related to dose, the alpha component) and lethal lesions produced from two radiation tracks (which are quadratically related to dose, the beta component); i.e. the latter quantifies the interaction of sublethal events. The dose at which these two components of cell killing are equal is the alpha/beta ratio. Since sublethal lesions can be repaired prior to resulting in a lethal event, the beta component is modified by the LeaCatcheside time factor (G_{N}) to take into account dependence on dose protraction or fractionation^{21}. Protracting the exposure time potentially allows the first lesion to be repaired before the second is produced, and the LQ approach quantifies this effect^{21}. In general, G_{N} is determined by the rate of sublethal damage repair and the particular fractionation pattern with which the dose is delivered. It is a dimensionless quantity that can take values from 0 to 1. For acute exposures ({G}_{N}to 1), and for very long exposures ({G}_{N}to 0) (“acute” and “long” are defined relative to the half time for sublethal damage repair).
Several studies addressing the issue of G_{N} factor calculation for PDRBT appear in literature^{19,22}, The computation of the G_{N} factor is based on the temporal characteristics of the dose (number of pulses, N, pulse duration, t, interpulse interval, x,) and the irradiated cells’ repair halftime T_{1/2}. Repair is assumed to follow firstorder kinetics and is modelled by a monoexponential function with rate constant μ (see also equation (7) below). The derivation of the following equation, used within CERONCO after each successive pulse to compute the survival fraction, is presented in the Supplementary Material (Section SA):
$$S{F}_{i}(d)=exp (alpha d)exp [beta {d}^{2}(i{G}_{i}(i1){G}_{i1})],,{rm{i}}=1,,ldots ,,{rm{N}}$$
(2)
G_{i} is computed (,forall ,ige 2) by the equations:
$${G}_{i}(PDR)=frac{2}{mu t}[1frac{iYS{Y}^{2}}{imu t}]$$
(3)
$$Y=1{e}^{mu t}$$
(4)
$$S=frac{iKKi{K}^{2}{e}^{mu t}+{K}^{i+1}{e}^{mu it}}{{(1K{e}^{mu t})}^{2}}$$
(5)
$$K={e}^{mu x}$$
(6)
$$mu =frac{ln2}{{T}_{1/2}}$$
(7)
where t is the duration of each pulse, x is the time between pulses without irradiation, μ is the repair rate constant, and T_{1/2} is the half time for sublethal damage repair.
For the first pulse the following equation is used^{23}:
$${G}_{1}=frac{2,[exp (mu t)+mu t1]}{{(mu t)}^{2}}$$
(8)
Equation (8) can be derived from equations (3–7) for i = 1 and x → ∞.
Equation (2) is a modification of the LQ model that can be used after each successive pulse and takes into account the current number of living tumour cells. This number is defined by the competing processes of cell death (due to radiotherapy, apoptosis, and necrosis) and cell birth as incorporated in the cytokinetic model of Fig. 1.
Following tumour initialization (section “The tumour profile concept”), at each subsequent time step the mesh is scanned and the spatiotemporal evolution rules are applied. Each complete scan can be viewed as consisting of two sequential scans^{10}. The first one updates the state of each GC by applying the rules of the cytokinetic model of Fig. 1. The second one deals with the rules governing the movement of cells throughout the tumour region. Τhe nonuniform dose distribution of BT renders a spatial handling of the tumours imperative.
A concise description of CERONCO tumour dynamics parameters is given in Table 1. A literature review has been performed to retrieve typical parameter values and value ranges for cervical cancer tumours. Table 1 includes important literaturederived quantitative information about CERONCO parameters and other tumour features whose values result from the selection of model parameter values.
Patient data
Eight patients with squamous cervical carcinoma have been included in our study (Supplementary Material Section SB.1). The patients were treated as part of the EMBRACE clinical study^{24}. The therapeutic protocol involves EBRT with concomitant cisplatin, followed by two PDRBT fractions. Followup data were not available.
The patientspecific imaging data included T2 weighted MRIderived 3D reconstructions of the Gross Tumour Volume for up to five time points:

Pretherapy (before start of EBRT)

Midterm (during EBRT)

BT0 (before start of BT)

BT1 (start of first BT fraction)

BT2 (start of second BT fraction)
These 3Dreconstruction files supply the model with the tumour’s spatial information and correspond to the region of interest onto which the discretizing mesh is superimposed. Each GC of the mesh is labeled as tumour or non/tumour. BT1 and BT2 spatial dose distribution files (total dose per GC) are provided as well. For a short outline of the procedure used to create the above files see Supplementary Material (Section SB.2).
Patientspecific treatment data included:

EBRT schedule: total dose, number of fractions, fractionation scheme (dates for the 5 fractions per week − 1 fraction per day, no irradiation during weekends)

Cisplatin administration schedule (number of cycles − once per week)

PDRBT schedule (two fractions of 20 pulses each, interpulse interval, pulse duration, date of each fraction administration). The GC pulse dose is derived from the total dose distribution file, by dividing the total dose to the GC by the number of pulses. The interval between successive pulses is 1 hour and the pulse duration is variable (0.2–0.3 hours).
All patient information was given in the context of the EUFP7 project DrTherapat, Grant agreement no. 600852. Patient information was obtained with due observance of the rights of all patients involved and in compliance with all applicable laws and regulations, including the Declaration of Helsinki as revised by the World Medical Assembly, as well as the applicable procedures and the internal guidelines of the institution. Prior to the disclosure of the patient information the clinical institution has obtained appropriate informed consents from all the patients involved, or approval from the applicable ethical review board has been obtained, all in compliance with Applicable Patient Regulations. The patients were accrued at Aarhus University Hospital and the name of the ethical committee is “Videnskabsetisk komité, Region Midt, Denmark”.
The tumour profile concept
A new virtual tumour initialization and parameter estimation workflow has been designed.
The first step is to assign values to the following three salient characteristics of the initial tumour, hereafter called “the tumour profile”:

1.
Growth Fraction
(GF, percentage of proliferating cells over all living tumour cells)

2.
Hypoxic Fraction
(HF, percentage of hypoxic cells, residing outside of the active cell cycle, over all tumour cells)

3.
Dead Fraction
(DF, percentage of dead cells over all tumour cells)
Additionally, the user selects:

4.
A cisplatin cell kill rate (CKR) value, reflecting the chemosensitivity of cancer cells to cisplatin

5.
A tumour volume doubling time (T_{d}) value or a series of T_{d} values.
Subsequently, the algorithm automatically determines sets of model parameter values that conform to these initial tumour characteristics. There exists a onetomany correspondence between these features {GF, HF, DF, Td, CKR} and possible sets of model parameter values. In the present study, the algorithm determines 100 sets of parameter values; it excludes parameter values lying outside the acceptable value ranges. For the mathematical relationships between tumour features and model parameter values see^{11} and Supplementary Section SC.
In this way, different sets of CERONCO parameter values can be grouped to different solution families, defined by the selected values of tumour features. This simulation workflow creates in essence an interface between the clinical and the mathematical reality of tumour evolution.
The tumour profilebased search can consider various regions of the parameter space which imply different tumour behaviour. Different profiles translate into different initial tumour constitutions in terms of the various cell populations (proliferating cells, dormant cells etc.), which in turn are expected to exhibit variable response to therapy and longterm evolution.
If the user cannot acquire tangible patientspecific information about the tumour features, then they can test candidate scenarios. Otherwise, the available data should be used in order to refine the initialization procedure. GF and HF estimation in cervical cancer has been intensely researched, e.g. through Ki67 studies and polarographic electrodes, respectively^{5,25,26,27,28,29,30}, Recently, imagingbased methods have been also reported^{31,32,33,34}. DF estimation can be based on MRI data, as was the case in our study. Similarly, literature abounds with methods for estimation of tumour volume doubling time^{35,36}.
In sharp contrast, information about cisplatin chemosensitivity of cervical cancer cells is rather scarce; some efforts have addressed this issue as well^{14,15,16}. Previous model sensitivity analyses indicate that this tumour feature has a profound impact on the result of therapy^{37}. We have therefore decided to expose this crucial tumour characteristic as an adjustable feature of a simulation, in order to offer the possibility of testing several explicit scenarios.
Following initialization, tumour evolution is simulated according to the patientspecific treatment data. When the simulation is complete, the algorithm checks whether there is longitudinal volumetric agreement between the simulated and the clinical tumour, taking into account possible tumour delineation errors. A solution is a set of model parameter values for which the simulated tumour’s Volume Reduction Percentage (VRP)
$$VR{P}_{sim}=[frac{{V}_{initial}^{sim},{V}_{final}^{sim}}{{V}_{initial}^{sim}}]ast 100 % $$
(9)
differs up to a predefined threshold from the corresponding VRP_{clin}, the latter calculated based on the real tumour’s GTV data, at all timepoints for which volumetric data are available. ({V}_{initial}^{sim}) and ({V}_{final}^{sim}) are the simulated initial and final tumour volume, respectively. The initial tumour volume is the pretherapy volume. The final tumour volume can correspond to any of the subsequent timepoints for which GTV is available.
The deviation thresholds between the clinical and simulated tumour VRPs have been chosen so as to reflect the existence of the abovementioned possible tumour contouring errors. Since the exact magnitude of such errors is unknown^{38}, several different criteria for volumetric compliance were tested:
$$bullet ,{bf{V}}{bf{R}}{bf{P}},{bf{5}}:VR{P}_{simulated}{mathrm{VRP}}_{{rm{clinical}}}le ,5 % ,,mathrm{for},mathrm{Midterm},{rm{BT}}0,{rm{BT}}1,{rm{BT}}2,$$
(10)
$$bullet ,{bf{V}}{bf{R}}{bf{P}},{bf{10}}:,VR{P}_{simulated}{{rm{VRP}}}_{{rm{clinical}}}le 10 % ,,{rm{for}},{rm{Midterm}},,{rm{BT}}0,,{rm{BT}}1,,{rm{BT}}2,$$
(11)
$$begin{array}{c}bullet ,{bf{M}}{bf{i}}{bf{x}}{bf{e}}{bf{d}}:,VR{P}_{simulated},,{{rm{VRP}}}_{{rm{clinical}}}le ,5 % ,,{rm{for}},{rm{Midterm}},,{rm{BT}}0\ ,and,VR{P}_{simulated},,{{rm{VRP}}}_{{rm{clinical}}}le 10 % ,,{rm{for}},{rm{BT}}1,,{rm{BT}}2end{array}$$
(12)
$$begin{array}{c}bullet ,{bf{40}} % ,{bf{v}}{bf{o}}{bf{l}}{bf{u}}{bf{m}}{bf{e}},{bf{d}}{bf{e}}{bf{v}}{bf{i}}{bf{a}}{bf{t}}{bf{i}}{bf{o}}{bf{n}}({bf{40}} % ,{bf{d}}{bf{V}}):\ ,GT{V}_{simulated}GT{V}_{clinical}le 0.4ast GT{V}_{clinical},,{rm{for}},{rm{Midterm}},,{rm{BT}}0,,{rm{BT}}1,{rm{BT}}2end{array}$$
(13)
The last criterion has been based on clinical experience with regards to delineation errors for the case of cervical cancer^{38}. Depending on the specific value of a tumour’s volume at a particular timepoint, this criterion may be stricter or more lenient than VRP 10.
The described workflow ensures a multilevel compliance of the virtual tumour with: (a) the predefined tumour features, (b) longitudinal volumetric data, (c) any tumour characteristics for which clinical data are available (e.g. in the studied cases, the diameter of the necrotic component of the initial tumour, which dictates the tumour’s DF feature), and (d) biologically plausible value ranges of the model parameters and the derived therefrom tumour characteristics. These have been retrieved from literature for the specific tumour histological type whenever possible.
The use of a performance criterion serves for quantitatively representing the agreement of a derived solution with the real volumetric data; the errormeasure Mean Absolute Error (MAE) was chosen^{39}:
$$MA{E}_{solution}=frac{{sum }_{i}VR{P}_{clinical}VR{P}_{simulated}}{N},,i={Midterm,,BT0,,BT1,,BT2},( % )$$
(14)
N: the number of timepoints with tumour volumetric data.
Solutions with lower MAE values imply better agreement with the clinical data.
Predictive use of the model: a new methodological framework Mean value parameter sets: assigning a single representative parameter value set to each patient/tumour profile/CKR value combination
For each patient, and for each tumour profile and CKR value, the parameter estimation algorithm typically identifies a large number of solutions with tumour doubling times within the acceptable value range, which all belong to a particular solution family. We have observed that it is possible to use the members of a specific solution family in order to identify a single parameter value set that could be assigned to a specific patient, tumour profile, and CKRvalue combination. This characteristic set is created by assigning to each parameter the mean of the corresponding values in the solution family.
New simulation runs have been performed using these mean value parameter sets. Their performance in terms of volumetric agreement with the real tumour, and therefore their effectiveness in representing a particular tumour profile of a patient, can be quantified by using the Mean Absolute Error (MAE) of equation (14):
$$begin{array}{rcl}MA{E}_{patient,tumourprofile,CKRvalue} & = & frac{{sum }_{i}VR{P}_{clinical}VR{P}_{mean mbox{} value mbox{} set},}{N},,( % )\ i & = & {Midterm,,BT0,,BT1,,BT2}end{array}$$
(15)
where: N: the number of timepoints with tumour volumetric data, VRP_{meanvalueset}: the VRP for the simulated mean value parameter tumour.
The mean value derived from equation (15) for all CKR values characterizes a specific patient and tumour profile:
$$MA{E}_{patient,tumourprofile}=({sum }_{i}^{N},MA{E}_{patient,tumourprofile,CKRvalue})/M,( % )$$
(16)
where M is the number of distinct CKR values having retrieved adaptation solutions.
Figure 2 is a flowchart outlining the creation of the mean value parameter sets.
Predictive tests: patient pair methodology
Since the number of patients was limited, a formal evaluation of the predictive ability of CERONCO was out of the scope of our work. Nevertheless, a set of predictive tests has been performed. These tests are based on the consideration of patient pairs, wherein the mean value set assigned to a specific patient and tumour profile is used to predict the evolution of another patient’s tumour for the same tumour profile (and vice versa).
By running simulations for patient B, Profile I, and each CKR value separately, using the corresponding mean value parameter sets of patient A, we can compute the (MA{E}_{Bleftarrow A}) error using equation 16. Similarly, an (MA{E}_{Aleftarrow B}) error characterizes the handling of patient A using patient B parameter values. A total error can be assigned to the “A and B clinical case pair”:
$$MA{E}_{Aleftrightarrow B}=frac{MA{E}_{Aleftarrow B},+,MA{E}_{Bleftarrow A}}{2},( % )$$
(17)
When for a particular profile only one of the two clinical cases of a pair has retrieved solutions, then the MAE value of the single available run is assigned to the pair error. This simplification involves 2 out of the 8 patients, and has been adopted in order to preserve the generality of the pair methodology.
We can expect that the lower MAE_{A<−>B} is, the higher is the similarity in volumetric regression terms between A and B and the bilateral predictive ability of A and B. As exemplified in the following sections, through the pair methodology the similarity of the regression profiles of any two tumours is reflected quantitatively in their (MA{E}_{Aleftrightarrow B}) values.