Documentation under construction
A wide range of ODEs models for tumour growth (TG) and tumour growth inhibition (TGI) is available in the literature and correspond to different hypotheses on the tumor or treatment dynamics. In the absence of detailed biological knowledge, selecting the most appropriate model is a challenge. In MonolixSuite2020, we provide a modular TG/TGI model library that combines sets of frequently used basic models and possible additional features. This library permits to easily test and combine different hypotheses for the tumor growth kinetics and effect of a treatment, allowing to fit a large variety of tumor size data.
When available, analytical solutions are used.
Tumor growth models
Initial tumor size
The initial tumour size TS0 can be either be:
 a parameter to estimate,
 a regressor to read from the dataset.
It has to be noted that the time corresponding to the initial tumour size is not necessarily 0. It is the case for models implemented with an analytical solution, but for models implemented with an ODE system it corresponds to the initial integration time. This time is not fixed for most models from this library, which means that it is by default the first dose or observation time in the data set, and can vary between individuals. If one wishes to fix the initial integration time for all individuals, one can uncomment the line “t_0=0” in the model to fix the initial integration time at 0, and modify the line to choose another initial time, which must necessarily be before any dose or observation in the data set. Thus only models implemented with an analytical solution can capture observations at prior times than the time for TS0.
Kinetics of tumor growth
After a growing phase, tumors experience a saturation due to limits of nutrient supply. However, this saturation property is often never measured in patients in practice because the host dies in the majority of cases before this saturation phase begins. Also in preclinics, the experiments have to be canceled if a specific tumor size is reached due to ethical reasons. Thus, tumor growth models may be divided into two broad categories: those which are able to capture the saturation as the tumor grows (achieved by the introduction a carrying capacity or a spontaneous decay component) and those which do not.

Models without saturation

Linear

The linear tumor growth assumes a constant zeroorder growth rate.
ODE system  Analytical solution  
$$ \frac{dTS}{dt} = kgl $$ $$ TS(t_0)=TS0 $$ 
$$ TS=kgl*t+TS0 $$ 


Quadratic

The quadratric tumor growth combines linear and quadratic growth rates. Because time is involved in this model, the initial time is fixed to 0 in corresponding library models.
ODE system  Analytical solution  
$$ \frac{dTS}{dt} = kgl+2*kg2*t $$ $$ TS(t_0)=TS0 $$ 
$$ TS=kgl*t+kg2*t^2+TS0 $$ 


Exponential

The exponential growth assumes the growth rate of a tumor is proportional to tumor burden (firstorder growth).
ODE system  Analytical solution  
$$ \frac{dTS}{dt} = kge*TS $$ $$ TS(t_0)=TS0 $$ undefined 
$$ TS=TS0*e^{kge*t} $$ 


Power law

If 0 < gamma < 1, the power law model (also named generalized exponential) provides a description in terms of a geometrical feature of the proliferative tissue: the growth rate is proportional to the number of proliferative cells as a fraction of the full tumor volume. The case gamma=1 corresponds to proliferative cells uniformly distributed within the tumor and leads to exponential growth. The case gamma=2/3 represents a proliferative rim limited to the surface of the tumor, where the tumor radius grows linearly in time.
ODE system  Analytical solution  
$$ \frac{dTS}{dt} = kge*TS^\gamma $$ $$ TS(t_0)=TS0 $$ 
$$ TS=(kge*(1\gamma)*t+TS0^{1\gamma})^{\frac{1}{1\gamma}} $$ 


Exponentiallinear

Experimental data of tumor growth in untreated xenograft mice suggest that tumor growth dynamics follow two distinctively different phases of growth: an initial exponential growth, followed by a linear growth phase. Biologically, this is due to the fact that the growth of the tumor is almost unlimited at first, but that as the nutrients go scarce, its growth is then necessarily restricted.
The simple exponentiallinear model captures this with a transition between an exponential (firstorder) growth to a linear (zeroorder) growth. The transition time is computed to ensure that the solution is continuously differentiable, but this is possible only if the model is not combined with any treatment effect or additional features. In that case, the transition time corresponds to TS = kgl/kge.
In the library, the exponentiallinear model is provided with its analytical solution and can not be combined to any treatment effect or additional features.
ODE system  Analytical solution  
$$ \frac{dTS}{dt} = kge*TS, t\le\tau$$ $$ \frac{dTS}{dt} = kgl, t\ge\tau$$ $$ TS(t_0)=TS0 $$ $$\tau = \frac{1}{kge}ln(\frac{kgl}{kg*TS0})$$ 
$$ TS=TS0*e^{kge*t}, t\le\tau $$
$$ TS=kgl*(t\tau)+TS0*e^{kge*\tau}, t\ge\tau $$ $$\tau = \frac{1}{kge}ln(\frac{kgl}{kg*TS0})$$ 


Simeoni

The Simeoni model approximates the exponentiallinear model with a single differential equation. The power term is fixed to 20 to allow the system to pass from the firstorder to the zeroorder growth sharply enough, as in the original system. With this value, the growth rate is approximated by kge*TS when TS is smaller than kgl/kge, and by kgl when TS is larger than kgl/kge. The advantage of the Simeoni model over the exponentiallinear model is that it is differentiable even when combined with any type of treatment effect.
When used in combination with the description of nonproliferating subpopulations of cells (transit compartments with cell distribution, or two populations model), with TS the size of the proliferating subpopulation, the growth function is slowed down by the factor TS/TotalTS only when the tumor is in its linear growth phase, but not in the exponential growth phase.
ODE system  Reference  
$$ \frac{dTS}{dt} = \frac{kge*TS}{(1+(\frac{kge}{kgl}*TS)^\psi)^(\frac{1}{\psi}} $$ $$ TS(t_0)=TS0 $$ 
Simeoni et al. (2004).
Cancer Research, 64(3), 1094–1101. 


Koch

Similarly to the exponentiallinear and Simeoni models, the Koch tumor growth model is a nonlinear growth function that follows an initial exponential growth followed by a linear growth phase. Its particularity is a smooth transition between exponential and linear growth phase, rather than a rapid switch after a threshold reached for the time (exponentiallinear) or tumor size (Simeoni).
When used in combination with the description of nonproliferating subpopulations of cells (transit compartments with cell distribution, or two populations model), with TS the size of the proliferating subpopulation, the growth function is slowed down by the factor TS/TotalTS.
ODE system  Analytical solution  Reference  
$$ \frac{dTS}{dt} = \frac{2kge*kgl*TS}{kgl+2kge*TS} $$ $$ TS(t_0)=TS0 $$ 
$$ TS=TS0e^{2kge(t+\frac{1}{kgl}(TS0TS))} $$  Koch, G. et al. (2009).
Journal of Pharmacokinetics and Pharmacodynamics, 36(2), 179–197. 

Models with saturation

Logistic

This model assumes an exponential growth rate which decelerates linearly with respect to the tumor size. This results in sigmoidal dynamics – with an initial exponential growth phase followed by a growthsaturated phase as the tumor reaches its carrying capacity (TSmax).
ODE system  Analytical solution  
$$ \frac{dTS}{dt} = kge*TS(1\frac{TS}{TSmax}) $$ $$ TS(t_0)=TS0 $$ 
$$ TS=\frac{TSmax*TS0}{TS0+(TSmaxTS0)*e^{kge*t}} $$ 


Generalized logistic

In this model, tumor growth depends on tumor size with a generalized logistic function. It reduces to the logistic model if , or to the Gompertz model when is close to 0.
ODE system  Analytical solution  
$$ \frac{dTS}{dt} = kge*TS(1(\frac{TS}{TSmax})^{\gamma}) $$ $$ TS(t_0)=TS0 $$ 
$$ TS=\frac{TSmax*TS0}{(TS0^\gamma+(TSmax^\gammaTS0^\gamma)*e^{kge*\gamma*t})^{\frac{1}{\gamma}}} $$ 


Simeonilogistic

In [HaddishBerhane et al., 2013], a hybrid model derived from the Simeoni model is proposed, that combines exponential, linear and logistic growth.
ODE system  Reference  
$$ \frac{dTS}{dt} = \frac{kge*TS*(1\frac{TS}{TSmax})}{(1+(\frac{kge}{kgl}*TS)^\psi)^{\frac{1}{\psi}}} $$ $$ TS(t_0)=TS0 $$ 
HaddishBerhane et al., 2013 


Gompertz

The Gompertz model follows the observation of a deceleration of the growth rate over time, without any phase at which the growth rate remains constant. The volume converges asymptotically to .
Several parameterizations can be used, and two are displayed on the table below. The first parameterization is the one implemented in the library, and the one for which the impacts of the parameters are the most separate and easiest to control. However, the carrying capacity TSmax can be difficult to estimate and can lack biological meaning if no clear tumor size saturation is seen. With the second parameterization, alpha and beta have an impact on the carrying capacity and on the growth rate, while TS0 also has an impact on the carrying capacity.
ODE system  Analytical solution  
$$ \frac{dTS}{dt} = TS*\beta*ln(\frac{TSmax}{TS}) $$ $$ TS(t_0)=TS0 $$ 
$$ TS=TSmax*e^{e^{\beta*t}*ln(\frac{TS0}{TSmax})} $$  
$$ \frac{dTS}{dt} = \alpha*e^{\beta*t}*TS $$
or, equivalently, $$ \frac{dTS}{dt} = (\alpha\beta*ln(TS))TS $$ $$ TS(t_0)=TS0 $$ 
$$ TS=TSmax*e^{e^{\beta*t}*ln(\frac{TS0}{TSmax})} $$ 


ExponentialGompertz

An exponentialGompertz model can also be considered, built upon the assumption that the tumor follows at first an exponential growth, and is then akin to a Gompertz model once the nutrients start to go scarce.
ODE system  Refefence  
$$ \frac{dTS}{dt} = min(kge*TS,TS*\beta*ln(\frac{TSmax}{TS})) $$ $$ TS(t_0)=TS0 $$ 
Wheldon TE. Mathematical models in cancer research. Bristol: Hilger; 1988 


Von Bertalanffy

This model is based on balance equations of metabolic processes. The growth is limited with a loss term.
ODE system  Analytical solution  Reference  
$$ \frac{dTS}{dt} = kg*TS^{2/3}kd*TS $$ $$ TS(t_0)=TS0 $$ 
$$ TS= \left[ \frac{kg}{kd}+(TS0^{1/3}\frac{kg}{kd})*e^{1/3*kd*t} \right]^3 $$  Von Bertalanffy, 1957 


Generalized Von Bertalanffy

This is the generalized version of the Von Bertalanffy model. The von Bertalanffy model corresponds to the case gamma=⅔, where the growth is proportional to the surface of the tumor.
The loss term induces for some parameter values a saturation of the tumor, with . When the loss term is neglected, the generalized von Bertalanffy model reduces to a power law.
When the Von Bertalanffy model (or generalized Von Bertalanffy model) is used in combination with a tumor inhibition model which makes use of the Norton Simon killing hypothesis, the treatment effect is only applied to the growth term of the model (and not to its decay term).
ODE system  Analytical solution  Reference  
$$ \frac{dTS}{dt} = kg*TS^{\gamma}kd*TS $$ $$ 0 \leq \gamma \leq 1 $$ $$ TS(t_0)=TS0 $$ 
$$ TS= \left[ \frac{kg}{kd}+(TS0^{1\gamma}\frac{kg}{kd})*e^{(1\gamma)*kd*t} \right]^\frac{1}{1\gamma} $$  Von Bertalanffy, 1957 
Additional features
Additional features can be included to consider more complex tumor growth models:
 Dynamic carrying capacity due to angiogenesis
 Immune dynamics causing shrinkage or oscillations of tumor size
Treatment effect
Type of treatment
Flexible treatment definition is provided thanks to several possibilities:
 No treatment
This option corresponds to simple tumor growth models.
 PK model
Dosing information is encoded in the dataset, and the effect of the treatment is modeled via a simple PK model (one compartment, bolus administration, linear elimination) which can be easily adapted to a more complex model or to fix some parameters. The predicted concentration inhibits tumor size.
Note that in case of a high number of doses, it is preferable to approximate the exposure by a piecewise regressor or by a constant treatment, to reduce computation cost. Particularly when the dosing is daily and the disease lasts several years, modeling of the PK is not useful and very costly.
 Exposure as regressor
Values which describe the evolution of the treatment’s exposure/concentration with respect to time are already stored within the data set and are used within the TGI model under the variable EXPOSURE.
DATA FORMAT RESTRICTIONS: The values describing the treatment’s exposure with respect to time must be stored in a separate column of the data set labelled ‘REGRESSOR’.
 Treatment start at t=0
The tumour growth inhibition component is applied with a constant effect after t=0.
DATA FORMAT RESTRICTIONS: The data must be formatted as to have the treatment being added at time t=0 for all individuals. Measurements before that can make use of negative values for time. Note that in that case, the initial integration time (corresponding to the initial tumour size TS0) is the first measurement time and can vary between individuals.
 Treatment start time as regressor
The time at which the treatment is added is a variable named T in the model, which may vary between individual, and whose value is read from a regressor column in the data set. When time t<T, a simple tumour growth model is being implemented. After, a tumour growth inhibition component to the model is added.
DATA FORMAT RESTRICTIONS: The values for T must be stored within the data set in a separate column tagged as ‘REGRESSOR’.
 No treatment (0) vs treatment (1) regressor
A simple tumour growth model is being implemented if the variable Trt = 0. If Trt = 1, a tumour growth inhibition component to the model is added. Trt is read from a regressor column in the dataset.
DATA FORMAT RESTRICTIONS: A 0 – if the measurement was taken before or after treatment – or 1 – if it was taken during treatment – must be matched and stored within a separate column under the header ‘REGRESSOR’ for each of the measurements of the data set.
Killing hypothesis
There are two different most commonly used hypothesis relevant to the method of killing of any given treatment within the literature:


If the growth is exponential (), the logkill and NortonSimon hypotheses result in inhibition dynamics with equivalent shapes, but the inhibition from the NortonSimon model depends on the growth rate while the logkill inhibition does not: they are equivalent if K from the logkill model is equal to K/kge in the NortonSimon model.
On the figure on the right for example, the tumor size follows an exponentiallinear growth, a constant treatment effect is applied either in the exponential phase (at t=20) or in the linear phase (at t=60). The treatment kinetics is linear. The parameters are adapted such that , thus the treatment effect is the same in the exponential phase. However if the same treatment is applied later in the linear phase, a smaller proportion of tumor is removed with the NortonSimon hypothesis, following a linear killing while the tumor remains in the linear phase.
It enfolds that under the NortonSimon hypothesis the optimal treatment should deliver the most effective dose level of drug over as short a time as possible (high dose density), since tumors given less time to grow between treatments are more likely to be eradicated.
Dynamics of treatment effect
5 different exposuredependent treatment kinetics which we found to be most common within the literature can be combined with both of these killing hypotheses.
Below are plots showing the curves resulting from the addition of these TGI models to an exponential TG model, with a punctual drug exposure resulting from a single dose which is absorbed at time 0 and then eliminated. On the first figure, a linear range of different dose amounts has been applied in order to exhibit the linear and nonlinear relationships between exposure and treatment effect.
When the treatment exposure is not given, the treatment effect is constant.
Delay of treatment effect
Chemotherapeutic effects often appear days or weeks following drug exposure. To account for a more progressive tumour regression, a delay in either the treatment effect (i.e. signal distribution model) or cell death (i.e. cell distribution model) may be added via transit compartments.

Signal distribution model
Lobo and Balthasar applied to tumor growth the transit compartment model proposed by Sun and Jusko (1998) to characterize delayed drug effects that occur via a cascade or signal transduction.The operative rate function of tumor cells killing, K3, is related to K via a series of transit compartments (K1K3). Tau refers to the mean transit time in each transit compartment.

Cell distribution model
The cell distribution model proposed in [Simeoni et al., 2004] assumes that the anticancer treatment makes some cells nonproliferating. These cells pass through different stages (also called transit compartments), characterized by progressive degrees of damage, and, eventually, they die. Hence, the attacked tumor cells have a lifespan.
The number of stages of damage and the value of tau affect the shape of the distribution of the timetodeath of damaged tumor cells. The lifespan, or timetodeath, of the attacked tumor cells, is the mean transit time that it takes for the tumor cells, affected by the action of a drug, to go through the cascade of damaging events to cell death. For n transit compartments, the average lifespan is computed by Mtt = n*tau. This library includes an implementation of the model with three transit compartments, as in [Simeoni et al., 2004].
When the model considers two populations of cancer cells in the tumor (sensitive/resistant), the treatment induces the transfer or both types of cells into the first transit compartment, with different rates.
Below are plots showing the curves resulting from the addition of both of these delays to a model making use of exponential TG component and a logkill constant TGI component depending on constant treatment starting at time 30.
Emergence of resistance
A common feature of tumours is the emergence of treatmentresistance. This phenomena can be modelled in a variety of ways. Decrease in treatment efficacy can be simply added in a phenomenological approach. More mechanistic models try to represent the cell population in its heterogeneity, splitting it into at least two subpopulations: the proliferating and the quiescent cells.
In this library, we included two common methods used within the literature to achieve this:
 Claret exponential
This simple phenomenological model accounts for the loss of druginduced decay over time due to declining efficacy of the drug, which can be due to the emergence of resistance. It uses a single resistance parameter , and it can be combined with any treatment effect.
Equation  Refefence 
Figure for an exponential growth function and constant treatment effect starting at time t=0. 
$$ K’ = K*e^{\lambda*t} $$ 
Claret et al. (2009) 
 Resistant population
This model assumes that a fraction of the tumor is resistant to the treatment, thus being killed with a smaller rate than the sensitive part of the tumor. This model uses a parameter f as initial fraction of resistant tumour cells, and the number of parameters characterizing the treatment effect doubles to have populationspecific effects.
Equation  Refefence 
Figure for an exponential growth function and constant treatment effect starting at time t=0. 
Logkilll hypothesis:
$$ \frac{dTS_s}{dt} = growth – TS_s*K_{TS_s}$$ $$ \frac{dTS_r}{dt} = growth – TS_r*K_{TS_r}$$ NortonSimon hypothesis: $$ \frac{dTS_s}{dt} = growth * (1 K_{TS_s})$$ $$ \frac{dTS_r}{dt} = growth *(1K_{TS_r})$$ Initial conditions: $$ TS_s(t_0)=TS0*(1f) $$ $$ TS_r(t_0)=TS0*f $$ 
Hahnfeldt et al (2003) 
Below is a figure showing the curves resulting from the addition of these two resistance components to a model making use of exponential TG component and a logkill linear TGI component.
Effect of treatment on angiogenesis or immune dynamics
If more complexity was added to the TG model by means of modelling immune dynamics or a dynamic carrying capacity, it can be decided that these should also be influenced by the addition of a treatment.
 DePillis immune dynamics model
If the treatment added does not only target tumour cells but a more general cell population (such as chemotherapy), it may be relevant to include a killing term on the immune cell dynamics as well.
 Active inhibition of angiogenesis
Combination therapy is often used to maximise chances of tumour regression and minimise chances of treatmentresistant clones from emerging. As such, it may be wished to model the effect of an angiogenesis inhibitor, on top of the tumour inhibitor. Here, the effect of an angiogenesis inhibitor is modelled using a simple exposureindependent firstorder log cell kill killing term. This may be complexified manually according to individual needs.
Shortcut models
Lastly, we found that some recurring models within the literature. In order to ease access to these models, we created a shortcut tab which automatically leads one to the appropriate selections when applicable, or loads the dedicated model otherwise. In each case, the initial tumor size TS0 can be set as a parameter or as a regressor.

Claret model
This model corresponds to the following choices in the library:
The version with a constant treatment effect starting at t=0, corresponding to the choice “Treatment start at t=0”, has an analytical solution impemented in the model. It is therefire much faster to compute than other models implemented as ODEs.

Simeoni model

Twopopulation model

Stein regressiongrowth model

Wang model

Bonate model

Ribba model
References
Reviews
 Ribba, B., Holford, N. H., Magni, P., Trocóniz, I., Gueorguieva, I., Girard, P., Sarr, C., Elishmereni, M., Kloft, C., & Friberg, L. E. (2014). A Review of MixedEffects Models of Tumor Growth and Effects of Anticancer Drug Treatment Used in Population Analysis. CPT Pharmacometrics Syst. Pharmacol. https://doi.org/10.1038/psp.2014.12
 Benzekry, S., Lamont, C., Beheshti, A., Tracz, A., Ebos, J. M. L., Hlatky, L., & Hahnfeldt, P. (2014). Classical Mathematical Models for Description and Prediction of Experimental Tumor Growth. PLoS Computational Biology. https://doi.org/10.1371/journal.pcbi.1003800
 Yin, A., Moes, D. J. A. R., Hasselt, J. G. C., Swen, J. J., & Guchelaar, H. (2019). A Review of Mathematical Models for Tumor Dynamics and Treatment Resistance Evolution of Solid Tumors. CPT: Pharmacometrics & Systems Pharmacology, 720–737. https://doi.org/10.1002/psp4.12450
 AlHuniti, N., Feng, Y., Yu, J., Lu, Z., Nagase, M., Zhou, D., & Sheng, J. (2020). Tumor Growth Dynamic Modeling in Oncology Drug Development and Regulatory Approval: Past, Present, and Future Opportunities. CPT: Pharmacometrics and Systems Pharmacology, 1–9. https://doi.org/10.1002/psp4.12542