#1 #1

[1]Carlos M. Nieto

[1]Oscar M. Pimentel

[1]Fabio D. Lora-Clavijo

[1]Escuela de Física, Universidad Industrial de Santander, Carrera 27 Calle 9,
Bucaramanga, 680002, Santander, Colombia

# Novel second-order model for tumor evolution: description of cytostatic and
cytotoxic effects

<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>*

###### Abstract

Cancer is a disease that takes millions of lives every year. Then, to propose
treatments, avoid recurrence, and improve the patient’s life quality, we need
to analyze this disease from a biophysical perspective with a solid
mathematical formulation. In this paper we introduce a novel deterministic
model for the evolution of tumors under several conditions (untreated tumors
and treated tumors using chemotherapy). Our model is characterized by a
second-order differential equation, whose origin and interpretation are
presented by exploiting our understanding of fluid mechanics (via continuity
equations) and the theory of differential equations. Additionally, we show
that our model can fit various experimental data sets. Thus, we prove that our
nuanced and general model can describe accelerated growth, as well as
cytostatic and cytotoxic effects. All in all, our model opens up a new window
in the understanding of tumor evolution and represents a promising connection
between the macroscopic and microscopic descriptions of cancer.

###### keywords:

Tumor evolution, cytostatic effects, cytotoxic effects, second-order model,
Sturm-Liouville problem, continuity equation

##  1 Introduction

Cancer is a disease that has been affecting a large number of humans around
the world. According to the World Health Organization, ‘in 2022, there were an
estimated 20 million new cancer cases and 9.7 million deaths’ [20]. This
alarming number of cases has to be complemented by the number of patients
whose access to health services, quality of life during the treatment, and the
possibility of cancer recurrence. Therefore, more resources should be directed
to the understanding of this disease. Given the complexity of cancer, arriving
at full control of the disease has been a challenge [20]. Studies on the
dynamics of cancer span over a large number of directions. We believe that a
deep understanding of cancer requires the collective work of medical doctors,
patients, biologists, chemists, physicists, and mathematicians. Therefore, it
is crucial to involve all areas of science to tackle this severe disease. In
particular, a mathematical description involving physical laws, and biological
and chemical processes, is needed to predict the evolution of a tumor. This
description would allow physicians to determine the best treatment for a
patient depending on the health state of the patient, the location of the
tumor, and the time of diagnosis. It has been found that the evolution of a
tumor entails a crucial number of biological processes. Among them, we have
genetic mutations, deactivation of immune responses and growth suppressors,
angiogenesis, and resistance to cell death [11].

We naturally expect the full description of cancer to be extremely complex.
Its description requires concepts from molecular biophysics, fluid mechanics,
and statistical mechanics. Several studies have been performed in this
direction [5, 12]. However, another approach has also been taken in the last
decades [15, 2]. Given the nature of experiments regarding tumor evolution,
several studies report the evolution of the size of a tumor in time. The
information regarding the volume of a tumor can be seen as the macroscopic
result of all the processes involved inside a tumor and its microenvironment.
Therefore, tumor evolution can be understood as a population problem where
cancer cells, normal cells, the immune system, and chemotherapy are at play.
That approach has been used over the last decades and it has shown to be
suitable in several scenarios. Those approaches are called volumetric models
and they involve an ordinary differential equation of a variable $V$ (volume)
as a function of time. The simplest model corresponds to the Exponential
model, in which the volume grows without control up to large values
(mathematically, it goes to infinity). Other common models are the Logistic
model, the Gompertz model, and the Von Bertalanffy [2]. However, more complex
models have been proposed over the last years (see, for instance, the work of
[19]). It is important to stress that tumor growth models play a crucial role
in understanding the complex dynamics of cancer progression and guiding
therapeutic interventions. By simulating tumor growth, researchers can gain
valuable insights into the underlying mechanisms driving cancer development
and assess the effectiveness of potential treatments. The integration of
experimental data with these models allows for personalized predictions and
represents a promising direction in understanding cancer as a whole.

Having a close look at the models discussed in the literature, we encountered
two recurrent aspects: the models are derived from a first-order differential
equation, and these models are used to describe growth and quasi-static
evolution. It would be fruitful to have a model that can describe growth,
quasi-stationary evolution, and shrinkage. It is possible to achieve such a
goal by using a piecewise function containing the standard models and other
simple functions (such as linear functions). However, having a continuous
function would be more meaningful since it can give us information about the
physics of the system. Thus, we can get a hint about the origin of such
macroscopic evolution of a tumor. On the other hand, having a higher-order
differential equation would allow for a more nuanced evolution.

In this paper we propose a model that can describe tumor evolution under
several conditions. Our model has the advantage of being adapted to several
experimental data sets related to drug treatments. Namely, our model can
account for usual growth, as well as cytostatic effects and cytotoxic effects.
To give strong arguments validating our model, we organize the paper as
follows. In section 2 we discuss the most commonly known models used to
describe tumor evolution, in section 3 we introduce our general model and
discuss how it can be used to fit several experimental data sets containing a
variety of patterns. In section 4 we discuss the possible origin and
biophysical interpretation of our model; here, we talk about the
interpretation of the parameters present in our model and their possible
connection to medical parameters. We also argue why the second-order model
nature of our model is expected from physical arguments and we exploit the
mathematical properties of our formulation. Finally, we close with the Methods
in section 5 and the Conclusions in section 6. Additionally, we add two
important appendices to understand the power of our model. In Appendix A we
discuss an interesting submodel that follows from our general formulation. In
Appendix B we show the fits for all the $85$ data sets analyzed in our study.

To make our discussion understandable, we clarify the nomenclature used in two
of the references used to extract the experimental data points. We start with
the data sets found in [16]. In that article Figure 4 is crucial. There are
two plots in that Figure: A and B. Additionally, each plot contains four
curves: control, r (red), g (green), and v (violet). The case _Control_
corresponds to an untreated tumor, while the other data points are obtained
using drug treatments. On the other hand, in [7] some abbreviations are worth
clarifying. In that reference, three different drugs were employed:

  * •

DDP: Cisplatin

  * •

BEV: Bevacizumab

  * •

PTX: Paclitaxel

However, different particular doses and combinations were used in the study.
We report the doses used in [7] with the nomenclature employed in their plots
and our plots:

  * •

Control: No drug treatment whatsoever.

  * •

DDP: Cisplatin, 3 mg/Kg every 3 weeks.

  * •

BEV: Bevacizumab, 5 mg/Kg every e week.

  * •

PTXConv: Paclitaxel, 20 mg/Kg every three weeks.

  * •

PTX Equi: Paclitaxel, 8 mg/Kg weekly.

  * •

PTX High: Paclitaxel, 12 mg/Kg weekly.

  * •

DDP/PTX Conv: 3 mg/Kg of Cisplatin with 20 mg/Kg of PTX, both once every three
weeks.

  * •

DDP/PTX Equi: 3 mg/Kg of Cisplatin once every three weeks with 8 mg/Kg of PTX
weekly.

  * •

DDP/PTX High: 3 mg/Kg of Cisplatin once every three weeks with 12 mg/Kg of PTX
weekly.

  * •

BEV/DDP/PTX Conv: 3 mg/Kg of Cisplatin and 5mg/Kg of Bevacizumab with 20 mg/Kg
of PTX, all once every three weeks.

  * •

BEV/DDP/PTX Equi: 3 mg/Kg of Cisplatin and 5mg/Kg of Bevacizumab once every
three weeks, with 8 mg/Kg of PTX weekly.

  * •

BEV/DDP/PTX High: 3 mg/Kg of Cisplatin and 5mg/Kg of Bevacizumab once every
three weeks, with 12 mg/Kg of PTX weekly.

##  2 Classical models for tumor growth

In this section we discuss the most widely used classical models for tumor
evolution, starting from the simplest to the most complex model. The models
presented in this section are characterized by first-order differential
equations. It is important to highlight this property since our work presents
arguments to go beyond this assumption. The motivation for having first-order
models relies on their simplicity of interpretation. First-order differential
equations are certainly easy to understand in terms of increasing and
decreasing functions. However, this fact does not forbid other approaches to
the problem of tumor evolution.

###  2.1 Exponential model

The Exponential model is one of the fundamental classical models for
understanding tumor evolution. It assumes that tumor growth follows an
exponential pattern, where the growth rate is directly proportional to the
size of the tumor at a given time. This model can be used to describe
uncontrolled growth and is characterized by the simple first-order
differential equation

$\frac{dV}{dt}=\alpha V,$ (1)

whose solution for $V(t)$ takes the form

$V(t)=V_{0}e^{\alpha t},$ (2)

where $V_{0}$ is the tumor volume at $t=0$ and $\alpha$ characterizes the
growth of $V(t)$. The parameter $\alpha$ can be understood as the initial
fraction rate of change in $V$. Some authors call it the _birth rate_ [8].
However, this statement is true as long as $\alpha>0$. If the parameter
$\alpha$ is negative, we have a decreasing evolution. It is worth mentioning
that this model is also recognized for its limitations, as it may not fully
capture the complexities of tumor growth dynamics, such as interactions with
the tumor micro-environment or the development of resistance to treatment.
Despite these limitations, the exponential model serves as a valuable starting
point for studying tumor evolution.

###  2.2 Logistic model

The Logistic model is another classical model commonly used to describe tumor
growth dynamics. It assumes that tumor growth initially follows an exponential
pattern, which eventually slows down leading to a plateau. This late-time
behavior is often interpreted as a competition among proliferating tumor cells
for space or nutrients [17]. The first-order differential equation of this
model contains the previous factor $\alpha V$ but it has an additional term
with a negative sign that creates a transition towards a horizontal asymptote
[18, 2, 6, 8]. The differential equation characterizing this model is usually
written as

$\frac{dV}{dt}=\alpha V\left(1-\frac{V}{K}\right),$ (3)

where $\alpha$ and $K$ are two positive constants. We see that the rate of
change has a positive contribution (leading to growth) weighted by $\alpha$
and a negative contribution (leading to a decreasing behavior) weighted by
$1/K$. The solution of this differential equation takes the form

$V(t)=\frac{KV_{0}}{V_{0}+(K-V_{0})e^{-\alpha t}},$ (4)

where $V_{0}$ is the volume at $t=0$. It is clear that, given the location of
the exponential function, this expression goes asymptotically to a fixed value
$V^{\infty}$ when $t\rightarrow\infty$. Taking the limit of (4) we get

$V^{\infty}=\lim_{t\rightarrow\infty}V(t)=K.$ (5)

The constant $K$ corresponds to the carrying capacity of the tumor and
represents the maximum and asymptotic value that the tumor volume can reach.
It is important to note that this carrying capacity is a free parameter in the
differential equation. In other words, it is not a derived quantity, contrary
to what happens in our model discussed in the following sections. The logistic
model offers a more realistic portrayal of tumor growth dynamics, albeit
simplified compared to actual biological systems’ complexities.

###  2.3 Gompertz model

The Gompertz model is another classical mathematical framework commonly
utilized to describe tumor growth. Unlike the Exponential and Logistic models,
the Gompertz model acknowledges the phenomenon of tumor growth deceleration
over time. It was first introduced by [10] and it can be presented in
different ways. In this paper we adopt the convention given in [17].
Therefore, the differential equation has the following form

$\frac{dV}{dt}=V\left[\alpha-\beta\ln\left(\frac{V}{V_{0}}\right)\right],$ (6)

where $\alpha$ and $\beta$ are two positive constants and $V_{0}$ is the
initial tumor volume. We see that the presence of the parameter $\beta$
modifies the exponential model and it counteracts the uncontrolled growth
driven by $\alpha$. The solution of (6) has the following form

$V(t)=V_{0}e^{\frac{\alpha}{\beta}(1-e^{-\beta t})},$ (7)

where $V_{0}$ is the volume at $t=0$. This model describes a tumor that has an
initial growth that ends in a horizontal asymptote at the value

$V^{\infty}=\lim_{t\rightarrow\infty}V(t)=V_{0}e^{\frac{\alpha}{\beta}}.$ (8)

Therefore, we can say that there is also a carrying capacity whose value is
$V_{0}e^{\frac{\alpha}{\beta}}$. The Gompertz model is characterized by a
sigmoidal growth curve similar to the Logistic model but with more pronounced
deceleration as the tumor approaches its maximum size.

###  2.4 Von Bertalanffy model

We finish with a model that is also widely used in population growth studies.
This model is also characterized by a first-order differential equation. This
equation contains a positive term that leads to growth and a negative term
that leads to a horizontal asymptote. This model is also a modification of the
exponential model. However, there is a change in the sign of the linear term
and the positive term is proportional to $V^{2/3}$. It is said that this last
term is related to the area of the tumor and it leads to growth. The
conventional form of the differential equation for the Von Bertalanffy model
is [3, 2, 6, 8]

$\frac{dV}{dt}=\alpha V^{2/3}-\beta V,$ (9)

where $\alpha$ and $\beta$ are two positive constants. The solution of this
differential equation can be written as

$V(t)=\left[\frac{\alpha}{\beta}+\left(V_{0}^{1/3}-\frac{\alpha}{\beta}\right)e^{-\frac{\beta}{3}t}\right]^{3},$
(10)

where $V_{0}$ is the tumor volume at $t=0$. This model also exhibits initial
growth and a transition towards a horizontal asymptote. Therefore, there is
also a carrying capacity whose expression is given by

$V^{\infty}=\lim_{t\rightarrow\infty}V(t)=\left(\frac{\alpha}{\beta}\right)^{3}.$
(11)

It is worth noticing that this maximum value does not depend on the initial
value of the tumor volume. $V^{\infty}$ depends only on the parameters
$\alpha$ and $\beta$. The Von Bertalanffy model postulates that growth is not
only influenced by intrinsic factors but also by extrinsic environmental
factors.

Other models in the literature are extensions of the previously discussed
models [6, 8]. Such models attempt to account for other properties, e.g.,
volume shrinkage. There are different ways of modifying the classical models.
Some authors introduce new variables with the respective differential equation
such that they get a system of coupled differential equations. This system of
differential equations usually represents the interaction between the
evolution of cancer cells, the action of the immune system, and the effect of
cancer therapy. Other models include time-dependent parameters in the
differential equations [6, 8]. For instance, the carrying capacity in (3) can
be taken as a time-dependent term. Additionally, it is possible to define
piecewise functions that account for the different phases that can exist in
tumor evolution (e.g., growth phase, stationary phase, and shrinkage phase).
In the next section we introduce a new model that can account for different
phases observed in experimental data. Our model comes from a second-order
differential equation and contains parameters and time-dependent functions
that resemble those in the classical models.

##  3 Second-order model

###  3.1 General model

In this section we introduce our general model for tumor evolution, aiming to
describe a wide spectrum of experimental data sets. As discussed in the
Introduction, most models describe the exponential growth period and a
stationary phase. However, experimental studies have shown that cytostatic and
cytotoxic effects can also manifest in patients undergoing cancer treatments
(e.g., chemotherapy and radiotherapy). Cytostatic effects are associated with
a stationary phase or near-stationary phase, as observed in the traditional
Logistic and Gompertz models. However, physicians primarily seek cytotoxic
effects to shrink tumors. Therefore, we also require a model that accounts for
tumor shrinkage. Moreover, other treatments have shown tumor recurrence after
a stationary or shrinkage period. Hence, we need a new model that offers a
wider spectrum of tumor dynamics and evolution.

To describe a large number of experimental data regarding tumor growth
evolution, we propose the following function for the volume as a function of
time

$V=\frac{1}{Ae^{-a\phi(t)}+Be^{b\phi(t)}},$ (12)

where $a$ and $b$ are two positive constants, $A$ and $B$ are two integration
constants and $\phi(t)$ is a function of time. The presence of the integration
constants can be understood by introducing the differential equation leading
to the function in (12). The two exponential functions in the denominator of
the expression for $V(t)$ are associated with the increasing, stationary, and
decreasing behaviors of the tumor volume as long as the parameters $a$ and $b$
are always positive and the function $\phi(t)$ is increasing, as will be
clarified in the next sections. Equation (12) can be seen as the solution of
the second-order differential equation

$\frac{d^{2}V}{dt^{2}}-\frac{2}{V}\left(\frac{dV}{dt}\right)^{2}+\left[\dot{\phi}(a-b)-\frac{\ddot{\phi}}{\dot{\phi}}\right]\frac{dV}{dt}+(\dot{\phi})^{2}abV=0.$
(13)

The second-order nature of this model makes it different from the traditional
models discussed in the literature. The consequences of this fact will be
discussed in the following sections. On the other hand, we can treat $V$ as a
function of $\phi(t)$ (whose nature will be clarified soon). The differential
equation of our model becomes

$\frac{d^{2}V}{d\phi^{2}}-\frac{2}{V}\left(\frac{dV}{d\phi}\right)^{2}+(a-b)\frac{dV}{d\phi}+abV=0.$
(14)

Since we have a second-order differential equation, the constants $A$ and $B$
are determined using the initial conditions $V(t=t_{0})=V_{0}$ and
$\dot{V}(t=t_{0})=\dot{V}_{0}$. That is, we need to know the initial volume
and its rate of change to derive the future behavior of the volume. The
expressions for the integration constants $A$ and $B$ are

$A=\frac{e^{a\phi(t_{0})}}{(a+b)V_{0}}\left[b+\frac{\dot{V}_{0}}{V_{0}\dot{\phi}(t_{0})}\right],$
(15)
$B=\frac{e^{-b\phi(t_{0})}}{(a+b)V_{0}}\left[a-\frac{\dot{V}_{0}}{V_{0}\dot{\phi}(t_{0})}\right].$
(16)

To show the power of our model, we analyze particular cases in detail. The
choice of these cases was made to describe all the experimental data found in
[7] and [16].

###  3.2 Particular models

In this subsection we discuss the applicability of specific models derived
from (12). These models are effective in describing tumor evolution under
various medical conditions. However, before presenting all the selected
particular models in our study, we introduce a form of our model for a
specific choice of $\phi(t)$ that will be relevant throughout the paper. When
we take a power-law form for this function $\phi(t)$, we get

$V=\frac{1}{Ae^{-a(t-c)^{n}}+Be^{b(t-c)^{n}}},$ (17)

where $n$ is an integer number and $c$ is any real number. The value of $n$
will determine, for instance, how fast the volume increases or decreases, and
the possible existence of a near-stationary phase. On the other hand, $c$ will
determine the transition between phases as we will see in the next sections
and subsections. For practical reasons, we introduce a notation that will be
useful in our discussion. That is, we give a name to this choice for the
function $\phi$ by defining

$\phi_{n}=(t-c)^{n}.$ (18)

When we use this form for the function $\phi(t)$ in the expressions for the
constants $A$ and $B$ of Eqs. (15-16), we get contributions to the exponential
function in Eq. (17). For instance, when $t_{0}=0$, we obtain a modified
factor that depends on $c$ within the exponential functions in (17).
Therefore, it is suitable to define

$\tilde{\phi}_{n}=\left[(t-c)^{n}-(-c)^{n}\right].$ (19)

Now we define another quantity that will be important in this work: the
relative volume at a time $t$ concerning the initial volume $V_{0}=V(t_{0})$
as follows

$V_{R}(t)=\frac{V(t)}{V_{0}}.$ (20)

In general, the resulting behavior of $V$ depends on the choices of the
parameters, the initial conditions, and the function $\phi$. Therefore, we
split the particular models into two big different sets: models with terminal
shrinkage and models without terminal shrinkage. Since one of our goals is to
understand cytostatic and cytotoxic effects from a quantitative perspective,
this two-set separation allows us to highlight the most relevant aspects of
each particular model. The two sets are characterized by the value of the
parameter $b$. If $b=0$, we obtain models without terminal shrinkage. If
$b\neq 0$, we obtain models with terminal shrinkage.

In this section we demonstrate the power of this separation and its
application to fit experimental data. In the following sections we explore the
biophysical implications of having $b=0$ or $b\neq 0$. We refer to the set
with $b=0$ as _models with terminal shrinkage_ because they all satisfy the
following property

$\lim_{t\rightarrow\infty}V_{R}(t)\rightarrow 0.$ (21)

In other words, the relative volume goes to zero at large times. This result
implies that the final volume will be much smaller than the initial volume
$V_{0}$. Essentially, there is no recurrence at that particular location in
the body. The phenomenon of recurrence or regrowth will be observed in some
models in the other set: _models without terminal shrinkage_. The models in
this big set do not have a $V$ approaching zero for large $t$. Instead, they
either exhibit a divergent behavior

$\lim_{t\rightarrow\infty}V_{R}(t)\rightarrow\infty,$ (22)

or the volume reaches a finite maximum value $V_{R}^{\infty}$

$\lim_{t\rightarrow\infty}V_{R}(t)\rightarrow V_{R}^{\infty}.$ (23)

####  3.2.1 Models without terminal shrinkage

Before starting the detailed discussion of the innovative models derived from
Eq. (17), we stress that our general Eq. (12) contains the traditional
exponential model, the Gompertz model and a modification of the Logistic
model. The _exponential model_ is obtained with the following choice:
$\phi(t)=t$, $b=0$ and $a=\dot{V}_{R0}$. The _Gompertz model_ is obtained with
the following choice: $b=0$, $\phi(t)=-e^{-\beta t}$ and
$a=\dot{V}_{R0}/\beta$ (with $\beta$ a positive constant). On the other hand,
the model that represents a _modification of the Logistic model_ is obtained
with the choice: $b=0$ and $\phi(t)=t$. In the following lines, we discuss in
detail some properties of this _modified logistic model_ and compare it with
two traditional models.

Modified logistic model

We refer to this case as the _modified logistic model_ because it shares a
similar structure with the traditional logistic model in Eq. (4). Although the
traditional and the modified logistic models have some similarities, they are
not equivalent. The explicit expression of our modified logistic model is

$V_{R}(t)=\frac{a}{\left(\dot{V}_{R0}\right)e^{-at}+\left(a-\dot{V}_{R0}\right)},$
(24)

where we have set $t_{0}=0$. It is important to note the presence of
$(a-\dot{V}_{R0})$ in the denominator. The value and sign of this term define
the overall behavior of the function $V_{R}(t)$. In the following lines, we
want to show the properties of our modified logistic model by discussing three
distinct cases and by comparing this new model with three widely used models
in the literature (exponential, logistic, and Gompertz). We rewrite here the
equations for the relative volume in each formulation for completeness

$\displaystyle V_{R}^{E}(t)=e^{\alpha t},\ \hskip
71.13188ptV_{R}^{G}(t)=e^{\frac{\alpha}{\beta}\left(1-e^{-\beta t}\right)},$
$\displaystyle V_{R}^{L}(t)=\frac{K_{R}}{1+\left(K_{R}-1\right)e^{-\alpha
t}},\ \
V_{R}^{ML}(t)=\frac{a}{\left(\dot{V}_{R0}\right)e^{-at}+\left(a-\dot{V}_{R0}\right)},$
(25)

where $K_{R}=K/V_{0}$ is the relative carrying capacity and we have used the
following labels for the models: simple exponential (E), Gompertz (G),
logistic (L), modified logistic (ML). The simplest situation in the modified
logistic model occurs when $a=\dot{V}_{R0}$. This relation leads to the usual
exponential model, as explained before. We use the expression _simple
exponential_ because we introduce extensions of this model in the following
pages. In the middle panel of Fig. 1 we show the four models for one data set
found in [16]. In this case, our model predicts that $a\approx\dot{V}_{R0}$
and all the four models have similar behavior, which is an exponential
evolution during the time period of interest. The relative carrying capacity
predicted by the logistic and modified logistic models is of the order $50$,
while the Gompertz model predicts a relative carrying capacity of the order
$8\times 10^{8}$. Thus, the evolution seen in all models is practically an
exponential growth over the first 100 days of the evolution.

In the case $a>\dot{V}_{R0}$, we get an evolution that is similar to the
logistic and Gompertz models because there is an initial fast growth and then
the volume reaches a plateau after some time. The maximum value reached by
$V_{R}$ in our modified logistic model is

$\lim_{t\rightarrow\infty}V_{R}=V_{R}^{\infty}=\frac{a}{a-\dot{V}_{R0}}.$ (26)

This stationary region brings us the notion of carrying capacity. In the case
at hand, the carrying capacity is not a parameter in the differential
equation, it is rather a quantity depending on the initial growth rate
$\dot{V}_{R0}$ and the parameter $a$. It is important to note that this
carrying capacity can go to infinity when $a=\dot{V}_{R0}$. In such a case, we
recover the exponential model. That is, an infinite carrying capacity
corresponds to uncontrolled exponential growth. In the left panel of Fig. 1 we
show the comparison of four fits using our modified logistic model, the
logistic model, the Gompertz model, and the exponential model. We see that our
modified logistic model, the logistic model, and the Gompertz model have all
similar behavior. According to the value of $R^{2}$, our modified logistic
model and the traditional logistic model show the best fit.

On the other hand, the case $a<\dot{V}_{R0}$ is both novel and interesting.
When $a<\dot{V}_{R0}$, a singularity arises that splits the function into two
branches. Focusing on the interval going from $t_{0}=0$ to the asymptote,
there is fast growth. The increase in $V_{R}$ occurs faster than in the
exponential case. If we analyze (26), we find that the value of $V_{R}$
becomes negative when $t$ approaches infinity. This result might suggest a
lack of validity of the model. However, we have to keep in mind that we are
not interested in the branch that lives beyond the asymptote. That part of the
function exists from a mathematical point of view, but it does not have any
biophysical meaning. The expression associated to the vertical asymptote is
given by

$t_{asymp}=\frac{1}{a}\ln\left[\frac{\dot{V}_{R0}}{\dot{V}_{R0}-a}\right].$
(27)

In the right panel of Fig. 1 we observe that the logistic, the Gompertz, and
the exponential models have all a similar evolution. On the contrary, our
modified logistic model stands out as the best description both visually and
quantitatively. This $a<\dot{V}_{R0}$ case is a novel evolution that can be
understood in terms of the relation between $a$ and $\dot{V}_{R0}$. As will be
explained in more detail in the sub-section 4.1, having a parameter $a$ that
is smaller than the initial growth rate is insufficient to stop the growth and
the tumor gets larger and larger with time.

Figure 1: Comparison of the models shown in Eq. (25) using data from [16] in
the left and middle panels, and data from [7] in the right panel. The dots
represent the experimental data points and the lines represent the fits as
explained in each legend. On top of each panel we wrote in magenta the
relation between $a$ and $\dot{V}_{R0}$ predicted by our modified logistic
model. Moreover, we used the label _ModiLogi_ for our modified logistic model
and _SimpExp_ for the simple/traditional exponential model

From the form of Eq. (24) and the magenta fits in Fig. 1, we can make some
conclusions regarding the role of the parameters $a$ and $\dot{V}_{R0}$. As
expected, $\dot{V}_{R0}$ determines how fast the volume will start increasing
with time, while $a$ is related to the stationary transition. There is only a
stationary region or plateau when $a>\dot{V}_{R0}$. That is, there is a
cytostatic effect when $a$ is larger than the initial growth rate. If $a$ is
smaller or equal to $\dot{V}_{R0}$, there is no stationary phase; that is,
there is no carrying capacity and the volume increases to large values when
$t$ is also large. However, there is a difference between these two cases.
When $a$ is equal to $\dot{V}_{R0}$, there is a balance that is not enough to
stop the growth of the tumor; this relation leads to traditional exponential
growth. Instead, when $a$ is smaller than $\dot{V}_{R0}$, a vertical asymptote
arises in the expression for the relative volume. This asymptote makes $V_{R}$
go quicker to very large values.

In the left panel of Fig. 1, we see that our modified logistic model predicts
a plateau and then a maximum value for $V_{R}$ that is a sign of cytostatic
effects. Taking $t_{0}=0$, the value of $t$ associated with this transition
from the growth phase to the stationary phase is equal to

$t_{tr}=\frac{1}{a}\ln\left[\frac{\dot{V}_{R0}}{a-\dot{V}_{R0}}\right].$ (28)

If the value of this time and the maximum relative volume are very high, the
tumor will have harmful effects on the patient. Therefore, one of the goals of
cancer treatment is to make $t_{tr}$ and $V_{R}^{\infty}$ as small as
possible. A specific treatment might make $t_{tr}$ and $V_{R}^{\infty}$
smaller or they might produce a different behavior. The outcome will depend on
the treatment and the living organism. From Fig. 1 we show that our modified
logistic model is an alternative option to the traditional models. However,
our original formulation goes beyond the traditional models as we explicitly
explain in the following lines.

_Generalized modified logistic model of degree $n$_:

We now discuss the straightforward generalization of our modified logistic
model. It is also obtained by taking $b=0$ in (17). We call this function the
_generalized modified logistic model of degree $n$_ because it has the same
structure as (24) but it allows for different values of the parameter $n$. The
explicit expression for this model is

$V_{R}(t)=\frac{a}{\left(\frac{\dot{V}_{R0}}{n(-c)^{n-1}}\right)e^{-a((t-c)^{n}-(-c)^{n})}+\left(a-\frac{\dot{V}_{R0}}{n(-c)^{n-1}}\right)},$
(29)

for $c\neq 0$. As clear from our explicit form of the function, this solution
is only valid when $c\neq 0$. If $c=0$, we derive an alternative equation from
(17). Thus, when $c=0$, the correct expression to use is given by

$V_{R}(t)=\frac{1}{\tilde{A}e^{-at^{n}}+\left(1-\tilde{A}\right)},$ (30)

When $n=3$, we get the form

$V_{R}(t)=\frac{a}{\left(\frac{\dot{V}_{R0}}{3c^{2}}\right)e^{-a\left[(t-c)^{3}+(c)^{3}\right]}+\left(a-\frac{\dot{V}_{R0}}{3c^{2}}\right)},$
(31)

for $c\neq 0$, or

$V_{R}(t)=\frac{1}{Ae^{-at^{3}}+\left(1-\tilde{A}\right)},$ (32)

for $c=0$. These models belong to the set of _models without terminal
shrinkage_ because they never go to zero when $t\rightarrow\infty$. According
to the value of $\Delta=\left(a-\frac{\dot{V}_{R0}}{3c^{2}}\right)$ in (31),
we can have three types of behavior. When $\Delta$ is equal to zero we get
uncontrolled growth; this case will be analyzed in detail in the following
pages and will be called _generalized exponential model of degree $n$_. When
$\Delta<0$, we get a vertical asymptote in the function; this behavior was not
found in the analyzed experimental data. However, this case might be
interesting in future studies. Finally, when $\Delta>0$ we get a function that
resembles the logistic model in the sense that $V_{R}$ reaches a plateau when
$t$ approaches infinity. However, in this case, there is interestingly an
additional stationary region. In Fig 2 we can see this function fitting an
experimental data set. We observe that the function reaches a near-stationary
region at some point. However, $V_{R}$ resumes its growth until it reaches a
definite stationary region that leads to a maximum $V_{R}^{\infty}$ when
$t\rightarrow\infty$. From a biophysical point of view, the first stationary
region of our model corresponds to a cytostatic effect on the tumor that lasts
for a short period. At some point, the volume increases again until it reaches
a second and definite stationary region. It is interesting to see that the
existence of two stationary regions is predicted by our generalized model and
that there is experimental evidence of this effect. The general expression for
the maximum volume with a general $n$ is

$V_{R}^{\infty}=\frac{a}{\left(a-\frac{\dot{V}_{R0}}{n(-c)^{n-1}}\right)}.$
(33)

Figure 2: Result of the fitting process for the data set _DDP6_ from [7] using
Eq. (31): the dots represent the experimental data and the continuous line
represents the best fit. The dashed vertical lines show the value of $c$

The model that follows corresponds to a particular case of the previous model.
However, it is important to discuss it in detail because it is connected to
interesting interpretations in terms of treatment effects. We call it
_exponential_ in the sense that it predicts a relative volume $V_{R}$ that
goes to infinity when $t\rightarrow\infty$. The way it goes to infinity
differs from the traditional exponential model but we observe uncontrolled
growth that is associated with the inefficacy or absence of treatment at large
times.

_Generalized exponential model of degree $n$_:

Now we go beyond the traditional exponential model. To generalize the
exponential model, we start with the generalized modified logistic model of
degree $n$ of Eq. (29) and take

$a=\frac{\dot{V}_{R0}}{\dot{\phi})t_{0}}=\frac{\dot{V}_{R0}}{n(t_{0}-c)^{n-1}}=\frac{\dot{V}_{R0}}{n(-c)^{n-1}},$
(34)

where we have assumed that $t_{0}=0$. Thus, the model takes the form

$V_{R}(t)=e^{\left(\frac{\dot{V}_{R0}}{n(-c)^{n-1}}\right)\left[(t-c)^{n}-(-c)^{n}\right]},$
(35)

when $c\neq 0$. If $c=0$, we have to derive a new equation starting from our
original form (17). We can also start from (32) and take $\tilde{A}=1$. Thus,
we arrive at

$V_{R}(t)=e^{at^{n}},$ (36)

These generalized exponential models allow us to have an intermediate near-
stationary period followed by fast growth. That is, they represent an
extension of the traditional exponential model, where there is only fast
growth. Several data sets from the sources used [7, 16] can be fitted with
this general model. All these fits are found with $n=2$ or $n=3$. There are
some cases in which $c=0$, and other cases in which $c\neq 0$. In the case of
a non-zero $c$, we have the function

$V_{R}(t)=e^{\left(\frac{\dot{V}_{R0}}{3(c)^{2}}\right)\left[(t-c)^{3}+(c)^{3}\right]},$
(37)

for $n=3$, and

$V_{R}(t)=e^{\left(\frac{-\dot{V}_{R0}}{2c}\right)\left[(t-c)^{2}-(c)^{2}\right]},$
(38)

for $n=2$. Clearly, in this last case $\dot{V}_{R}(0)$ must be less than zero
such that $a$ remains positive as defined since the beginning of our work. In
Fig. 3 we show one of the data sets fitted with the generalized exponential
model with $n=3$ and $c\neq 0$. There is a near-stationary region and fast
growth afterward. The value of $c$ indicates the point where there is a change
in the concavity of the function. This fact can be understood by noting that
both the first and second derivatives of $V_{R}$ vanish at $t=c$. Different
values of $n$ would modify the _near-stationary_ region. In Subsection 4.1 we
discuss the role of $n$ in more detail. Additionally, in Appendix B we show
fitted data with different values of $n$ and $c$.

Figure 3: Result of the fitting process for the data set _DDP4_ from [7] using
Eq. (37): the dots represent the experimental data and the continuous line
represents the best fit. The vertical line marks the position of $c$

####  3.2.2 Models with terminal shrinkage

We now turn our attention to models exhibiting terminal volume shrinkage. All
these models are characterized by non-zero parameters $a$ and $b$. As we did
in the previous subsection, we start from the simplest possible model and then
we move to more complex scenarios. The complexity of such models is linked to
the rich dynamics involved in the evolution of a given tumor. According to our
knowledge of chemotherapy, volume shrinkage is associated with a cytotoxic
effect of drugs on the tumor. Thus, in this subsection we deal with models
that can describe volume shrinkage due to cytotoxic drugs and the effect on
the immune system.

_Asymmetric shrinkage model of first degree_ :

We start with the simplest model. That is, we start with the case in which
$\phi(t)=t$. If we take Eq. (17) and use this linear function, we arrive at

$V_{R}(t)=\frac{(a+b)}{(b+\dot{V}_{R0})e^{-at}+(a-\dot{V}_{R0})e^{bt}},$ (39)

where we have taken $t_{0}=0$. If $a\neq b$, we obtain an asymmetric function.
That is, the growth phase and the shrinkage phase have different behavior. On
the other hand, given the increasing and decreasing regions of the function
$V_{R}$, there is a maximum at

$t_{cri}=\frac{1}{(a+b)}\ln\left(\frac{a(b+\dot{V}_{R0})}{b(a-\dot{V}_{R0})}\right).$
(40)

The position of the maximum depends on $a$, $b$, and the initial growth rate
$\dot{V}_{R0}$. Most treatments aim at reaching this maximum in the smallest
possible time. In medical terms, a smaller $t_{crit}$ generally implies a more
effective treatment, although it is important to take side effects into
account [13, 14]. In the right panel of Fig. 4, we show a data set that is
fitted with this model. There is a growth phase and a shrinkage phase.

_Symmetric shrinkage model of first degree_ :

If we take $a=b$ in (39), we get the symmetric model

$V_{R}(t)=\frac{2a}{(a+\dot{V}_{R0})e^{-at}+(a-\dot{V}_{R0})e^{at}}.$ (41)

This model is called symmetric because it presents a symmetry around the
critical time

$t_{crit}=\frac{1}{2a}\ln(A/B)=\frac{1}{2a}\ln\left(\frac{a+\dot{V}_{0}/V_{0}}{a-\dot{V}_{0}/V_{0}}\right)=\frac{1}{2a}\ln\left(\frac{a+\dot{V}_{R0}}{a-\dot{V}_{R0}}\right).$
(42)

In the left panel of Fig 4 we show a data set that is fitted using this
symmetric model. We observe a mirror symmetry concerning the position of the
maximum. The presence of such a symmetry might not be evident from (41).
However, we can show that this model presents a time-inversion symmetry.
First, we introduce a new time variable $\tau$ such that

$t=\tau+t_{crit},$ (43)

where $t_{crit}$ is the critical time given in (42). Thus, we obtain the
following expression for the relative volume

$V_{R}(\tau)=\frac{2a}{(a-\dot{V}_{R0})^{1/2}(a+\dot{V}_{R0})^{1/2}}\cdot\frac{1}{e^{-a\tau}+e^{a\tau}}$
(44)

From this result, it is clear that

$V_{R}(-\tau)=V_{R}(\tau).$ (45)

We conclude then that the function $V_{R}$ is symmetric around the critical
point (42).

Figure 4: Result of the fitting process: the dots represent the experimental
data and the continuous line represents the best fit. Left: Data set
_BEV/DDP/PTXHigh4_ from [7] fitted using Eq. (39). Right: Data set _BEV1_ from
[7] fitted using Eq. (41). The dotted vertical lines represent the points
where each function has its maximum value

It is important to note that these first two models do not describe long-term
cytostatic effects. That is, there is a sharp transition from growth to
shrinkage. We argue that this fact is due to the linear function $\phi(t)$
that we chose. Therefore, we postulate that cytostatic effects can be present
if we take a power-law form for $\phi(t)$. Thus, both cytostatic and cytotoxic
effects can be obtained by keeping $b\neq 0$ and choosing a function
$\phi(t)=(t-c)^{n}.$ We discuss this case in the next lines.

_Generalized shrinkage model of degree $n$_:

We end our discussion of particular models with a function that contains both
cytostatic and cytotoxic effects. As discussed above, we already know that a
non-zero $b$ induces a terminal-shrinkage effect and a power-law $\phi$
induces a near-stationary phase. Therefore, we take the expression

$V_{R}(t)=\frac{(a+b)}{\left(b+\frac{\dot{V}_{R0}}{n(-c)^{n-1}}\right)e^{-a\tilde{\phi}_{n}}+\left(a-\frac{\dot{V}_{R0}}{n(-c)^{n-1}}\right)e^{b\tilde{\phi}_{n}}},$
(46)

where $\tilde{\phi}_{n}$ was defined in (19). It is important to stress that
the equation is valid as long as $c\neq 0$. If $c$ is equal to zero, we have
to derive an expression that is valid in that case only. Starting from (17),
we obtain that the relative volume will depend on four numbers: $a$, $b$,
$\tilde{A}$ and $n$

$V_{R}(t)=\frac{1}{\tilde{A}e^{-at^{n}}+\left(1-\tilde{A}\right)e^{bt^{n}}},$
(47)

where the parameter $\tilde{A}$ will depend on the lowest non-zero derivative
of $V_{R}$. That is, $\tilde{A}$ will depend on $V^{(n)}(0)$.

On the other hand, since there is a growth phase and a shrinkage phase, there
should be at least a maximum value of $V_{R}$. Working with the general
equation (17), we get the following expression for the critical time

$t_{cri}=\left[\frac{1}{(a+b)}\ln\left(\frac{a\left[b+\frac{\dot{V}_{R0}}{n(-c)^{n-1}}\right]e^{a(-c)^{n}}}{b\left[a-\frac{\dot{V}_{R0}}{n(-c)^{n-1}}\right]e^{-b(-c)^{n}}}\right)\right]^{\frac{1}{n}}+c.$
(48)

Depending on the value of $n$ we can have more than one maximum If $n$ is odd,
we get only one $t_{crt}$; if $n$ is even, we get two values for $t_{cri}$. In
our fitting process, we found that $n=3$ is an optimal value for one of the
data sets in [16]. Therefore, we use the following specific model

$V_{R}(t)=\frac{(a+b)}{\left(b+\frac{\dot{V}_{R0}}{3c^{2}}\right)e^{-a\tilde{\phi}_{3}(t)}+\left(a-\frac{\dot{V}_{R0}}{3c^{2}}\right)e^{b\tilde{\phi}_{3}(t)}},$
(49)

for $c\neq 0$, or

$V_{R}(t)=\frac{1}{\tilde{A}e^{-at^{3}}+\left(1-\tilde{A}\right)e^{bt^{3}}},$
(50)

for $c=0$ In Fig 5 we see a data set that is fitted using this model with a
zero $c$. The initial growth is less tilted compared to the previous models.
In fact, the first derivative of $V_{R}$ is zero at $t_{0}=0$. This zero
derivative is the result of the form (50). Even though $\dot{V}_{R}(0)=0$, the
volume increases for $t>0$ because the third derivative is not zero. The
behavior seen with this model contains a near-stationary period, a growth
phase, and a shrinkage phase. Therefore, it can be used to describe tumors
having these three important phases. In Subsection 4.1 we discuss the role of
the parameters in this general model and their interpretation.

Figure 5: Fit of the data set _4B_ v from [16] using Eq. (50). We show the fit
in a larger domain to appreciate the form of the function

##  4 Origin and biophysical interpretation of our model

###  4.1 Parameters of the model

In this section we want to analyze the role of each parameter of the model in
Eq. (17). In subsection 3.2 we already gained an understanding of the role
played by each of the main parameters in that model. However, in the present
section, we want to explore the effect of each parameter in more detail. Thus,
we are open to the possibility of connecting each parameter with a biophysical
process in tumor evolution. Since our model contains a free function
$\phi(t)$, we stick to the models discussed in subsection 3.2.

We start with the quantity $\dot{V}_{R0}$, although is formally an initial
condition. In our discussion of the particular models, we stated that the
traditional exponential model is obtained when $n=1$ and $a=\dot{V}_{R0}$.
Therefore, the value of the initial growth rate has to be taken as a parameter
to be determined from a fitting process. If we plot the exponential model
$V_{R}(t)=e^{\dot{V}_{R0}t}$ using different values for $\dot{V}_{R0}$, we get
Fig. 6. As expected from the nature of this quantity, a high positive value
produces fast growth at large times. Given that $\dot{V}_{R0}>0$, the initial
growth factor tells us how fast the tumor will grow. It is also important to
note that the tumor will always have an increasing behavior. This situation
represents the case of an untreated tumor that is not even _attacked_ by the
immune system. To sum up, $\dot{V}_{R0}$ determines how the tumor starts
growing.

Figure 6: Effect of $\dot{V}_{R0}$ our derived traditional exponential model
$V_{R}(t)=e^{\dot{V}_{R0}t}$, which only contains the initial growth rate as a
parameter

If we introduce the parameter $a$, as done in our Modified Logistic model, we
get Eq. (24). In this model we have the two quantities: $\dot{V}_{R0}$ and
$a$. We already know what $\dot{V}_{R0}$ does. Therefore, we analyze the
effect of the parameter $a$. In Fig. 7 we observe that $a$ flattens the curve
and takes the relative volume to a constant value, as predicted by Eq. (28).
In this figure we have also marked the _transition time_. This time represents
the moment when there is a transition from fast growth and to a stationary
phase. This time is found by looking for the point at which $\ddot{V}_{R}=0$.
Fig. 7 also shows us that a large $a$ leads to a faster transition to the
stationary phase and a smaller asymptotic relative volume. The transition time
from growth to the stationary phase, as well as the value of the asymptotic
volume, depends on the relation between $\dot{V}_{R0}$ and $a$. The
expressions for both quantities can be seen in Eq. (28) and (26). Thus, we
conclude that $a$ contains the effects of the immune system and a specific
treatment leading to a global cytostatic terminal phase. In other words, the
parameter $a$ avoids uncontrolled growth.

Before moving to the other parameters, we want to show again that our Modified
Logistic model can lead to growth characterized by the presence of a vertical
asymptote. This asymptote makes the growth fast and corresponds to the case in
which $\dot{V}_{R0}>a$. In other words, this situation happens when the value
of the parameter $a$ containing information about cytostatic effects is not
enough to counteract the initial growth driven by $\dot{V}_{R0}$. In Fig. 8 we
see three curves showing tumor evolution when $\dot{V}_{R0}>a$.

Figure 7: Behavior of the Modified Logistic model of Eq. (24) for three
different values of the parameter $a$ and $\dot{V}_{R0}=1\ [1/d]$. The dotted-
dashed vertical lines represent the transition times for each case: $32.189$
hours for the magenta curve, $11.090$ hours for the black curve, and $1.330$
hours for the blue curve Figure 8: Behavior of our Modified Logistic model of
Eq. (24) when $\dot{V}_{R0}>a$ for three different values of $a$ and
$\dot{V}_{R0}=3\ [1/d]$. The three vertical dashed lines represent the
asymptotes of each function. In this case, the model is clearly valid from
$t=0$ to the position of each asymptote

Now we talk about the effect of the parameter $b$. The model that allows us to
easily understand the effect of $b$ is the one given in Eq. (39) and called
_model with terminal asymmetric shrinkage_. As understood from this name, this
model will show a late-times decreasing evolution. The reason for this
decreasing phase can be understood from the form of Eq. (39) since there is an
exponential function in the denominator that grows with time. In Fig. 9 we
observe the behavior of the model for three values of the parameter $b$. A
large $b$ leads to an evolution that approaches zero in a shorter time. Thus,
we can say that $b$ has a cytotoxic nature. This behavior was already clear in
Subsection 3.2, where we separated our model into two groups: those with
terminal shrinkage and those without terminal shrinkage. All models with
terminal shrinkage are characterized by a non-zero $b$. We want to stress that
the actual decreasing behavior depends on all the parameters since the
exponential function $e^{bt}$ is multiplied by a factor $B$. To sum up, we can
argue undoubtedly that $b$ has a cytotoxic effect as long as it is positive.

Figure 9: Behavior of the model of Eq. (28) for three different values of the
parameter $b$, $\dot{V}_{R0}=0.5\ [1/d]$ and $a=2\ [1/d]$. The dotted lines
represent the position of the maximum for each function: $1.427$ hours for the
magenta curve, $55.828$ minutes for the black curve, and $41.426$ minutes for
the blue curve

Continuing with our analysis, we address the role of the parameter $c$.
Looking at the general model of Eq. (35) we see that it represents a shift of
the function when $n>1$. However, it can induce more changes depending on the
particular values of $a$, $b$, and $n$. For instance, in Fig. 10 we see that
$c>0$ shifts the function to the right, and that a larger $c$ leads to a
longer stationary phase before the function grows rapidly to infinity. This
result is due to the $c$ dependence of the Generalized Exponential model. In
Eq. (35) we can see that the function $(t-c)^{n}$ is multiplied by the factor
that goes like $1/(-c)^{n-1}$. Then, a larger $c$ leads to a smaller factor in
the exponent.

The generalized model with terminal shrinkage of Eq. (46) also allows us to
see the effect $c$. We get a similar behavior as before. First, there is a
horizontal shift depending on the value of $c$. Additionally, the near-
stationary phase is modified by the chosen value for $c$.

Figure 10: Behavior of the model of Eq. (35) for three different values of $c$
and $\dot{V}_{R0}=0.1\ [1/d]$, $n=3$. Each dashed vertical line represents the
position of $t=c$. We observe that a large $c$ creates a longer near-
stationary phase

We also discuss the effect of $c$ in the Generalized Modified Logistic model
of Eq. (29). An interesting aspect of this model relies on the existence of
two _plateaux_ in the evolution of $V_{R}(t)$. We see in Fig. 12 how the
properties of those plateaux change with the value of $c$. Namely, the clear
separation and existence of the two plateaux are imperceptible as $c$ grows.

Figure 11: Behavior of the model of Eq. (46) for three different values of $c$
and $\dot{V}_{R0}=1\ [1/d]$, $a=2\ [1/d^{3}]$, $b=4\ [1/d^{3}]$, $n=3$. Each
dashed vertical line represents the position of $t=c$, and each dotted line
represents the position of the maximum of each function. We clearly see that
the value of $c$ modifies the evolution of the function and the properties of
the region between growth and shrinkage Figure 12: Behavior of the model of
Eq. (29) for three different values of $c$ and $\dot{V}_{R0}=3\ [1/d]$, $a=4\
[1/d^{3}]$ and $n=3$. Each dashed vertical line represents the position of
$t=c$. We clearly see that the value of $c$ modifies the evolution of the
function and the properties of the stationary regions

Now we discuss the effect of the parameter $n$, the only dimensionless
parameter in the model of Eq. (17). In Fig. 13 we see that $n$ creates a near-
stationary phase around the value of $c$ in the model (35). For $n=1$ there is
only growth. However, as long as $n>1$, there is a near-stationary phase.
During this phase, the change in $V_{R}$ is imperceptible and hence very
small. Subsequently, the function grows to infinity or goes to zero depending
on the value of $n$. We see that a large $n$ indicates the sharpest growth at
late times when the parameter $n$ is odd. However, the evenness of $n$ is also
important. If $n$ is even, the function decreases at large times. If $n$ is
even in the model (35), we get a negative value for the parameter $a$, which
is inconsistent with our initial assumptions. We can observe something similar
in Fig. 14 for the model in Eq. (29). For this generalized modified logistic
model the terminal behavior of $V_{R}$ and its intermediate near-stationary
region are also affected by the evenness of the parameter $n$.

On the other hand, in Fig. 15 we see the impact of the parameter $n$ in the
model described by Eq. (46). Once again, we see that there is a near-
stationary region as long as $n>1$. Moreover, this phase tends to last longer
when $n$ increases. Regarding the shrinkage phase, the decay is sharper when
$n$ is large. Therefore, we can say that the parameter $n$ might both
_cytostatic_ and _cytotoxic_ effects on tumor evolution.

Figure 13: Behavior of the model of Eq. (35) for seven different values of $n$
and $\dot{V}_{R0}=0.1\ [1/d]$, $c=2\ [d]$. The dashed vertical line represents
the point $t=c=2\ [d]$. A parameter $n>0$ creates a near-stationary region.
The evenness of the parameter $n$ also affects the evolution, as can be seen
for the cases $n=2,4,6$. However, these last cases are inconsistent with our
initial assumption of a positive parameter $a$ when $\dot{V}_{R0}>0$ Figure
14: Behavior of the model of Eq. (29) for six different values of $n$ and
$\dot{V}_{R0}=1\ [1/d]$, $a=4\ [1/d^{n}]$ and $c=0.6\ [d]$. The dashed
vertical line represents the point $t=c=0.6\ [d]$. A parameter $n>1$ creates a
near-stationary region. The evenness of t$n$ also affects the evolution, as
can be seen for the cases $n=2,4,6$ Figure 15: Behavior of the model of Eq.
(46) for five different values of $n$ and $\dot{V}_{R0}=0.1\ [1/d]$, $a=4\
[1/d^{n}]$, $b=8\ [1/d^{n}]$. The dashed vertical line represents the point
$c=2\ [d]$, and the dotted lines represent the position of the maxima for each
function; the case $n=4$ displays two symmetric maxima around $t=c$. It is
clear that a parameter $n>0$ creates a near-stationary region, and its
duration increases with the value of $n$

###  4.2 Second-order nature of the model

In this subsection we discuss the second-order nature of our model in more
detail. As we showed before, most models for tumor evolution come from first-
order differential equations. This fact makes the interpretation of the
differential equation clear in terms of rates. Moreover, those models require
only the knowledge of the volume at an initial time $t_{0}$. However, we know
that dynamical models often come from second-order differential equations.
Since the evolution of tumors is in-depth related to the motion of molecules
and dynamical processes, we expect the resulting volume to be described by a
function that is derived from a second-order differential equation. Being the
volume a macroscopic descriptor, it is reasonable to think that it can be
derived from the microscopic dynamics. Therefore, a natural form for the
differential equation governing $V$ would have the form

$\frac{d^{2}V}{dt^{2}}=F\left(V,\frac{dV}{dt},t\right),$ (51)

where $F$ is any function of the variables shown in the equation. This form is
precisely the one we introduced in Eq. (13).

On the other hand, it is also reasonable to expect that the evolution of the
tumor volume depends on its initial rate of change. Any dynamical problem
requires the knowledge of the initial rate of change of its variables. Hence,
our model is compatible with this physically-motivated argument. Moreover, the
necessity of knowing $\dot{V}(t_{0})$ urges us to design more accurate
experiments that allow us to get information about this initial rate of
change. Thus, not only our model is novel but also lights the way for new
experiments in the area of cancer research. In the following subsections we
show how Eq. (13) can be written to understand its origin and give
interpretations to our findings in the previous sections.

###  4.3 Sturm-Liouville problem

For a given variable $y$ depending on $x$, the Sturm’Liouville differential
equation is given by

${\cal L}y=\frac{d}{dx}\left[p(x)w(x)\frac{dy}{dx}\right]+q(x)y=-\lambda
w(x)y,$ (52)

where $p(x)$ and $q(x)$ are functions of the variable $x$, over a region of
interest $x\in[x_{i},x_{f}]$. Here the constant $\lambda$ corresponds to the
eigenvalue and $w(x)>0$ is known as the weight function. Additionally, the
operator ${\cal L}y$ is called a self-adjoint operator [1].

In order to write our Eq. (13) as a Sturm-Liouville problem, we start from Eq.
(14)

$\frac{d^{2}V}{d\phi^{2}}-\frac{2}{V}\left(\frac{dV}{d\phi}\right)^{2}+(a-b)\frac{dV}{d\phi}+abV=0,$
(53)

which was obtained by changing the independent variable from $t$ to $\phi$.
Then, we introduce the variable $y=\frac{1}{V}$. Thus, we arrive at the
differential equation

$\frac{d}{d\phi}\left[\frac{dy}{d\phi}+(a-b)y\right]=aby.$ (54)

It is straightforward to show that this equation can be written in the
following form

$\frac{d}{d\phi}\left[e^{(a-b)\phi}\frac{dy}{d\phi}\right]=abe^{(a-b)\phi}y,$
(55)

where the constant $\lambda=-ab$ corresponds to the eigenvalue.

###  4.4 Continuity equation

Besides having the property of being written as a Sturm-Liouville differential
equation, our model (13) can also be written as a continuity equation. This
form is particularly interesting it opens up the possibility of understanding
our model in terms of fluid dynamics. If we start from (55) and go back to the
variable $t$, we get the form

$\frac{d}{dt}\left[\frac{dV}{d\phi}\frac{e^{(a-b)\phi}}{V^{2}\dot{\phi}}\right]+\frac{ab\dot{\phi}e^{(a-b)\phi}}{V}=0.$
(56)

Since the expression inside the derivative of the first term is a function of
time, we can define

$Q(t)=\frac{dV}{d\phi}\frac{e^{(a-b)\phi}}{V^{2}\dot{\phi}}.$ (57)

Moreover, we can also introduce the function

$\Phi=\frac{ab\dot{\phi}e^{(a-b)\phi}}{V}.$ (58)

Thus, it is clear that the differential equation (56) can be written as

$\frac{\partial Q}{\partial t}+\Phi=0.$ (59)

This form resembles a continuity equation where $Q$ and $\Phi$ come from
integrals of the form

$Q=\int_{V}\rho\ dV,\ \ \Phi=\int_{S}\mathbf{J}\cdot\mathbf{dS}.$ (60)

Here, the quantity $\rho$ represents the volume density of a given biophysical
concept, and $\mathbf{J}$ represents the associated current density. As a
consequence, our model can be interpreted in terms of a conservation law. In
the context of cancer, we have the dynamics of fluids coming in and out of the
tumor. Hence, the existence of a continuity equation is unsurprising. It is
also interesting to note that the _flux_ $\Phi$ is multiplied by the parameter
$b$. As we saw in the previous sections, the parameter $b$ is related to
terminal-shrinkage effects as long as $\phi(t)$ is positive at large times.
Hence, if $b=0$ there is no outward flux and the volume does not decrease at
large times. In that case, we either have uncontrolled growth or a maximum
volume $V^{\infty}$ at large times. Additionally, the case $b=0$ shows us that
$Q$ is a conserved quantity or a constant of motion since

$\frac{\partial Q}{\partial t}=0.$ (61)

Therefore, we say that this quantity is a constant $Q(t)=C$. In fact, for the
case of a vanishing $b$, Eq. (56) becomes

$\frac{d}{dt}\left[\frac{dV}{d\phi}\frac{e^{a\phi}}{V^{2}\dot{\phi}}\right]=0\
\longrightarrow\frac{dV}{d\phi}\frac{e^{a\phi}}{V^{2}\dot{\phi}}=C,$ (62)

which can be easily solved and has the solution

$V(t)=\frac{1}{D+(C/a)e^{-a\phi(t)}},$ (63)

where $D$ is an integration constant. This Eq. (63) leads all the exponential
and modified logistic models discussed in the previous sections. We conclude
then that all the models without terminal shrinkage come from a conserved
quantity $C$. This interesting result will be analyzed in future work where we
will study the spatiotemporal microscopic dynamics of tumor evolution. We are
certain that this conclusion opens up questions that must be studied
theoretically and experimentally by more research groups within the areas of
Biology, Mathematics, Physics and Chemistry.

We want to close this section by discussing the opposite case, the one in
which $a=0$. In such a case, there is also a conserved quantity

$\frac{d}{dt}\left[\frac{dV}{d\phi}\frac{e^{-b\phi}}{V^{2}\dot{\phi}}\right]=0\
\longrightarrow\frac{dV}{d\phi}\frac{e^{-b\phi}}{V^{2}\dot{\phi}}=\tilde{C}.$
(64)

With this form, we find the expression for the volume

$V(t)=\frac{1}{\tilde{D}-(\tilde{C}/b)e^{b\phi(t)}},$ (65)

where $\tilde{D}$ is an integration constant. This function has the property
of being decreasing for large times. Therefore, it corresponds to a terminal-
shrinkage model as long as the pole (vertical asymptote) in the function lies
at ‘negative times’ (before the initial time $t_{0}$). This result supports
our findings regarding the nature of the parameter $b$: it promotes the
shrinkage of $V$.

##  5 Methods

In this _Methods_ section we want to explain how we obtained the parameters
from the experimental data sets used throughout the paper. As stated in
Appendix B, we used the `curve_fit` tool in Python. Since the cases discussed
in subsec. 3.2 present all characteristic features, we paired data sets with
models having similar properties from a graphical point of view. Such
properties are, for instance, terminal shrinkage or the existence of plateaux.
In our analysis we used reasonable initial-guess values for the parameters
given the discussion in subsection 4.1 and the trend observed in each data
set. Additionally, we used a large `maxfev` number of the order of thousands
to have a large number of iterations. For each of the fitted data, we made
sure that $R^{2}\gtrsim 0.7$ and the standard deviations of the parameters
were small (a relative error smaller than around $0.5$). These requirements
led us to favor some models against others. For instance, in the
BEV/DDP/PTXHigh group, we used the simple exponential model instead of our
symmetric shrinkage model for two of the data sets. The reader is invited to
check the fits shown in this paper and test our model using other experimental
studies as well.

##  6 Discussion and Conclusions

We introduced and analyzed a new model that can be used to describe the
evolution of biological populations. In this paper we showed an application to
tumor evolution. The model that we introduced is innovative in several
aspects. First, we hypothesize that a second-order model can be used to
account for different stages in tumor evolution: growth, quasi-stationary
phases, and shrinkage. Second, the parameters in our model have a
straightforward interpretation and can be used to link medical/therapeutic
concepts to the mathematical description (see subsec. 4.1, where we study the
effect of each parameter in detail). The parameters of our model combine to
produce specific patterns; however, we can undoubtedly say that the number $c$
marks the point of transition between two regions with different concavity,
$n$ can create longer regions of quasi-stationary behavior, and $b$ determines
the existence of terminal shrinkage. Third, we explore interesting properties
of our mathematical formulation that exhibit a possible origin of our model in
terms of a continuity equation. Fourth, we show that some classical models can
be derived from our general formulation. Moreover, we explicitly show that our
model(s) can be used to fit experimental data sets with diverse patterns. Our
model does not restrict to growth or stationary regions as found in the
traditional models; instead, we can account for a variety of phases and
transitions that are present in tumor evolution. In Appendix B we included the
plots of all the 85 fits of the experimental data sets used to test our
general model. We find particularly interesting those sub-models presenting
growth and shrinkage. Since most treatments attempt to reduce the volume of a
tumor, we find it beneficial to have a model describing such a phenomenon. As
discussed in subsection 3.2, terminal shrinkage is present when the parameter
$b$ is different from zero. Thus, we can connect the effectiveness of a given
treatment with $b$ and other parameters of our model. Another important aspect
that we want to highlight is the interpretation of the function $\phi(t)$. We
showed that its form drastically determines the evolution of $V_{R}$. So, we
argue that this function $\phi(t)$ is connected to the response of the
biological system to some specific treatment and external factors, and we
might see it as a _biological clock_ (a similar idea in different contexts was
discussed in [9] and [4]).

Even though our model proves to be successful in the description of tumor
evolution, further work is needed to test our formulation with more
experimental data sets. It is necessary to show that our model can describe
experiments performed with different systems and under different conditions.
Moreover, a detailed study using experiments, data analysis, and theoretical
considerations is needed to establish a clear link among the parameters of our
model and medical concepts involved in cancer treatments. In particular, it
would be extremely impactful to check the predictivity of our formulation.
That is, given a few initial experimental data points, it would be remarkable
to predict the late-times state of a tumor. Thus, physicians would be able to
assess the effectiveness of a specific treatment. Relatedly, future
experimental studies should also include a detailed analysis of the
uncertainties present in the measurements. These uncertainties would restrict
the spectrum of models that are suitable for fitting experimental data and a
more accurate description would be possible.

Another aspect that we should consider is the deterministic nature of our
formulation. We know that there are stochastic processes in these biological
systems. Therefore, small fluctuations should also be considered in future
studies. We claim that our formulation can be used as a background model in
the study of tumors, cell colonies, and other types of biological populations.
Work in this direction might be fruitful in the understanding of cancer.
Additionally, another interesting study that can be carried out entails the
possible values of the parameter $n$. In our study, we focused on integer
values; however, there are no restrictions on this parameter as long as it is
not of the form $\frac{p}{2q}$, with $p$ any odd number and $p\in\mathbb{N}$.

Acknowledgments

We would like to thank Professor Paolo Ubezio for sharing the experimental
data sets from the references [16] and [7]. F.D.L-C is supported by the
Vicerrectoría de Investigación y Extensión - Universidad Industrial de
Santander, under Grant No. 3703.

## Declarations

Conflict of interest: The authors declare that they have no conflict of
interest.

##  Appendix A Polynomial time dependence

In this Appendix we discuss a further case for the choice of $\phi(t)$. As we
exemplified above, the form of $\phi(t)$ is a key factor in the evolution of
the tumor volume. So far we have focused on the case of a power-law dependence
$\phi(t)=(t-c)^{n}$. This choice proved to be effective in the description of
the majority of data sets found in [7] and [16]. However, there is no
restriction on the form of $\phi(t)$ and we are free to choose a more complex
form. In particular, we can take a polynomial in $t$. The choice of a
polynomial function allows us to describe small variations of $V$ during its
global evolution.

Given the number of free parameters that might be involved in the case of a
polynomial $\phi(t)$, we decided to discuss such a case in this appendix.
Instead, we focused on models with a small number of free parameters with a
more straightforward interpretation. Although the polynomial case is not the
main result of our work, we found it important to present this situation in
this article.

_Polynomial exponential model:_

In our analysis of the experimental results in [7] and [16] we found that the
sets fitted with the generalized exponential model were also fitted with an
exponential $\phi(t)$. That is, we might have a _polynomial exponential
model_. However, the fit looked quite fine-tuned and some of the standard
deviations of the fitted parameters were extremely large. Therefore, we took
the simplest models in our fitting process for the majority of experimental
data sets. However, three data sets show particular evolution. These three
data sets were fitted with the case that we called _polynomial exponential
model_. This model can be derived from the main equation (12) and it can be
written as

$V_{R}(t)=e^{\phi_{pol_{n}}(t)},$ (66)

where

$\phi_{pol_{n}}(t)=\alpha
t+\alpha_{2}t^{2}+\alpha_{3}t^{3}+...+\alpha_{n}t^{n},$ (67)

and the parameter $a$ has been absorbed in $\phi_{pol_{n}}(t)$

_Interpretation of the polynomial exponential model:_

Although the possibility of having a polynomial function $\phi(t)$ to describe
the evolution of a tumor looks fine-tuned and perhaps unrealistic, we want to
propose an interpretation that opens up new directions in the study of tumors.
To avoid the presence of unexplained coefficients with very small values in
the polynomial, we rewrite Eq. (67) as

$\phi(t)=\beta_{1}\left(\frac{t}{t_{1}}\right)+\beta_{2}\left(\frac{t}{t_{2}}\right)^{2}+\cdots+\beta_{i}\left(\frac{t}{t_{i}}\right)^{i}+\cdots+\beta_{n}\left(\frac{t}{t_{n}}\right)^{n},$
(68)

where $\beta_{i}$ can be either $+1$ or $-1$, and $t_{i}$ are characteristic
times. The presence of the parameters $\beta_{i}$ removes the fine-tuning
problem of adjusting very small numbers. Then, the relevance of any term is
determined by the characteristic time $t_{i}$. If we have that
$t_{n}>t_{n-1}>\cdots>t_{2}>t_{1}$, we can associate the relevance of each
term in the polynomial to a specific time when a given effect starts taking
place. Thus, we can link the evolution of $V$ to some medical events that are
known or expected by physicians. In Fig. 16 we show an example of a data set
that was fitted with our polynomial exponential model. In this case, the
degree of the polynomial is equal to $4$. In Appendix B we present other two
cases where we used the polynomial exponential model.

Figure 16: Result of the fitting process: the dots represent the experimental
data _PTXEqui3_ from [7] and the continuous line represents the best fit using
(66). The vertical lines mark the characteristic times $t_{i}$ defined in (68)

##  Appendix B Fitting all the data sets

In this Appendix we want to show all the fits of the experimental data sets
found in references [7] and [16]. Here we use the same acronyms found in those
references, whose meaning is explained at the end of the Introduction. Our
purpose with this appendix is twofold: on the one hand, we want to show the
power of our proposal; on the other hand, we want to make clear the different
phases that can be observed in tumor evolution. This last aspect is extremely
relevant in the understanding of tumors and their response to several
treatments. All fits shown throughout the paper were found using the
`curve_fit` option in _Python_. To make the plots look more organized, we
defined some labels that stand for each model described in this paper. We
defined this notation as follows:

  * •

ModiLogi: Modified Logistic model of Eq. (24).

  * •

SimpExp: Traditional exponential model.

  * •

GenExpn: Generalized exponential model of degree $n$ of Eqs. (35) or (36)
depending on the value of $c$.

  * •

GenModiLogin: Generalized modified logistic model of degree $n$ of Eqs. (29)
or (30) depending on the value of $c$.

  * •

SymShr: Symmetric shrinkage model of Eq. (41).

  * •

AsymShr: Asymmetric shrinkage model of Eq. (39).

  * •

GenShr: Generalized shrinkage model of degree $n$ of Eqs. (46) or (47)
depending on the value of $c$.

  * •

GenExpPoln: Generalized exponential polynomial model of degre $n$ of Eq. (66).

It is important to stress that the word _degree_ stands for the power used in
$\phi(t)=(t-c)^{n}$. The only case in which a full polynomial of degree $n$ is
used corresponds to the GenExpPol model.

The following figures show the versatility of our formulation and its
usefulness in describing different phases in tumor evolution. There are two
findings that we want to highlight. First, all the _control_ sets are fitted
with our modified logistic model of Eq. (24). Second, data sets within the
same _family_ do not necessarily have identical evolution. As expected, the
response to a specific treatment depends on the living system under study; in
medical terms, the effectiveness of drug treatment will depend on the patient.
Thus, our formulation opens up the possibility of connecting some parameters
of our model to the human body response in terms of biological processes.

Figure 17: Fits of the data sets named 4A in [16]. The points represent the
experimental values and the continuous lines represent the best fits. The
dashed vertical lines in the last panel (4Av) marks the position of $c$ Figure
18: Fits of the data sets named 4B in [16]. The points represent the
experimental values and the continuous lines represent the best fits. The
dashed vertical line in 4Bg marks the position of $c$ Figure 19: Fits of the
data sets named Control in [7]. The points represent the experimental values
and the continuous lines represent the best fits. All these data points are
fitted using our modified logistic model of Eq. (24) Figure 20: Fits of the
data sets named DDP in [7]. The points represent the experimental values and
the continuous lines represent the best fits. The dashed vertical lines in
some plots mark the position of $c$. We observe that none of the models
predict a vanishing volume in the far future Figure 21: Fits of the data sets
named BEV in [7]. The points represent the experimental values and the
continuous lines represent the best fits. The dashed vertical line in BEV4
marks the position of $c$. We observe that only in BEV1 and BEV3 there is
volume shrinkage when $t\rightarrow\infty$ Figure 22: Fits of the data sets
named BEV/DDP in [7]. The points represent the experimental values and the
continuous lines represent the best fits. The dashed vertical lines in some
plots mark the position of $c$. We observe that none of the models predict a
vanishing volume in the far future Figure 23: Fits of the data sets named
PTXConv in [7]. The points represent the experimental values and the
continuous lines represent the best fits. The dashed vertical lines in some
plots mark the position of $c$. We observe that none of the models predict a
vanishing volume in the far future Figure 24: Fits of the data sets named
PTXEqui in [7]. The points represent the experimental values and the
continuous lines represent the best fits. The dashed vertical lines in some
plots mark the position of $c$. We observe that none of the models predict a
vanishing volume in the far future Figure 25: Fits of the data sets named
PTXHigh in [7]. The points represent the experimental values and the
continuous lines represent the best fits. The dashed vertical lines in some
plots mark the position of $c$. We observe that only in PTXHigh4 there is
volume shrinkage when $t\rightarrow\infty$ Figure 26: Fits of the data sets
named DDP/PTXConv in [7]. The points represent the experimental values and the
continuous lines represent the best fits. The dashed vertical lines in some
plots mark the position of $c$. We observe that none of the models predict a
vanishing volume in the far future Figure 27: Fits of the data sets named
DDP/PTXEqui in [7]. The points represent the experimental values and the
continuous lines represent the best fits. The dashed vertical lines in some
plots mark the position of $c$. We observe that none of the models predict a
vanishing volume in the far future Figure 28: Fits of the data sets named
DDP/PTXHigh in [7]. The points represent the experimental values and the
continuous lines represent the best fits. The dashed vertical lines in some
plots mark the position of $c$. We observe that there are 3 data sets for
which volume shrinkage is predicted when $t\rightarrow\infty$ Figure 29: Fits
of the data sets named BEV/DDP/PTXConv in [7]. The points represent the
experimental values and the continuous lines represent the best fits. The
dashed vertical line in the set BEV/DDP/PTXConv1 marks the position of $c$. We
observe that the majority of our fits predict volume shrinkage when
$t\rightarrow\infty$. On the other hand, the polynomial used for
BEV/DDP/PTXConv4 does not contain the liner term Figure 30: Fits of the data
sets named BEV/DDP/PTXEqui in [7]. The points represent the experimental
values and the continuous lines represent the best fits. We observe that all
our fits predict volume shrinkage when $t\rightarrow\infty$ Figure 31: Fits of
the data sets named BEV/DDP/PTXHigh in [7]. The points represent the
experimental values and the continuous lines represent the best fits. We
observe that all our fits predict volume shrinkage when $t\rightarrow\infty$

## References

  *   * Arfken et al [2012] Arfken GB, Weber HJ, Harris FE (2012) Mathematical Methods for Physicists (Seventh Edition). Academic Press, Boston, 10.1016/C2009-0-30629-7
  * Benzekry et al [2014] Benzekry S, Lamont C, Beheshti A, et al (2014) Classical mathematical models for description and prediction of experimental tumor growth. PLoS computational biology 10(8):1–19. 10.1371/journal.pcbi.1003800
  * von Bertalanffy [1957] von Bertalanffy L (1957) Quantitative laws in metabolism and growth. The Quarterly Review of Biology 32(3):217–31. 10.1086/401873
  * Castorina et al [2006] Castorina P, Delsanto PP, Guiot C (2006) Classification scheme for phenomenological universalities in growth problems in physics and other sciences. Phys Rev Lett 96(18):188701. 10.1103/PhysRevLett.96.188701
  * Cotner et al [2023] Cotner M, Meng S, Jost T, et al (2023) Integration of quantitative methods and mathematical approaches for the modeling of cancer cell proliferation dynamics. American Journal of Cell Physiology 324(2):C247–C262. 10.1152/ajpcell.00185.2022
  * Evain and Benzekry [2016] Evain S, Benzekry S (2016) Mathematical modeling of tumor and metastatic growth when treated with sunitinib. [Research Report] Inria Bordeaux Sud-Ouest https://hal.science/hal-01377994
  * Falcetta et al [2017] Falcetta F, et al (2017) Modeling cytostatic and cytotoxic responses to new treatment regimens for ovarian cancer. Cancer research 77(23):6759–6769. 10.1158/0008-5472.CAN-17-1099
  * Ghaffari et al [2022] Ghaffari L, et al (2022) Classical mathematical models for prediction of response to chemotherapy and immunotherapy. PLOS Computional Biology 18(2):e1009822. 10.1371/journal.pcbi.1009822
  * Gillooly et al [2002] Gillooly JF, Charnov EL, West GB, et al (2002) Effects of size and temperature on developmental time. Nature 417(6884):70–73. 10.1038/417070a
  * Gompertz [1825] Gompertz B (1825) On the nature of the function expressive of the law of human mortality, and on a new mode of determining the value of life contingencies. Philosophical Transactions of the Royal Society of London 115:513–583. 10.1098/rspl.1815.0271
  * Hanahan and Weinberg [2022] Hanahan D, Weinberg RA (2022) Hallmarks of cancer: new dimensions. Cancer Discovery 12(1):31–46. 10.1158/2159-8290.CD-21-1059
  * Jørgensen et al [2023] Jørgensen A, et al (2023) Data-driven spatio-temporal modelling of glioblastoma. Royal Society Open Science 10:221444. 10.1098/rsos.221444
  * Markman [1996] Markman M (1996) What does tumor shrinkage mean to the patient receiving chemotherapy? Cleveland Clinic Journal of Medicine 63(5):301. https://www.ccjm.org/content/ccjom/63/5/301.full.pdf
  * Ooki et al [2022] Ooki A, et al (2022) Impact of early tumor shrinkage on quality of life in patients treated with frst-line cetuximab plus chemotherapy for unresectable metastatic colorectal cancer: results of phase ii quack trial. BMB Cancer 22(1):711. 10.1186/s12885-022-09811-x
  * Sachs et al [2001] Sachs R, Hlatky L, Hahnfeldt P (2001) Simple ode models of tumor growth and anti-angiogenic or radiation treatment. Mathematical and Computer Modelling 33:1297–1305. 10.1016/S0895-7177(00)00316-2
  * Ubezio [2019] Ubezio P (2019) Beyond the t/c ratio: Old and new anticancer activity scores in vivo. Cancer research 11:8529–8538. 10.2147/CMAR.S215729
  * Vaghi et al [2020] Vaghi C, Rodallec A, Fanciullino R, et al (2020) Population modeling of tumor growth curves and the reduced gompertz model improve prediction of the age of experimental tumors. PLoS computational biology 16(2):e1007178. 10.1371/journal.pcbi.1007178
  * Verhulst [1838] Verhulst PF (1838) Notice sur la loi que la population poursuit dans son accroissement. Correspondance Mathématique et Physique 10:113–121
  * West and Newton [2019] West J, Newton PK (2019) Cellular interactions constrain tumor growth. Proceedings of the National Academy of Sciences of the United States of America 116(6):1918–1923. 10.1073/pnas.1804150116
  * WorldHealthOrganization [2024] WorldHealthOrganization (2024) Global cancer burden growing, amidst mounting need for services. https://www.who.int/news/item/01-02-2024-global-cancer-burden-growing--amidst-mounting-need-for-services [Accessed: 31st March 2024]