# Cancer recurrence times from a branching process model

Stefano Avanzini and Tibor Antal  School of Mathematics, University of
Edinburgh, Edinburgh, EH9 3JZ, UK

###### Abstract

As cancer advances, cells often spread from the primary tumor to other parts
of the body and form metastases. This is the main cause of cancer related
mortality. Here we investigate a conceptually simple model of metastasis
formation where metastatic lesions are initiated at a rate which depends on
the size of the primary tumor. The evolution of each metastasis is described
as an independent branching process. We assume that the primary tumor is
resected at a given size and study the earliest time at which any metastasis
reaches a minimal detectable size. The parameters of our model are estimated
independently for breast, colorectal, headneck, lung and prostate cancers. We
use these estimates to compare predictions from our model with values reported
in clinical literature. For some cancer types, we find a remarkably wide range
of resection sizes such that metastases are very likely to be present, but
none of them are detectable. Our model predicts that only very early
resections can prevent recurrence, and that small delays in the time of
surgery can significantly increase the recurrence probability.

###### Contents

  1. 1 Introduction
  2. 2 Results
    1. 2.1 Model setup
    2. 2.2 Time to reach detectable size
    3. 2.3 Exponential population growth
    4. 2.4 Relapse time with resection
    5. 2.5 Metastasis classification
  3. 3 Discussion
    1. 3.1 Parameter estimation
    2. 3.2 Model predictions
    3. 3.3 Conclusions
  4. 4 Materials and methods
    1. 4.1 Single type process
    2. 4.2 Scaled relapse time distribution
    3. 4.3 Explicit results for exponential primary growth
      1. 4.3.1 Scaled relapse time
      2. 4.3.2 Small initiation limit
      3. 4.3.3 Mean relapse time
      4. 4.3.4 Numerical computation

##  1 Introduction

Metastases develop as cancer cells disseminate from a primary tumor and
establish new malignant lesions in the surrounding tissue or at other sites
[1]. However, the full process of metastasis formation is much more complex
and many related aspects are not yet fully understood. In particular, it is
still unclear whether metastases are initiated during early or late stages of
carcinogenesis (see e.g. [2, 3, 4]). These details, however, affect the
chances of a patient presenting detectable or undetectable metastases at
diagnosis, which in turn influence treatment strategies and prognosis. For
these reasons, different authors (see e.g. [5, 6] and the references therein)
have proposed mathematical models to improve our understanding of the dynamics
of metastasis formation.

Metastases frequently arise in cancer patients, and their occurrence greatly
diminishes the chances of effective treatment. In fact, even when a therapy is
initially successful, metastases often lead to relapse and are responsible for
an estimated $90\%$ of cancer related deaths [7]. Despite this common disease
progression, reliable predictions for cancer recurrence rates and times are
still lacking [8].

Recently, many generalizations of the Luria-Delbrück model [9] have been
employed to study specific traits of tumor evolution, such as the development
of drug resistance [10, 11, 12, 13], the role of driver mutations [14, 15] and
metastasis formation [5, 6, 16, 17]. Another line of research focused on
temporal features, after the first stochastic model for the time to tumor
onset was proposed by Armitage and Doll in their pioneering work on
carcinogenesis [18]. A few decades later authors began to investigate
stochastic models of tumor latency time. In particular, these works led to
mathematical descriptions of optimal schedules of cancer surveillance [19,
20], cure rates [21] and cancer recurrence [22]. These models are studied in
the context of survival analysis and reviewed in the excellent book of
Yakovlev and Tsodikov [23].

In this paper we build a model for cancer recurrence by joining these two
approaches, that is we use Luria-Delbrück type models to study cancer relapse
times. In particular, we consider a deterministically growing tumor seeding
metastases at a rate depending on its size [24], and model the evolution of
each metastasis (or clone) as independent birth-death branching processes. A
similar setup was used by Lea and Coulson to mimic mutations occurring in a
growing bacterial population [25]. In our model though we interpret these
mutation events (from wild-type cells to mutants) as metastasis initiation
events. The distribution of mutant close sizes was studied with an
exponentially growing wild-type population [26] and with more general wild-
type growth function [16]. Kendall [27] also allowed the wild-type population
to grow stochastically, but this extension left the mutant behavior unchanged
for small initiation (mutation) rates [28, 29]. Hence in this paper we model
the size of the primary tumor as a deterministic function (focusing on
exponential and logistic growth as examples), while allow the seeded
metastases to evolve stochastically according to branching processes.

Within this framework we study the time to cancer relapse, defined as the
interval between the primary onset and the first time that any of the
metastases reaches a fixed detectable size. Similar characterizations are
employed in the threshold models described in [23, 22].

The rest of the paper is organized as follows: In Results we present our
mathematical model of metastases initiation and growth, and derive an explicit
formula for the probability distribution of the time to relapse. We then
extend our model to include the resection of the primary tumor at a given time
and distinguish between synchronous and metachronous metastases. In Discussion
we report parameter estimates for five different cancer types (namely breast,
colorectal, headneck, lung and prostate) and analyze the corresponding
predictions yielded by our model. Quantitative results are compared with data
collected from clinical literature. In Material and Methods we present details
about the mathematical formulation of our model and related derivations.

##  2 Results

Our mathematical characterization of the time to cancer recurrence is based on
a stochastic model of metastasis formation. We first present the fundamental
assumptions and features of this model, and then use them to derive the
probability distribution of the time to relapse.

###  2.1 Model setup

Figure 1: Sample realisation of the model obtained by simulations. The primary
tumor grows according to a deterministic exponential function $n(t)$ \-
depicted by the blue line. It initiates distant metastases at rate $\nu n(t)$,
and each of them grows as an independent branching process (only the first
five are plotted). The first time $\tau$ that any of these metastases reaches
a minimal detectable size $M$ is defined as the time to cancer relapse. Also,
the primary tumor is surgically removed at a given time $T$, when it is made
of $N=n(T)$ cells. In the realisation shown, the third established metastases
(green curve) is the first to reach detectable size, and hence determines the
time to cancer relapse $\tau$. Based on clinical data (summarized in Table 1),
we estimated model parameters (summarized in Table 2), and here we use those
for colorectal cancer, with $N=2\times 10^{11}$. Note that a similar
illustration for metastasis formation appears in [30].

We model the number of cells in the primary tumor as a deterministic function
of time $n(t)$. The tumor initiates metastases at rate $\nu n(t)$, where $\nu$
is constant. We implicitly assume that all tumor cells can metastasize at the
same rate. Since we make no assumptions on $n(t)$, one can define initiation
at rate $\nu n(t)^{\gamma}$ to model scenarios where only a fraction of the
primary tumor can metastasize, for example only the cells near its surface or
close to blood vessels (see e.g. [6]). The initiated metastases are then
modelled as independent branching birth-death processes [31], all with the
same birth rate $\alpha$ and death rate $\beta$. We assume that they are
supercritical, that is they have a positive net growth rate
$\lambda=\alpha-\beta>0$, and consequently grow exponentially for large times
[31]. Exponential growth has also been observed in clinical studies [32] for
untreated human lung metastases, which supports our modelling choice.

Under these assumptions each metastasis will eventually go extinct with
probability $q=\beta/\alpha<1$. The surviving ones instead grow unboundedly
and will reach any given size [31]. Let $M$ be a fixed number of cells
representing the minimal detectable size of a cancerous lesion. We aim to
describe the time to cancer recurrence, defined as the first time $\tau$ that
any metastasis reaches the detectable size $M$.

The minimal detectable size $M$ is typically very large, with estimates over
$10^{6}$ (see parameter estimations in Discussion). As the probability that a
large supercritical population goes extinct is negligibly small, we assume
that each metastasis survives indefinitely if it reaches $M$. Then, due to the
splitting property of Poisson processes, the surviving metastases that
eventually reach the detectable size are initiated as a non-homogeneous
Poisson process $(K_{t})_{t\geq 0}$ with rate $\nu(1-q)n(t)$. Here $K_{t}$
denotes the number of metastases established by $t$, conditioned on survival.
The expected number of established metastases at time $t$ is thus

$a_{t}=\mathrm{E}[K_{t}]=\nu(1-q)\int_{0}^{t}n(s)ds$

and the probability that at least one is present at $t$ is equal to

$\,\mathrm{P}(K_{t}\geq 1)=1-e^{-a_{t}}$ (1)

Surviving metastases are initiated at times $\sigma_{i}:=\inf\\{t\geq
0:K_{t}=i\\}$ and are described by i.i.d. birth-death processes
$(S_{i}(s))_{s\geq 0}$, where $S_{i}(s)$ is the number of cells in the $i$-th
metastasis at time $s$ after its establishment. In particular, we have
$S_{i}(0)=1$ for every $i$. For each of these processes we define
$\Theta_{i}:=\inf\\{s\geq 0:S_{i}(s)=M\\}$ as the time needed by the $i$-th
established metastasis to grow to the detectable size $M$, counting again from
its initiation. Since the processes $S_{i}(s)$ are independent, the hitting
times $\Theta_{i}$ are also independent and identically distributed. As shown
in Material and Methods, for $M$ large their distribution asymptotically
satisfies

$\,\mathrm{P}(\Theta_{i}\leq t\mid\Omega_{\infty}^{(i)})\sim G(t)\equiv
e^{-(1-q)Me^{-\lambda t}}$ (2)

where $\Omega_{\infty}^{(i)}$ denotes the eventual survival for the $i$-th
metastasis. Interestingly, the distribution $G(t)$ is of a Gumbel type, which
generally describes the maximum of independent random variables with
exponential (right) tail. This Gumbel type has two parameters, $a$ and $b>0$,
and distribution function $\exp(-e^{-\frac{x-a}{b}})$. Hence, conditioned on
survival, we asymptotically have
$\Theta_{i}\sim\mathrm{Gumb}_{\mathrm{max}}\left(\frac{\log
M(1-q)}{\lambda},\frac{1}{\lambda}\right)$ for every $i$.

###  2.2 Time to reach detectable size

Given the definitions in the previous section, we have that the $i$-th
metastasis reaches detectable size at time $\tau_{i}:=\sigma_{i}+\Theta_{i}$,
measured from primary onset. Metastases are initiated at time $s$ at rate
$\nu(1-q)n(s)$ and then reach the detectable size before $t$ with probability
$G(t-s)$ for $s\leq t$. Hence, the thinning property of Poisson processes
yields that metastases which become detectable by time $t$ are initiated at
time $s$ at rate $\nu(1-q)n(s)G(t-s)$. Consequently, the number of metastases
detectable by $t$ follows a Poisson process $(S_{s})_{0\leq s\leq t}$ with
respect to time $s$ for a fixed $t$. In particular, the number $S_{t}$ of such
metastases established by $t$ is thus a Poisson random variable with mean

$b_{t}=\mathrm{E}[S_{t}]=\nu(1-q)\int_{0}^{t}n(s)G(t-s)ds$ (3)

The relapse time is defined as the first time that any metastasis reaches the
detectable size, $\tau:=\min_{i}\\{\tau_{i}\\}$. Hence, $\tau$ is smaller than
$t$ if by that time at least one metastasis that becomes detectable before $t$
is initiated, and so

$\,\mathrm{P}(\tau\leq t)=\,\mathrm{P}(S_{t}\geq 1)=1-e^{-b_{t}}$ (4)

A sample realisation of our model, including the relapse time $\tau$, is
depicted in Figure 1. In the large detectable size $M$ limit, the relapse time
distribution converges to a simpler form (see Material and Methods)

$\tau-\frac{1}{\lambda}\log M\xlongrightarrow{d}\overline{\tau}$

where the random variable $\bar{\tau}$ is distributed as

$\,\mathrm{P}(\bar{\tau}\leq
t)=1-\exp\left({-\nu(1-q)\int_{-\infty}^{t}n(t-s)e^{-(1-q)e^{-\lambda
s}}ds}\right)$ (5)

Hence for large $M$ the relapse time decomposes as
$\tau\approx\frac{1}{\lambda}\log M+\overline{\tau}$ into a deterministic part
which depends only on $\lambda$ and $M$, plus random fluctuations described by
$\bar{\tau}$. This decomposition leads to the estimate
$\mathrm{E}[\tau]\sim\frac{1}{\lambda}\log M+C$ for the expected value of the
relapse time, where the constant $C=\mathrm{E}[\bar{\tau}]$ can be obtained
from equation 5.

Figure 2: Relapse time densities
$f_{\tau}(t)=\frac{d}{dt}\,\mathrm{P}(\tau\leq t)$ computed from equation 4
for logistic and exponential primary growths and
$\nu=10^{-10},10^{-11},10^{-12},10^{-13},10^{-14}$ cells/day from left to
right. Using parameter estimates for colorectal cancer (see Table 2), the
logistic densities (dashed lines) converge to the corresponding exponential
ones as the initiation rate increases. Furthermore, in the exponential case
and for all the above values of $\nu$, the densities derived from equation 4
and their approximation obtained from equation 6 are indistinguishable.

###  2.3 Exponential population growth

Two commonly employed primary growth functions are the exponential and the
logistic ones (see e.g. [33]). These are given by $n(t)=e^{\delta t}$ and
$n(t)=\frac{Ke^{\delta t}}{K+e^{\delta t}-1}$, respectively, where $\delta$
denotes the primary tumor net growth rate and $K$ a carrying capacity. Relapse
time densities for these two growth types and different initiation rates are
shown in Figure 2. We observe that as $\nu$ increases, the logistic
distributions converge to the exponential ones (see Material and Methods for
more details). Moreover, for all our parameter estimates our model predicts
the same results with these two growth types. The reason is that the
metastases determining the time to relapse are initiated during the early
phase of tumor evolution which is almost exponential even for a logistic
growth. Therefore, from now on we will focus on exponentially growing primary
tumors. Exponential growth has the additional advantage that if only a portion
of primary cells can metastasize and their number is proportional to
$n(t)^{\gamma}$ (say only cells close to the surface of a spherical tumor for
$\gamma=2/3$), then this would be equivalent to changing the primary net
growth rate, that is using $n(t)=e^{\gamma\delta t}$.

Since the initiation rate $\nu$ is by far the slowest rate in our model, here
we study in detail the most relevant case, that is the small $\nu$ limit for
an exponentially growing tumor. The deterministic part of the relapse time
remains $\frac{1}{\lambda}\log M$, but interestingly the fluctuations
$\bar{\tau}$ are distributed as

$\bar{\tau}\sim\mathrm{Gumb}_{\mathrm{min}}\left(-\frac{1}{\delta}\log\frac{\nu(1-q)^{1-\frac{\delta}{\lambda}}\Gamma\left(\frac{\delta}{\lambda}\right)}{\lambda},-\frac{1}{\delta}\right)$
(6)

This Gumbel distribution describes the minimum of independent random variables
with exponential (left) tail, has two parameters $a$ and $b<0$ and
distribution $1-\exp(-e^{-\frac{x-a}{b}})$. Parameter $a$ describes a shift in
the distribution, and since $a\sim\log\nu$, it explains the equal spacing
between the densities in Figure 2 for logarithmically-spaced values of the
initiation rate. Also notice that these curves are left skewed, as it is
expected from the Gumbel for the minimum. On the other hand, the Gumbel for
the maximum - which describes the fluctuations of the time to detection
starting from a single initial cell - is right skewed. For small initiation
rates $\nu$ and large detectable sizes $M$, the mean relapse time is
approximately given by

$\mathrm{E}[\tau]\approx\frac{1}{\lambda}\log
M+\frac{1}{\delta}\log\frac{\delta}{\nu}+C\,,\qquad
C=-\frac{1}{\delta}\left(\log\frac{\delta(1-q)^{1-\frac{\delta}{\lambda}}\Gamma\left(\frac{\delta}{\lambda}\right)}{\lambda}+\gamma_{E}\right)$
(7)

where $\gamma_{E}\approx 0.5772$ denotes the Euler-Mascheroni constant. As
shown by Figure 3, this expression fits simulations even for relatively large
values of $\nu$ and small values of $M$.

Figure 3: Relapse time distribution for an exponentially growing primary tumor
using parameter estimates for colorectal cancer (see Table 2). Symbols
represent simulation results for a deterministic (diamonds) or a stochastic
primary growth (circles, see discussion at the end of the parameter estimation
section), while solid lines correspond to the theory in the small $\nu$ \-
large $M$ limit. On the left, each starred dot denotes the mean of $1000$
simulations, while lines represent the theoretical expectation given by
equation 7. These match the simulated means well for values of $\nu=10^{-6}$
or less. On the right, the relapse time densities derived from equation 14
yield a good approximation of the simulated data ($10000$ simulations per
curve) for $M=100$ or greater for both deterministic and stochastic primary
growth.

Equation 7 highlights a simple dependence of the mean relapse time $E[\tau]$
on $M$ and $\nu$. In Material and Methods we also compute the mean time to
detectability of the first established metastasis, $\mathrm{E}[\tau_{1}]$,
where $\tau_{1}=\sigma_{1}+\Theta_{1}$ is equal to the sum of the first
initiation time and the hitting time to $M$. Interestingly
$\mathrm{E}[\tau_{1}]$ has the same $M$ and $\nu$ dependence shown in equation
7, but the constant term $\tilde{C}$ is different. For example, using the
parameter estimates for colorectal cancer (see Table 2) we find $C\approx 250$
and $\tilde{C}\approx 309$. The reason for this difference is that even in the
small $\nu$ \- large $M$ limit, later established metastases can outrun the
earlier ones in reaching $M$ first.

While the mean relapse time $\tau$ shows logarithmic increase in terms of $M$
and $\delta/\nu$, its variance stays constant,
$\mathrm{Var}(\tau)=\pi^{2}/(6\delta^{2})$, see equation 15. Hence, due to the
slow logarithmic growth of the mean, the fluctuations of the relapse time stay
relevant even for large detection sizes and small mutation rates.

###  2.4 Relapse time with resection

Surgery is still the most common and effective type of treatment for solid
tumors, although often used in combination with other kind of therapies (see
e.g. [34]). However, how the time of resection affects prognosis, and in
particular the estimation of the time to relapse, is still unclear. In order
to investigate this question in a theoretical framework, we now embed surgery
in our model and study how it changes the distribution of the time $\tau$ to
relapse. Let us assume that at a given moment after detection a primary solid
tumor is surgically removed. This event can be mathematically implemented in
our model by considering a resection time $T$ such that $n(t)\equiv 0$ for
$t\geq T$. In particular, this implies that after $T$ no metastases can be
initiated. The number of metastases already established at resection is equal
to $K_{T}$, and their size distribution is given in [16]. The distribution of
the time $\tau$ to relapse can then be expressed exactly as in equation 4,
however here $\tau$ is not a proper random variable. In fact, as
$\int_{0}^{\infty}n(s)ds=\int_{0}^{T}n(s)ds<\infty$, there is a positive
probability that no metastasis will ever occur (notice that from this point of
view our framework can be seen as a cure model - see e.g. [35]) and in this
case we set $\tau=\infty$. The distribution of the relapse time conditioned on
at least one metastasis being established by resection is simply

$\,\mathrm{P}(\tau\leq t\mid K_{T}\geq 1)=\frac{\,\mathrm{P}(\tau\leq
t)}{\,\mathrm{P}(K_{T}\geq 1)}=\frac{1-e^{-b_{t}}}{1-e^{-a_{T}}}$ (8)

where we used that a relapsing metastasis had to be initiated before
resection, that is $\\{\tau\leq t\\}\subset\\{K_{T}\geq 1\\}$. This
conditional distribution for different resection times is depicted in Figure
4. In this and following figures, the resection time is shown at the bottom of
the figure, and the corresponding resection size $N=e^{\delta T}$ is shown on
the top. As $T\to 0$ all metastases have to be initiated close to time zero,
so the relapse time becomes the time to reach size $M$ from a single cell,
which has the Gumbel distribution for the maximum given by equation 2. If we
then increase the resection time, the conditional densities shift to the right
by the same amount. Finally, as $T\to\infty$ the relapse time distribution
converges to the case without resection

$\,\mathrm{P}(\tau\leq t\mid K_{T}\geq 1)\to\,\mathrm{P}(\tau\leq t)$

The fluctuations for the unconditional distribution follow a Gumbel type for
the minimum, as per equation 6. Hence, as time increases, the relapse time
distribution turns from a right-skewed Gumbel to a left-skewed Gumbel.

Note that the densities in Figure 4 become indistinguishable from the large
time limit as $\,\mathrm{P}(K_{T}\geq 1)$ approaches one. The reason is that
by this time metastases have probably already been initiated and one of the
early established ones is likely to relapse first. This suggests that only
early enough resection times change the behaviour of the model. For example in
the case of colorectal cancer, according to Figure 4, only resections of
tumors smaller than $10^{9}$ cells affect the time to recurrence.

Right skewed densities are often chosen to fit probability distributions
arising in survival analysis. This is due to the fact that most survival data
suffer from right censoring [36], where only a lower bound is known for data
points. Looking at the densities in Figure 4, though, we can see both left and
right skewed distributions. While a few survival datasets are negatively
skewed [37], cancer relapse times are typically right censored as a
consequence of limited follow-up and patients decease before relapse (see e.g.
[38]). However, our model does not take into account any of these events.
Furthermore [39] recently proposed a model for the estimation of screening
times for colorectal cancer based on the observation that some datasets suffer
from left censoring as well.

Figure 4: Relapse time densities $f_{\tau}(t\mid K_{T}\geq 1)$ conditioned on
at least one metastasis initiated by the time of resection $T$. For different
values of $T$, marked with ticks of corresponding colors, these densities are
computed by differentiating equation 8. As $T$ becomes larger, the probability
of metastases being established before resection (see equation 1) increases
and the conditional relapse time densities converge to the red limit one. Here
we have used parameter estimates for colorectal cancer (see Table 2),
$n(t)=e^{\delta t}$ and $7$ equally spaced resection times between $0.25$y and
$16.15$y. The curves for $T>15$y look identical to the limit density.

###  2.5 Metastasis classification

If the resection is successful and the primary tumor is completely removed,
the therapy can still fail due to the formation of metastases. For this
reason, it is common practice to start looking for detectable metastases
several weeks before the surgery. In this section we thus want to characterize
the metastases which are detectable at a given time and those which are not.

In general, for a fixed time $t$, the metastasizing process $(K_{s})_{0\leq
s\leq t}$ can be split into two independent Poisson processes $(S_{s})_{0\leq
s\leq t}$ and $(M_{s})_{0\leq s\leq t}$ describing the initiation of
metastases which reach size $M$ before or after $t$, respectively. Following
the same argument we used to derive the relapse time distribution, we obtain
that

$S_{t}\sim\mathrm{Pois}\left(b_{t}\right),\qquad
M_{t}\sim\mathrm{Pois}\left(c_{t}\right)$

where

$c_{t}=a_{t}-b_{t}$

In particular, we have that the events $\\{\tau>t\\}$ and $\\{S_{t}=0\\}$ are
equivalent. We also stress that the definitions above naturally extend to the
case of a primary resection, by simply redefining $n(t)$ to be zero after the
resection time $T$.

Now, despite an ongoing discussion on the following nomenclature (see e.g.
[40]), in the rest of the paper we will call a metastasis synchronous if it
reaches the detectable size $M$ before or up to the time of resection, and
metachronous otherwise (hereby the choice of notation $S_{t}$ and $M_{t}$).
These characterizations immediately allow us to estimate the probability of
some clinically relevant events. The probability of no synchronous metastases
is equal to

$\,\mathrm{P}(S_{T}=0)=\,\mathrm{P}(\tau>T)=e^{-b_{T}}$ (9)

Also, under this condition, relapse is not certain: the probability that at
least one metastasis was initiated given that there are no visible ones at
resection is

$P(K_{T}\geq 1\mid S_{T}=0)=\,\mathrm{P}(M_{T}\geq 1)=1-e^{-c_{T}}$ (10)

since $S_{T}$ and $M_{T}$ are independent. In next section we will study the
above and related quantities in greater detail.

##  3 Discussion

In this section we compare the predictions provided by our model with clinical
data collected for different cancer types. To this purpose, we first need to
estimate the parameter values for each of these cancer types.

###  3.1 Parameter estimation

The net growth rates of the primary and metastatic tumors, $\delta$ and
$\lambda$, are inferred from the corresponding tumor volume doubling times
(denoted $DT_{pt}$ and $DT_{m}$, respectively) as

$\delta=\log 2/DT_{pt}\,,\quad\lambda=\log 2/DT_{m}$

These times have been studied by many authors, starting from the influential
papers of [32, 41, 42]. Many authors still refer to these early works,
although in some case more recent estimates are available. Colorectal, breast
and lung cancers are the most frequently studied. Furthermore, more papers
focus on primary doubling times than on metastatic ones.

Similarly, the birth rate $\alpha$ is derived from the potential doubling time
$T_{pot}$, defined as the average time between cell divisions in the absence
of cell death [43, 44, 45]. In this case we simply use the estimation

$\alpha=1/T_{pot}$

Note that some authors (see e.g. [46]) define instead $T_{pot}$ as the tumor
doubling time in absence of cell death. While in this paper we employ the
former definition, the latter would simply yield a factor $\log 2$ in the
formula above.

As for the primary tumor size $N$ at resection, many studies report data on
the primary maximum diameter, allowing for ellipsoidal forms. However, given
the relatively small tumor volume and the wide interpatient variability, we
assume a spherical shape and estimate $d_{pt}$ from the corresponding typical
range. By also assuming $10^{9}$ cells per cm3, the primary size at resection
(expressed in number of cells) is thus estimated as $N=\frac{1}{6}\pi
d_{pt}^{3}10^{9}$.

Table 1 summarizes typical ranges of these quantities for five different
cancer types, together with the estimates we picked for our model and the
corresponding literature references. Difficulties in distinguishing between
primary and secondary tumors or in tracking down the primary origin of a
metastatic cancer could in principle affect some of these data, but the wide
range and multiple references reported reduce the potential impact of this
effect.

Table 1: Typical ranges of volume doubling times for the primary tumor ($DT_{pt}$) and metastasis ($DT_{m}$), tumor potential doubling time ($T_{pot}$) and tumor diameter at resection ($d_{pt}$) for breast, colorectal, headneck, lung and prostate cancer. Cancer type |  Parameter |  Typical range |  Estimate |  References
---|---|---|---|---
Breast |  $DT_{pt}$ (days) |  $103-353$ |  $210$ |  [47, 48, 49, 50, 51, 52]
$DT_{m}$ (days) |  $85-199$ |  $105$ |  [53, 54]
$T_{pot}$ (days) |  $8-35$ |  $15$ |  [44, 55]
$d_{pt}$ (cm) |  $1.4-3$ |  $2.5$cm |  [56, 57, 58]
Colorectal |  $DT_{pt}$ (days) |  $130-438$ |  $175$ |  [59, 60, 61]
$DT_{m}$ (days) |  $45-155$ |  $105$ |  [62, 54, 63, 64]
$T_{pot}$ (days) |  $3-4$ |  $4$ |  [65, 55]
$d_{pt}$ (cm) |  $3.5-5.1$ |  $4.5$cm |  [59, 66, 61, 67]
Headneck |  $DT_{pt}$ (days) |  $15-256$ |  $84$ |  [68, 69]
$DT_{m}$ (days) |  $9.5-320$ |  $56$ |  [70, 71]
$T_{pot}$ (days) |  $1-14$ |  $4$ |  [65, 72]
$d_{pt}$ (cm) |  $1.3-4$ |  $2.8$cm |  [73, 74]
Lung |  $DT_{pt}$ (days) |  $22-269$ |  $168$ |  [75, 76, 54, 77, 78]
$DT_{m}$ (days) |  $32-98$ |  $56$ |  [42, 79]
$T_{pot}$ (days) |  $2-17.5$ |  $2.5$ |  [75, 80]
$d_{pt}$ (cm) |  $1.7-4.1$ |  $2$cm |  [81, 82, 77]
Prostate |  $DT_{pt}$ (days) |  $36-1080$ |  $392$ |  [83, 84, 85]
$DT_{m}$ (days) |  $29-213$ |  $98$ |  [86, 85]
$T_{pot}$ (days) |  $15.2-97.8$ |  $34$ |  [87, 84]
$d_{pt}$ (cm) |  $0.1-2.9$ |  $1.2$cm |  [88, 89]

Notice that by estimating the rates $\lambda$ and $\alpha$ we also infer
values for the death rate $\beta=\alpha-\lambda$ and the extinction
probability $q=1-\lambda/\alpha$. For the two remaining parameters, namely the
initiation rate $\nu$ and the minimal detectable size of a metastasis $M$, we
use common estimates across different cancer types. Various studies report a
lowest detectable tumor diameter of $0.2$cm for different cancer types (see
e.g. [90, 91, 92]), corresponding to $M\approx 4.19\times 10^{6}$ cells.
Moreover, several papers argue that the first metastases are likely to be
established long before the detection of the primary tumor (see for example
[54] and the references therein). In particular, the review of the progression
model for metastases formation in [24] reports that dissemination starts when
the primary diameter is between $0.1$ and $0.4$cm. We thus consider the
primary tumor size at the expected time of the first metastasis initiation and
estimate it to be $e^{\delta\mathrm{E}[\sigma_{1}]}=10^{8}$ cells,
corresponding to a diameter of about $0.58$cm. Hence, by using the results in
Material and Methods, we set

$\nu=\frac{\delta e^{-\gamma_{E}}}{1-q}e^{-\delta\mathrm{E}[\sigma_{1}]}$

Finally, the carrying capacity for the logistic primary growth studied in
Figure 2 is set to $K=10^{12}$ [24, 93]. Overall, we thus found estimates for
the following input vector

$(DT_{pt},DT_{m},T_{pot},d_{pt},d_{m},e^{\delta\mathrm{E}[\sigma_{1}]})$

and used them as described above to derive values for our model parameters,
i.e.

$(\delta,\lambda,\nu,q,N,M)$

Such estimates are summarized in Table 2.

Before we proceed to study our model predictions, let us further discuss the
assumption of a deterministic primary tumor growth function. Firstly, as we
just showed, the only data we found to infer the rate of growth of a primary
tumor refer to doubling times, whose notion implicitly assumes an exponential
growth. For this reason we focus here on growth functions that (at least in
their early stages) show an exponential behaviour. Other growth functions, for
example $n(t)=ct^{3}$ for spherically growing tumours or
$n(t)=c^{\prime}t^{2}$ for tumors with active cells only around the surface,
could be studied when more data becomes available. Secondly, one could model
not just the metastases but also the primary tumor growth as a branching
process to account for further stochastic effects. However, due to the large
tumor size at resection a branching process model would predict an almost
perfect exponential growth around resection time. For this reason we set
$n(t)=e^{\delta t}$, which then determines the initial time $t=0$. Note that
the tumor is not initiated precisely at $t=0$, but that time is distributed
according to a Gumbel distribution, analogously to the results in the Single
type process section in Materials and methods. Since the initiation time is
not accessible experimentally anyway, for simplicity we use this above
definition for $t=0$. In order to justify the exponential deterministic
approximation for the primary size, we performed simulations where we modelled
the primary tumor as a branching process as well. We found that for initiation
rates of $\nu=10^{-5}$ or less (and all other parameters set for colorectal
cancer) the exponential approximation of the primary causes less than a few
percent error in the relapse time distribution (see Figure 3). The
relationship between stochastic and deterministic wild type populations has
been studied rigorously in [29].

Table 2: Parameter estimates for the primary net growth rate $\delta$, the metastatic net growth rate $\lambda$, the initiation rate $\nu$, the extinction probability $q$, the primary tumor size at resection $N$ and the minimal detectable size $M$. |  Breast |  Colorectal |  Headneck |  Lung |  Prostate
---|---|---|---|---|---
$\delta$ (cells/day) |  $0.0033$ |  $0.0040$ |  $0.0083$ |  $0.0041$ |  $0.0018$
$\lambda$ (cells/day) |  $0.0066$ |  $0.0066$ |  $0.0124$ |  $0.0124$ |  $0.0071$
$\nu$ (cells/day) |  $1.87\times 10^{-10}$ |  $8.42\times 10^{-10}$ |  $9.36\times 10^{-10}$ |  $7.49\times 10^{-10}$ |  $4.13\times 10^{-11}$
$q$ |  $0.9010$ |  $0.9736$ |  $0.9505$ |  $0.9691$ |  $0.7595$
$N$ (cells) |  $8.18\times 10^{9}$ |  $4.77\times 10^{10}$ |  $1.15\times 10^{10}$ |  $4.19\times 10^{9}$ |  $9.05\times 10^{8}$
$M$ (cells) |  $4.19\times 10^{6}$ |  $4.19\times 10^{6}$ |  $4.19\times 10^{6}$ |  $4.19\times 10^{6}$ |  $4.19\times 10^{6}$

###  3.2 Model predictions

Now that we have estimated the parameters of our model in Table 2, we are in
position to study its predictions and compare them to clinical data.

Figure 5: Probability of extant metastases, $\,\mathrm{P}(K_{T}\geq 1)$,
dashed curve computed from equation 1, and synchronous metastases,
$\,\mathrm{P}(S_{T}\geq 1)$, solid curve computed from equation 4. These
probabilities are plotted as functions of the resection time $T$ for five
different cancer types. The primary tumor size at resection is $N=e^{\delta
T}$ and thus depends on the primary net growth rate. These resection sizes are
discussed in Table 3. For each cancer type, the shaded areas highlight
resection time intervals leading to a probability higher than $85\%$ of
established and all undetectable metastases. Using the parameter estimates
from Table 2, the widths of these intervals are $3.41,3.17,1.92,0.94,1.19$
years for breast, colorectal, headneck, lung and prostate cancer respectively.

Let us start by analyzing the simplest predictions of the model, which are
about the presence of synchronous and metachronous metastases. Figure 5 shows
the probability of initiated metastasis by resection $\,\mathrm{P}(K_{T}\geq
1)=1-e^{-a_{T}}$ (equation 1), and the probability of visible metastasis by
resection $\,\mathrm{P}(S_{T}\geq 1)=1-e^{-b_{T}}$ (equation 4) as functions
of the resection time $T$, for five different cancer types. From the figure we
notice that, obviously, the probability of having initiated metastasis is
always higher than the probability of having visible metastasis at resection.
For all five cancer types considered, one or more metastases have likely been
initiated by the time the primary tumor reaches about $8.2\times 10^{8}$ cells
(diameter $1.16$cm). While this value is similar across different primary
types (as a consequence of the parameters estimation procedure), the results
for the probability of synchronous metastases vary widely. For breast,
colorectal, headneck, lung and prostate cancer, Table 3 reports primary tumor
sizes at which synchronous metastases might start to appear and are likely to
be present, respectively (expressed both in terms of number of cells and tumor
diameter). By comparing these values to typical resection sizes in Table 1, we
find that detecting metastases at resection is very likely for lung and
prostate cancer and rare for headneck primary tumors.

Table 3: Resection sizes of the primary tumor which yield a $1\%$ and $99\%$ probability of synchronous metastases, respectively. For each cancer type considered, these sizes are computed with the parameter values in Table 2 and expressed both in terms of number of cells, $N$, and tumor diameter, $d$. |  |  Breast |  Colorectal |  Headneck |  Lung |  Prostate
---|---|---|---|---|---|---
$\,\mathrm{P}(S_{T}\geq 1)>0.01$ |  $N$ |  $1.32\times 10^{9}$ |  $2.13\times 10^{9}$ |  $7.03\times 10^{9}$ |  $1.03\times 10^{8}$ |  $6.27\times 10^{7}$
$d$ |  $1.36$ |  $1.60$ |  $2.38$ |  $0.58$ |  $0.49$
$\,\mathrm{P}(S_{T}\geq 1)>0.99$ |  $N$ |  $6.03\times 10^{11}$ |  $9.88\times 10^{11}$ |  $3.22\times 10^{12}$ |  $4.65\times 10^{10}$ |  $2.89\times 10^{10}$
$d$ |  $10.48$ |  $12.36$ |  $18.32$ |  $4.46$ |  $3.81$

One of the most challenging scenarios for the development of an effective
treatment is when there are only undetectable metastases present. In our
framework this scenario corresponds to the event

$U_{T}:=\\{K_{T}\geq 1,S_{T}=0\\}$ (11)

which has probability (see equations 10 and 9)

$\begin{split}\,\mathrm{P}(U_{T})&=\,\mathrm{P}(M_{T}\geq
1,S_{T}=0)=\,\mathrm{P}(M_{T}\geq 1)\,\mathrm{P}(S_{T}=0)\\\
&=e^{-b_{T}}-e^{-a_{T}}=\,\mathrm{P}(K_{T}\geq 1)-\,\mathrm{P}(S_{T}\geq
1)\end{split}$ (12)

Because of the last identity, the probability of established and all
metachronous metastases can be read out from Figure 5 as the difference of the
two curves. There, the shaded areas highlight intervals of resection times
yielding $\,\mathrm{P}(U_{T})>85\%$. These intervals, often referred to as
high-risk period [94], are especially wide for breast, colorectal and headneck
cancers. The reason is that these cancer types have a lower ratio of
metastatic over primary net growth rates, so metastases take longer to grow to
visible size. Hence, although for these cancer types metastases grow slower,
which improves prognosis, they stay undetectable for longer, which poses a
challenge for diagnosis. The estimated resection sizes given in Table 1 fall
within or close to these ranges ($\,\mathrm{P}(U_{T})$ equal to
$93.87\%,79.83\%,98.35\%,66.04\%$ and $85.85\%$ for the five primary tumor
types studied, respectively). In general, by assuming that the primary tumor
diameter at resection fits a normal distribution (with mean computed as the
mean of $d_{pt}$ and variance set so that 95% of the observations belong to
its typical range, see Table 1) we estimate that resections for breast and
headneck cancers fall in the high-risk window $99.8\%$ and $99.58\%$ of the
times respectively, followed by colorectal ($24.67\%$), prostate ($13.41\%$)
and lung ($0.69\%$) cancers.

Figure 6: Probability of established and all metachronous metastases -
$P(U_{T})$, as given by equation 12 \- plotted as a function of $T$ and
$\delta$ (left panel) and of $T$ and $\nu$ (right panel). The parameter
estimates used are those for colorectal cancer reported in Table 2. The plots
show that the width of the high-risk interval - the range of resection times
such that $\,\mathrm{P}(U_{T})$ is high - stays roughly constant for most
parameter values. This width (about $3$ years) shrinks only for metastases
growing significantly faster than the primary tumor that initiated them.

In order to check how robust the presence of a wide high-risk interval is, we
plotted in Figure 6 the probability of having only undetectable metastasis at
detection, $P(U_{T})$, for different values of the primary net growth rate
$\delta$ and of the initiation rate $\nu$. Other parameters are taken for
colorectal cancer. The width of the high-risk interval is constant with
respect to $\nu$, and shrinks only as the ratio between the primary and
metastatic net growth rate becomes very small. The same qualitative behaviour
can be obtained with the parameter estimates for the other cancer types. As
most metastases grow up to two times faster than the primary tumor they
originated from [24], our model suggests that for a wide choice of parameters
there is a substantial range of resection sizes that lead to a high
probability of established and all undetectable metastases.

Next, we ask how such a probability, $\,\mathrm{P}(U_{T})$, influences the
time to cancer recurrence. The conditional distribution of the relapse time
$\tau$ becomes

$\,\mathrm{P}(\tau\leq t\mid U_{T})=\frac{\,\mathrm{P}(T<\tau\leq t\mid
K_{T}\geq 1)}{\,\mathrm{P}(T<\tau\mid K_{T}\geq
1)}=\frac{e^{-b_{T}}-e^{-b_{t}}}{e^{-b_{T}}-e^{-a_{T}}}$

for $t\geq T$, where we used the definition of $U_{T}$ (see equation 11) and
equation 8. From this distribution we compute the expected relapse time
measured from resection and conditioned on $U_{T}$, $\mathrm{E}[\tau-T\mid
U_{T}]$. This expectation and the probability $\,\mathrm{P}(U_{T})$ are
plotted in Figure 7. We see that for resection sizes smaller than $10^{8}$
cells the relapse occurs on average between $4$ and $5$ years after resection,
independently of the primary size. For resection sizes around $10^{8}$ cells,
undetectable metastases become likely to be present and $\mathrm{E}[\tau-T\mid
U_{T}]$ starts to decrease with tumor size. At about $19$ years the
probability of only undetectable metastases present and the conditional mean
relapse time both approach zero.

Let us stress that while some clinical studies report data on the whole
distribution of recurrence times, these are usually measured from a varying
time of surgery, which corresponds to different primary tumor resection sizes.
Therefore, unless the distribution of relapse times is reported together with
the corresponding resection sizes, we cannot compare it directly to the
predictions of our model. However, we expect the variability of primary sizes
at resection to average out across large cohorts of patients, which is why we
analyzed the expected value of the time to recurrence.

Figure 7: Expected relapse time measured from resection, conditioned on extant
but all undetectable metastases (blue curve). The dashed line and the light
blue shaded area show $\,\mathrm{P}(U_{T})$ and how spread is the conditional
relapse time distribution, respectively. The parameter estimates used are
those for colorectal cancer reported in Table 2. For resection times close to
zero this conditional expectation coincides with that of the Gumbel
distribution given by equation 14, at about $5$ years. As $T$ starts to
increase $\mathrm{E}[\tau-T\mid U_{T}]$ reflects the convergence highlighted
for Figure 4, first slightly decreasing and then staying constant around $4.4$
years. Finally, when the resection time falls into the high-risk window, the
expected relapse time drops to zero. This suggests that the bigger the primary
tumor size is at resection, the faster relapse will occur. Figure 8: Disease-
free curves for different resection times. The earlier the primary tumor is
resected the higher is the probability that no metastases will arise, or cure
probability, represented by the value of the final plateaus. The resection
times are chosen so that $\,\mathrm{P}(K_{T}=0)=0.75,0.6,0.45,0.3,0.15,0.001$
respectively. With the parameter estimates for colorectal cancer (see Table 2)
these times range from $12.28$ to $14.48$ years, corresponding to sizes
between $5.12\times 10^{7}$ and $1.23\times 10^{9}$ cells (diameter
$0.46-1.33$cm), respectively.

Using the values from Table 2 we tested our model by computing the probability
of synchronous metastases and the mean relapse time conditioned on established
but all undetectable metastases. The predictions from our model, typical
ranges and references for each cancer type considered are summarized in Table
4. Notice that our predictions for the mean relapse time fall on the lower end
of the respective typical ranges. This is expected since we compute the time
to recurrence $\tau$ based on the minimal detectable size $M$, while in
practice metastases are often detected only at larger sizes.

In general, for different cancer types it is observed that metastases can grow
up to $2$ times faster than the primary tumor they originated from [24],
although values as high as $4$ has been proposed [95]. Our estimates fall
within this range ($\delta/\lambda=4$ for prostate cancer, $3$ for lung and
between $1.5$ and $2$ for the others).

As per the time interval from primary onset to surgery, the typical range is
$15-25$ years [43]. The high variability in our estimates of $DT_{pt}$ make
$T$ fall outside that range for headneck ($T=7.69$y), lung ($T=14.71$y) and
prostate ($T=32$y) cancers, classifying the first two as fast growing tumors
and the latter as a slow growing one. The singular features that the model
predicts for prostate cancer are in accordance with clinical studies (see e.g.
[96, 86]).

Table 4: Typical ranges of $\,\mathrm{P}(S_{T}\geq 1)$ and $\mathrm{E}[\tau-T\mid U_{T}]$, predicted value from the model and literature references for each cancer type. Cancer type |  Output |  Range from clinical data |  Theoretical prediction |  Reference
---|---|---|---|---
Breast |  $P(S_{T}\geq 1)$ ($\%$) |  $5-10$ |  $6.13$ |  [97, 98, 99]
$\mathrm{E}[\tau-T\mid U_{T}]$ (days) |  $590-1022$ |  $725$ |  [100, 101, 102]
Colorectal |  $P(S_{T}\geq 1)$ ($\%$) |  $15-25$ |  $20.17$ |  [103, 104, 105, 106, 107, 99]
$\mathrm{E}[\tau-T\mid U_{T}]$ (days) |  $353-760$ |  $356$ |  [108, 109, 106, 110, 107]
Headneck |  $P(S_{T}\geq 1)$ ($\%$) |  $1-16.8$ |  $1.65$ |  [111, 112]
$\mathrm{E}[\tau-T\mid U_{T}]$ (days) |  $219-623$ |  $435$ |  [113, 114, 115]
Lung |  $P(S_{T}\geq 1)$ ($\%$) |  $30-55.39$ |  $33.96$ |  [99, 116]
$\mathrm{E}[\tau-T\mid U_{T}]$ (days) |  $210-602$ |  $249$ |  [117, 118, 119]
Prostate |  $P(S_{T}\geq 1)$ ($\%$) |  $10-34$ |  $13.53$ |  [120, 121, 122]
$\mathrm{E}[\tau-T\mid U_{T}]$ (days) |  $730-1131$ |  $969$ |  [123, 124]

The last trait of cancer recurrence that we are going to examine is disease-
free rates. These generally correspond to the survival function of the relapse
time, $\,\mathrm{P}(\tau>t)$. However, following the previous discussion we
will condition this probability on no synchronous metastases, obtaining

$\,\mathrm{P}(\tau\leq t\mid S_{T}=0)=1-e^{-(b_{t}-b_{T})}$ (13)

for $t\geq T$. In this case we do not observe any convergence to the density
without resection, because if $T\to 0$ then no metastasis can be initiated and
if $T\to\infty$ the condition $S_{T}=0$ pushes the relapse time to infinity.
Let us also stress that our model does not provide information on survival
rates, as no modelling of the time to decease is incorporated. Furthermore,
notice that $P(\tau>t)$ yields a good description of the disease-free rates in
terms of metastases detectability, but not necessarily with respect to cancer
symptomaticity.

The relapse time distribution in case of no synchronous metastasis,
$\,\mathrm{P}(\tau>t\mid\tau>T)$, for different resection times is shown in
Figure 8, studying again the case of colorectal cancer. As we are not
conditioning on at least one metastasis being initiated, there is always a
positive probability that relapse will not occur, that is $\tau=\infty$. The
resection times are thus chosen so to yield cure probabilities -
$\,\mathrm{P}(K_{T}=0)$, corresponding to the final plateaus - equal to
$0.75,0.6,0.45,0.3,0.15$ and $0.001$, respectively. These times span across a
total range of about $2.2$ years. Furthermore, excluding the latest resection
time considered, the difference between two consecutive of these $T$ values is
between $0.28$ and $0.4$ years. Hence, our model suggests that delays of the
order of months in the time of primary resection can lead to a significant
decrease in the cure probability.

Figure 9: Probability of the first metastasis being initiated during surgery
delay. This probability $\,\mathrm{P}(K_{T+\Delta T}\geq 1,K_{T}=0)$ \- where
$T$ is the set time of resection and $\Delta T$ the surgery delay - is plotted
as a function of the resection size $N=e^{\delta T}$ (x-axis) and of the delay
$\Delta T$ (y-axis). With the parameter estimates for colorectal cancer (see
Table 2) we see that if a primary tumor is resected at a critical size (around
$2\times 10^{8}$ cells, diameter $\approx 0.725$cm), surgery delays of 2-3
months can decrease the cure probability of more than $10\%$.

To quantify more precisely the implications of surgery delays, we study the
probability that the first metastasis is initiated in the time interval
$(T,T+\Delta T)$

$\,\mathrm{P}(K_{T+\Delta T}\geq 1,K_{T}=0)=\,\mathrm{P}(K_{T+\Delta T}\geq
1)-\,\mathrm{P}(K_{T}\geq 1).$

This probability is depicted in Figure 9 using our parameter estimates for
colorectal cancer. We see in the figure that there is a middle range of
resection sizes where the recurrence probability can be significantly affected
by small surgery delays. For colorectal cancer we estimate that if the primary
resection is originally planned for a tumor of diameter between $0.44$ and
$0.9$ centimeters ($4.39\times 10^{7}$ and $3.89\times 10^{8}$ cells,
respectively), then a surgery delay of 60 days would decrease the cure
probability by $5-9\%$. Conversely, tumors smaller than this critical range
are less likely to metastasize during the surgery delay, while larger tumors
likely metastasized already, so that the effect of surgery delay for these
sizes is smaller.

###  3.3 Conclusions

We introduced a model of metastasis formation where metastases are initiated
at a time dependent rate, in the simplest case proportional to the size of a
growing primary tumor. All initiated metastases then evolve as independent
supercritical branching processes. Parameters of the model were estimated for
five different cancer types from the clinical literature. We studied the
relapse time $\tau$, that is the earliest time when any of the metastases
becomes detectable. We obtained the distribution of $\tau$ for a general
primary tumor growth and focused in particular on logistic and exponential
growth functions. For clinically relevant initiation rates the metastases
which relapse first are typically initiated in the early phase of the primary
tumor development, which is exponential for both growth functions considered.
Hence the distributions of $\tau$ for exponential and logistic primary growths
are practically identical unless the initiation rate is unrealistically small
($\nu\approx 10^{-13}$ or smaller) and we can thus exploit the much simpler
formulas for the exponentially growing tumor.

We model the resection of the primary tumor by introducing a cut-off for the
growth function $n(t)$. If metastases are likely already established at
surgery, their time of relapse is not influenced by the resection timing. We
categorized all metastases into synchronous and metachronous and computed
corresponding occurrence probabilities. With our estimated parameters we found
that the probability of synchronous metastases and the mean relapse time after
resection falls in the typical clinical range for all five different cancer
types we study.

A challenging scenario for treatment is that of patients with established but
all undetectable metastases. For all five cancer types we considered, the
probability of this event is high within a significant range of resection
sizes. Unfortunately, the typical size of a resected tumor falls in or near
this range for all cancer types. We found that relatively small delays in
these resection times can cause significant decrease in the cure probability.
Within our model, surgery only prevents recurrence if it is done before the
onset of the first surviving metastases.

The parameter estimates summarized in Table 2 yield realistic predictions for
several quantities of clinical interest. Although in principle we can explore
our model predictions across the whole range of parameters, this would often
lead to unrealistic outcomes. In this sense the quantitative predictions of
our model are quite sensitive to the parameter values, but we have been able
to find a combination of parameters that yields realistic results. On the
other hand, the qualitative features of our model are more robust to parameter
changes, as demonstrated for example in Figure 6.

In particular, our estimate for the initiation rate of metastases, $\nu$, is
based on the assumption of early dissemination at primary size,
$e^{\delta\mathrm{E}[\sigma_{1}]}=10^{8}$ cells. However, since the metastatic
net growth rate $\delta$ and extinction probability $q$ are estimated
independently from data on tumor volume doubling times (see Table 1), by
changing the early dissemination assumption our model predictions could fall
outside their typical ranges. For example, for colorectal cancer, assuming
$e^{\delta\mathrm{E}[\sigma_{1}]}=10^{9}$ cells would lead to the unrealistic
values $\mathrm{E}[\tau-T\mid U_{T}]=836$ days and $P(S_{T}\geq 1)=2.23\%$.
Thus, indirectly, our model supports the hypothesis of early metastatic
dissemination.

Note that in this paper we focused on the presence or absence of synchronous
or metachronous metastasis at resection as these events determine if there is
ever a relapse ($K_{T}\geq 1$) or if relapse has already occurred by resection
($S_{T}\geq 1$). Our model also provides estimates for the number and sizes of
metastases at resection, but these are less relevant for the study of the time
to cancer relapse, and have already been studied in detail in [16, 29]. A
general feature the model predicts is that the cumulative distribution of
metastases sizes at resection has a power law tail with exponent
$\delta/\lambda$. This power law tail was observed in [16] using data on 21
patients with colorectal cancer from [13], and the exponent was found to
typically be in the range $0.3-0.8$. Our estimate $\delta/\lambda=0.61$ falls
in this range, supporting our parameter inference. The paper above also
reports data on the number of visible metastases at surgery. In our model, for
a given primary resection size $e^{\delta T}$, this number is a Poisson random
variable $S_{T}$. However, since the primary tumor sizes are not published and
likely different for all patients in the data, we could not use this quantity
reliably for our parameter estimation. For example, if we infer the initiation
rate $\nu$ for colorectal cancer from the probability of visible metastasis at
resection (given by $P(S_{T}\geq 1)=1-e^{-b_{T}}$ and with estimate $0.2$ from
the data reported in Table 4) we would get essentially the same estimate as in
Table 2. But by using instead the mean number of visible metastases (expressed
as $E(S_{T}|S_{T}\geq 1)=b_{T}/(1-e^{-b_{T}})$, with estimated lower bound
$1.4$ from [107]), we would infer a 3 times greater estimate for $\nu$. Again,
a possible cause of this discrepancy lies in the different resection sizes for
patients which we have no data for.

Metastases are seeded and establish colonies via a specific and complex
process called metastatic cascade (for details see e.g. [125]). Since this is
known to be a multi-stage process, some authors (see for example [126, 6, 127]
and references therein) have described metastases initiation through two-type
stochastic models, where a cell needs to gain the ability to metastasize
before it can establish a new metastatic lesion. We did not choose that route
for several reasons: (i) the precise details of how and when cells reach this
ability are not clear [43, 128], (ii) in our model we can think of $n(t)$ as
the number of cells which can metastasize and so tailor the two approaches,
and (iii) if we assume that an acquired metastatic ability lowers the primary
net growth rate and that the seeding rate is sufficiently small (at most
$10^{-4}$ according to simulations), a branching process model would predict
the same exponential growth for the cells with this ability [29, 129], and
hence this would only change the estimate of the initiation rate in our model.

We did not include into our model a mechanism for metastasis seeding other
metastasis, although this phenomena has been observed in clinical studies
[130]. The main reason for this omission was the lack of reliable data for the
estimation of the secondary seeding rate. By assuming the same primary and
secondary seeding rates, however, we would expect metastases to initiate
secondary ones when they reach around $10^{8}$ cells, at which size they are
already detectable. Hence, by considering this scenario our predictions for
the time to cancer relapse would not change.

We aim to compare our model in the future to data where relapse times are
given jointly with primary tumor sizes at resection. Tumor size is of course
not the only relevant factor in predicting relapse times, so the model should
be extended to involve other features like a measure of malignancy, possibly
as in [131]. Many of the parameters of the model can differ between patients,
and also between each metastasis. Therefore, including a probability
distribution for the parameters could also make our model more realistic,
provided that such distributions can be estimated from data. Other possible
extensions could include interactions among metastatic cells and among
metastatic lesions, effects of the immune system, allowing metastases to seed
other metastases, and providing an estimate for the fraction of cells which
can metastasize, perhaps through modelling angiogenesis.

##  4 Materials and methods

In this section we provide more details about the mathematical foundations of
our model.

###  4.1 Single type process

Let $(Z_{t})_{t\geq 0}$ be a birth-death branching process, i.e. a Markov
chain on non-negative integers with transition rates

$\displaystyle i\mapsto\begin{cases*}i+1&\text{ rate } $i\alpha$\\\ i-1&\text{
rate } $i\beta$\end{cases*}$

The two positive constants $\alpha$ and $\beta$ are called birth and death
rate, respectively. In our model we employ this process to describe the
evolution of each metastasis. We assume that all metastases have the same
birth and death rate and that they are supercritical, that is they have
positive net growth rate $\lambda=\alpha-\beta>0$. Moreover, since we only
want to model surviving metastasis, we condition on the eventual survival of
the process, that is on the event
$\Omega_{\infty}=\\{\omega:Z_{t}(\omega)>0\text{ for all }t\geq 0\\}$. The
probability of such event is equal to $\,\mathrm{P}(\Omega_{\infty})=1-q$,
where $q=\beta/\alpha$ [31].

We define the first passage time to size $M$ as
$T_{M}:=\inf\\{t>0:Z_{t}=M\\}$. A well known property of branching processes
is that $e^{-\lambda t}Z_{t}\to W$ almost surely as $t\to\infty$, and
conditioned on survival and a single initial cell
$W\sim\mathrm{Expo}(\lambda/\alpha)$ [31]. Since $W$ and $T_{M}$ are connected
by $\lim\limits_{M\to\infty}Me^{-\lambda T_{M}}=W$, an immediate consequence
is that

$\,\mathrm{P}(T_{M}\leq t\mid\Omega_{\infty}\cap\\{Z_{0}=1\\})\sim
e^{-(1-q)Me^{-\lambda t}}\equiv G(t)\qquad\text{ as }M\to\infty$ (14)

Early derivations of this result already appear in [132, 133]. Interestingly,
$T_{M}$ follows the Gumbel distribution
$\mathrm{Gumb}_{\mathrm{max}}\left(\frac{\log
M(1-q)}{\lambda},\frac{1}{\lambda}\right)$, where

$Y\sim\mathrm{Gumb}_{\mathrm{max}}(a,b)\quad\Longleftrightarrow\quad\,\mathrm{P}(Y\leq
y)=e^{-e^{-\frac{y-a}{b}}},\quad a\in\mathbb{R},\;b>0$

The Gumbel type is an extreme value distribution. If $M_{n}$ denotes the
maximum of $n$ IID random variables $X_{i}$, the Gumbel distribution above
generally describes the limit of $M_{n}$ as $n\to\infty$, when $X_{i}$ have an
exponential (right) tail. A similar definition can be given for the reverse
Gumbel distribution, i.e. the limit of minimum of IID random variables with an
exponential (left) tail

$Y\sim\mathrm{Gumb}_{\mathrm{min}}(a,b)\quad\Longleftrightarrow\quad\,\mathrm{P}(Y\leq
y)=1-e^{-e^{-\frac{y-a}{b}}},\quad a\in\mathbb{R},\;b<0$

For both of these distributions we have

$Y\sim\mathrm{Gumb}(a,b)\quad\Longrightarrow\quad\mathrm{E}[Y]=a+b\gamma_{E},\quad\mathrm{Var}(Y)=\frac{\pi^{2}}{6}b^{2}$
(15)

where $\gamma_{E}\approx 0.5772$ denotes the Euler-Mascheroni constant. Hence
the mean hitting time to $M$ cells grows logarithmically with $M$, while its
variance remains constant

$\mathrm{E}[T_{M}]=\frac{\log
M(1-q)+\gamma_{E}}{\lambda},\quad\mathrm{Var}\left(T_{M}\right)=\frac{\pi^{2}}{6\lambda^{2}}$

Thus, for sizes $M\approx\alpha/\lambda$ the standard deviation is
approximately equal to the mean, but since the mean only grows logarithmically
with $M$, fluctuations of $T_{M}$ stay relevant even for much larger values of
$M$.

###  4.2 Scaled relapse time distribution

In Results we derive the general expression for the relapse time distribution,
whose full expression is obtained by combining equations 14, 3 and 4. Here we
show how to scale the detectable size $M$ out of this expression, so to split
the distribution into a deterministic part and a stochastic term. Let us focus
on the integral

$\int_{0}^{t}n(s)G(t-s)ds=\int_{0}^{t}n(s)e^{-(1-q)Me^{-\lambda(t-s)}}ds$

and apply the change of variables $z:=t-s-\frac{1}{\lambda}\log M$ to find

$\int_{-\frac{1}{\lambda}\log M}^{t-\frac{1}{\lambda}\log
M}n\left(t-z-\frac{1}{\lambda}\log M\right)e^{-(1-q)e^{-\lambda z}}dz$

By plugging this expression back into equation 4, at time
$t+\frac{1}{\lambda}\log M$ we get

$\,\mathrm{P}\left(\tau-\frac{1}{\lambda}\log M\leq
t\right)=1-\exp\left(-\nu(1-q)\int_{-\frac{1}{\lambda}\log
M}^{t}n(t-z)e^{-(1-q)e^{-\lambda z}}dz\right)$

Hence, as $M$ tends to infinity we obtain

$\tau-\frac{1}{\lambda}\log M\xrightarrow[M\to\infty]{d}\bar{\tau}$

where

$\,\mathrm{P}(\bar{\tau}\leq
t)=1-e^{-\nu(1-q)\int_{-\infty}^{t}n(t-s)e^{-(1-q)e^{-\lambda s}}ds}$

From the last two equations we also see that asymptotically as $M\to\infty$

$\mathrm{E}[\tau]\sim\frac{1}{\lambda}\log M+C,\qquad
C=\mathrm{E}[\bar{\tau}]$ (16)

###  4.3 Explicit results for exponential primary growth

Two commonly employed growth functions for primary tumors are the exponential
$n_{e}(t)=e^{\delta t}$ and the logistic $n_{l}(t)=\frac{Ke^{\delta
t}}{K+e^{\delta t}-1}$ ones (see e.g. [33]). A logistic growth implies that
the primary tumor has a carrying capacity $K$. During the first stages of its
development $n_{l}(t)$ follows the same exponential trajectory of $n_{e}(t)$
and then approaches a constant as it gets closer to size $K$. As the carrying
capacity is typically large, this slowdown for $n_{l}(t)$ happens around
$\hat{t}=\log(K)/\delta$. The differences between the results provided by
these two growths functions thus depend on the probability of metastases being
initiated by time $\hat{t}$, i.e. $\,\mathrm{P}(K_{\hat{t}}\geq 1)\approx
1-e^{-\frac{\nu(1-q)K}{\delta}}$. Hence, if

$\frac{\nu(1-q)K}{\delta}\gg 1$ (17)

metastases are likely established in the first stages of the primary growth,
i.e. when $n_{l}(t)\approx n_{e}(t)$. Otherwise, metastases are initiated late
in the primary evolution, when the two growth functions are substantially
different. This feature is visualized in Figure 2, where $\tau$ densities for
a logistic growth are shown to converge to the exponential ones as $\nu$
increases and the other parameters are fixed.

Using the parameter values from Table 2, however, we observe that the
condition in equation 17 is satisfied for all cancer types considered. In
other words, our estimates for $\nu,q,K$ and $\delta$ yield no difference
between exponential and logistic growth functions. In light of this, we study
in greater detail the results obtained with $n_{e}(t)$.

####  4.3.1 Scaled relapse time

When $n(t)=n_{e}(t)=e^{\delta t}$, the relapse time distribution has an
expression in terms of special functions. To show this, let us consider the
distribution of the scaled relapse time $\overline{\tau}$ as given by equation
5 and focus on the integral

$\int_{-\infty}^{t}n(t-s)e^{-(1-q)e^{-\lambda s}}ds=e^{\delta
t}\int_{-\infty}^{t}e^{-(1-q)e^{-\lambda s}-\delta s}ds$

This can be equivalently written as

$e^{\delta
t}\int_{-\infty}^{t}\frac{1}{(1-q)^{\delta/\lambda}}\left[(1-q)e^{-\lambda
s}\right]^{\delta/\lambda}e^{-(1-q)e^{-\lambda s}}ds$

The last expression then suggests the change of variable $x=(1-q)e^{-\lambda
s}$, which leads to

$\frac{e^{\delta t}}{\lambda(1-q)^{\delta/\lambda}}\int_{(1-q)e^{-\lambda
t}}^{\infty}x^{\frac{\delta}{\lambda}-1}e^{-x}dx=\frac{e^{\delta
t}}{\lambda(1-q)^{\delta/\lambda}}\Gamma\left(\frac{\delta}{\lambda},(1-q)e^{-\lambda
t}\right)$

and where $\Gamma$ denotes the incomplete upper gamma function
$\Gamma(a,t)=\int_{t}^{\infty}x^{a-1}e^{-x}dx$. The scaled relapse time
distribution for $n(t)=e^{\delta t}$ is thus given by

$\,\mathrm{P}(\bar{\tau}\leq
t)=1-\exp\left(-\frac{\nu(1-q)^{1-\frac{\delta}{\lambda}}e^{\delta
t}}{\lambda}\Gamma\left(\frac{\delta}{\lambda},(1-q)e^{-\lambda
t}\right)\right)$ (18)

Since $\Gamma(1,t)=e^{-t}$, for $\lambda=\delta$ this simplifies to

$\,\mathrm{P}(\bar{\tau}\leq
t)=1-\exp\left(-\frac{\nu}{\lambda}e^{-(1-q)e^{-\lambda t}+\lambda t}\right)$

####  4.3.2 Small initiation limit

While the initiation rate can vary significantly across different cancer
types, $\nu$ is typically orders of magnitude smaller then all other
parameters. Hence, we now investigate $\tau$ distribution in the $\nu\to 0$
limit. Let us first consider the result given by equation 18 for the scaled
time to recurrence $\overline{\tau}$ and write it as

$\,\mathrm{P}(\bar{\tau}\leq
t)=1-\exp\left(-\frac{\nu(1-q)^{1-\frac{\delta}{\lambda}}\Gamma\left(\frac{\delta}{\lambda}\right)e^{\delta
t}}{\lambda}\right)\exp\left(\frac{\nu(1-q)^{1-\frac{\delta}{\lambda}}\phi(t)e^{\delta
t}}{\lambda}\right)$ (19)

where

$\phi(t)=\int_{0}^{(1-q)e^{-\lambda
t}}s^{\frac{\delta}{\lambda}-1}e^{-s}ds<\int_{0}^{(1-q)e^{-\lambda
t}}s^{\frac{\delta}{\lambda}-1}ds=\frac{\lambda}{\delta}(1-q)^{\frac{\delta}{\lambda}}e^{-\delta
t}$

Notice that the second exponential factor in equation 19 is bounded below by
$1$ and above by $e^{\frac{\nu(1-q)}{\delta}}$ for all $t\geq 0$. Therefore,
as $\nu\to 0$, the distribution of $\bar{\tau}$ asymptotically converges to

$\,\mathrm{P}(\bar{\tau}\leq t)\sim
1-\exp\left(-\frac{\nu(1-q)^{1-\frac{\delta}{\lambda}}\Gamma\left(\frac{\delta}{\lambda}\right)e^{\delta
t}}{\lambda}\right)$

Equivalently, for small initiation rates the scaled relapse time $\bar{\tau}$
asymptotically follows a Gumbel distribution for the minimum,
$\bar{\tau}\sim\mathrm{Gumb}_{\mathrm{min}}\left(-\frac{1}{\delta}\log\frac{\nu(1-q)^{1-\frac{\delta}{\lambda}}\Gamma\left(\frac{\delta}{\lambda}\right)}{\lambda},-\frac{1}{\delta}\right)$.

####  4.3.3 Mean relapse time

By combining the last result with equations 15 and 16 we find that

$\mathrm{E}[\tau]\approx\frac{1}{\lambda}\log
M+\frac{1}{\delta}\log\frac{\delta}{\nu}+C$

where
$C=-\frac{1}{\delta}\left(\log\frac{\delta(1-q)^{1-\frac{\delta}{\lambda}}\Gamma\left(\frac{\delta}{\lambda}\right)}{\lambda}+\gamma_{E}\right)$.
Intuitively, the time to relapse is likely to be determined by one of the
first established metastases. Given the simple dependence of
$\mathrm{E}[\tau]$ on $M$ and $\nu$, we now compare it with the mean time to
detectability of the first metastasis, $\mathrm{E}[\tau_{1}]$. Let us first
recall that $\tau_{1}=\sigma_{1}+\Theta_{1}$ is equal to the sum of the first
initiation time and the hitting time to $M$. As $\nu\to 0$, the distribution
of the first arrival $\sigma_{1}$, given in general by $1-e^{-a_{t}}$,
converges to a reverse Gumbel with parameters
$\frac{1}{\delta}\log\frac{\delta}{\nu(1-q)}$ and $-\frac{1}{\delta}$. This
implies in particular that

$\mathrm{E}[\sigma_{1}]=\frac{1}{\delta}\left(\log\frac{\delta}{\nu(1-q)}-\gamma_{E}\right)=\frac{1}{\delta}\log\frac{\delta}{\nu}+C_{1}$

where $C_{1}=-\frac{\log(1-q)+\gamma_{E}}{\delta}$. Moreover, the hitting
times $\Theta_{i}$ follow the Gumbel distribution $G(t)$ \- see equation 2 \-
and hence

$\mathrm{E}[\Theta_{i}]=\frac{1}{\lambda}\log M-C_{2}$

for every $i$, where $C_{2}=-\frac{\log(1-q)+\gamma_{E}}{\lambda}$. Joining
the last two results we get

$\mathrm{E}[\tau_{1}]=\mathrm{E}[\sigma_{1}+\Theta_{1}]=\frac{1}{\lambda}\log
M+\frac{1}{\delta}\log\frac{\delta}{\nu}+\tilde{C}$ (20)

where $\tilde{C}=C_{1}-C_{2}$. By comparing equation 20 with the expression
for $\mathrm{E}[\tau]$, we notice indeed the same $M$ and $\nu$ dependence,
but the constants $C$ and $\tilde{C}$ have different analytical forms.

####  4.3.4 Numerical computation

Finally, all the plots and computations reported in this paper have been
performed on Matlab R2018b. The lines of code below provide an efficient way
(in the example for the exponential case) to calculate the relapse time
distribution given by equation 4 for a vector of times tspan.

⬇

n = @(t)(exp(delta*t));

G = @(t)(exp(-(1-q)*M*exp(-lambda*t)));

F =
@(t)(1-exp(-nu*(1-q)*integral(@(s)(n(s).*G(t-s)),0,t,’ArrayValued’,true)));

x␣=␣arrayfun(@(t)F(t),tspan);

## Acknowledgments

We thank Ivana Bozic, David Cheek, Jasmine Foo, Kevin Leder, Michael Nicholson
and Johannes Reiter for helpful discussions.

## References

  * [1] Sahai E.  Illuminating the metastatic process.  Nature Reviews Cancer. 2007;7(10):737–749.  doi:10.1038/nrc2229.
  * [2] Naxerova K, Brachtel E, Salk JJ, Seese AM, Power K, Abbasi B, et al.  Hypermutable DNA chronicles the evolution of human colon cancer.  Proceedings of the National Academy of Sciences. 2014;111(18):E1889–E1898.  doi:10.1073/pnas.1400179111.
  * [3] Harper KL, Sosa MS, Entenberg D, Hosseini H, Cheung JF, Nobre R, et al.  Mechanism of early dissemination and metastasis in Her2$+$ mammary cancer.  Nature. 2016;540(7634):588–592.  doi:10.1038/nature20609.
  * [4] Reiter JG, Makohon-Moore AP, Gerold JM, Heyde A, Attiyeh MA, Kohutek ZA, et al.  Minimal functional driver gene heterogeneity among untreated metastases.  Science. 2018;361(6406):1033–1037.  doi:10.1126/science.aat7171.
  * [5] Michor F, Nowak MA, Iwasa Y.  Stochastic dynamics of metastasis formation.  Journal of Theoretical Biology. 2006;240(4):521–530.  doi:10.1016/j.jtbi.2005.10.021.
  * [6] Haeno H, Michor F.  The evolution of tumor metastases during clonal expansion.  Journal of Theoretical Biology. 2010;263(1):30–44.  doi:10.1016/j.jtbi.2009.11.005.
  * [7] Chaffer CL, Weinberg RA.  A Perspective on Cancer Cell Metastasis.  Science. 2011;331(6024):1559–1564.  doi:10.1126/science.1203543.
  * [8] Tsikitis VL, Larson DW, Huebner M, Lohse CM, Thompson PA.  Predictors of recurrence free survival for patients with stage II and III colon cancer.  BMC Cancer. 2014;14(1).  doi:10.1186/1471-2407-14-336.
  * [9] Luria SE, Delbrück M.  Mutations of bacteria from virus sensitivity to virus resistance.  Genetics. 1943;48(6):491–511.
  * [10] Iwasa Y, Nowak MA, Michor F.  Evolution of Resistance During Clonal Expansion.  Genetics. 2006;172(4):2557–2566.
  * [11] Komarova N.  Stochastic modeling of drug resistance in cancer.  Journal of Theoretical Biology. 2006;239:351–366.
  * [12] Foo J, Leder K.  Dynamics of cancer recurrence.  The Annals of Applied Probability. 2013;23(4):1437–1468.  doi:10.1214/12-AAP876.
  * [13] Bozic I, Reiter JG, Allen B, Antal T, Chatterjee K, Shah P, et al.  Evolutionary dynamics of cancer in response to targeted combination therapy.  eLife. 2013;2.
  * [14] Durrett R, Moseley S.  Evolution of resistance and progression to disease during clonal expansion of cancer.  Theoretical Population Biology. 2010;77:42–48.
  * [15] Durrett R, Foo J, Leder K, Mayberry J, Michor F.  Evolutionary dynamics of tumor progression with random fitness values.  Theoretical Population Biology. 2010;78(1):54–66.  doi:10.1016/j.tpb.2010.05.001.
  * [16] Nicholson MD, Antal T.  Universal Asymptotic Clone Size Distribution for General Population Growth.  Bulletin of Mathematical Biology. 2016;78(11):2243–2276.  doi:10.1007/s11538-016-0221-x.
  * [17] Dingli D, Michor F, Antal T, Pacheco JM.  The emergence of tumor metastases.  Cancer Biology & Therapy. 2007;6(3):383–390.  doi:10.4161/cbt.6.3.3720.
  * [18] Armitage P, Doll R. The Age Distribution of Cancer and a Multi-stage Theory of Carcinogenesis; 1954.
  * [19] Hanin LG, Tsodikov AD, Yakovlev AY.  Optimal schedules of cancer surveillance and tumor size at detection.  Mathematical and Computer Modelling. 2001;33(12-13):1419–1430.  doi:10.1016/s0895-7177(01)80023-6.
  * [20] Hanin L, Pavlova L.  Optimal screening schedules for prevention of metastatic cancer.  Statistics in Medicine. 2012;32(2):206–219.  doi:10.1002/sim.5474.
  * [21] Tsodikov AD, Ibrahim JG, Yakovlev AY.  Estimating Cure Rates From Survival Data.  Journal of the American Statistical Association. 2003;98(464):1063–1078.  doi:10.1198/01622145030000001007.
  * [22] Yakovlev AY.  Threshold models of tumor recurrence.  Mathematical and Computer Modelling. 1996;23(6):153–164.  doi:10.1016/0895-7177(96)00024-6.
  * [23] Yakovlev AY, Tsodikov AD, Asselain B.  Stochastic Models of Tumor Latency and Their Biostatistical Applications.  World Scientific; 1996.
  * [24] Klein CA.  Parallel progression of primary tumours and metastases.  Nature Reviews Cancer. 2009;9:302.
  * [25] Lea DE, Coulson CA.  The distribution of the numbers of mutants in bacterial populations.  Journal of Genetics. 1949;49(3):264–285.  doi:10.1007/BF02986080.
  * [26] Keller P, Antal T.  Mutant number distribution in an exponentially growing population.  Journal of Statistical Mechanics: Theory and Experiment. 2015;2015(1).
  * [27] Kendall DG.  Birth-and-death processes, and the theory of carcinogenesis.  Biometrika. 1960;47:13–21.
  * [28] Kessler DA, Levine H.  Scaling Solution in the Large Population Limit of the General Asymmetric Stochastic Luria–Delbrück Evolution Process.  Journal of Statistical Physics. 2014;158(4):783–805.  doi:10.1007/s10955-014-1143-3.
  * [29] Cheek D, Antal T.  Mutation frequencies in a birth–death branching process.  The Annals of Applied Probability. 2018;28(6):3922–3947.  doi:10.1214/18-aap1413.
  * [30] Tubiana M.  The growth and progression of human tumors: Implications for management strategy.  Radiotherapy and Oncology. 1986;6(3):167–184.  doi:10.1016/s0167-8140(86)80151-7.
  * [31] Athreya KB, Ney PE.  Branching Processes.  Dover Publications; 2004.
  * [32] Collins VP, Loeffler RK, Tivery H.  Observations on growth rates of human tumors.  The American journal of roentgenology, radium therapy, and nuclear medicine. 1956;76:988–1000.
  * [33] Preziosi L.  Cancer Modelling and Simulation.  Chapman & Hall/CRC Mathematical and Computational Biology. CRC Press, Taylor & Francis Group; 2003.
  * [34] Bolognese A, Izzo L.  Surgery in Multimodal Management of Solid Tumors.  Springer Milan; 2009.
  * [35] Peng Y, Taylor JMG.  Cure Models.  In: Handbook of Survival Analysis. Chapman & Hall; 2014. p. 113–134.
  * [36] Allison PD.  Survival analysis using SAS: a practical guide.  2nd ed. SAS Publishing; 2010.
  * [37] Meeker WQ, Escobar LA.  Statistical Methods for Reliability Data.  John Wiley & Sons Inc.; 1998.
  * [38] Singh R, Mukhopadhyay K.  Survival analysis in clinical trials: Basics and must know areas.  Perspectives in Clinical Research. 2011;2(4):145.  doi:10.4103/2229-3485.86872.
  * [39] Hagar YC, Harvey DJ, Beckett LA.  A multivariate cure model for left-censored and right-censored data with application to colorectal cancer screening patterns.  Statistics in medicine. 2016;35:3347–3367.
  * [40] Adam R, de Gramont A, Figueras J, Kokudo N, Kunstlinger F, Loyer E, et al.  Managing synchronous liver metastases from colorectal cancer: A multidisciplinary international consensus.  Cancer Treatment Reviews. 2015;41(9):729–741.  doi:10.1016/j.ctrv.2015.06.006.
  * [41] Schwartz M.  A biomathematical approach to clinical tumor growth.  Cancer. 1961;14:1272–1294.
  * [42] Spratt JS, Spratt TL.  Rates of Growth of Pulmonary Metastases and Host Survival.  Annals of Surgery. 1964;159(2):161–171.
  * [43] Jones S, Chen Wd, Parmigiani G, Diehl F, Beerenwinkel N, Antal T, et al.  Comparative lesion sequencing provides insights into tumor evolution.  Proceedings of the National Academy of Sciences. 2008;105(11):4283–4288.  doi:10.1073/pnas.0712345105.
  * [44] Haustermans K, Fowler J, Geboes K, Christiaens MR, Lerut A, van der Schueren E.  Relationship between potential doubling time (Tpot), labeling index and duration of DNA synthesis in 60 esophageal and 35 breast tumors: is it worthwhile to measure Tpot?  Radiotherapy and Oncology. 1998;46(2):157–167.  doi:10.1016/s0167-8140(97)00164-3.
  * [45] Denekamp J.  New Approaches to the Measurement of Proliferation Rates.  In: Maragoudakis ME, Gullino P, Lelkes PI, editors. Angiogenesis in Health and Disease. Boston, MA: Springer US; 1992. p. 333–337.
  * [46] Bertuzzi A, Gandolfi A, Sinisgalli C, Starace G, Ubezio P.  Cell loss and the concept of potential doubling time.  Cytometry. 1997;29(1):34–40.  doi:10.1002/(sici)1097-0320(19970901)29:1¡34::aid-cyto3¿3.0.co;2-d.
  * [47] von Fournier D, Weber E, Hoeffken W, Bauer M, Kubli F, Barth V.  Growth rate of 147 mammary carcinomas.  Cancer. 1980;45:2198–2207.
  * [48] Kuroishi T, Tominaga S, Morimoto T, Tashiro H, Itoh S, Watanabe H, et al.  Tumor Growth Rate and Prognosis of Breast Cancer Mainly Detected by Mass Screening.  Japanese Journal of Cancer Research. 1990;81(5):454–462.  doi:10.1111/j.1349-7006.1990.tb02591.x.
  * [49] Peer PGM, Van Dijck JAAM, Verbeek ALM, Hendriks JHCL, Holland R.  Age-dependent growth rate of primary breast cancer.  Cancer. 1993;71(11):3547–3551.  doi:10.1002/1097-0142(19930601)71:11¡3547::aid-cncr2820711114¿3.0.co;2-c.
  * [50] Ryu EB, Chang JM, Seo M, Kim SA, Lim JH, Moon WK.  Tumour volume doubling time of molecular breast cancer subtypes assessed by serial breast ultrasound.  European Radiology. 2014;24(9):2227–2235.  doi:10.1007/s00330-014-3256-0.
  * [51] Förnvik D, Lång K, Andersson I, Dustler M, Borgquist S, Timberg P.  Estimates of breast cancer growth rate from mammograms and its relation to tumour characteristics.  Radiation Protection Dosimetry. 2015;169(1-4):151–157.  doi:10.1093/rpd/ncv417.
  * [52] Zhang S, Ding Y, Zhou Q, Wang C, Wu P, Dong J.  Correlation Factors Analysis of Breast Cancer Tumor Volume Doubling Time Measured by 3D-Ultrasound.  Medical Science Monitor. 2017;23:3147–3153.  doi:10.12659/msm.901566.
  * [53] Kusama S, Spratt Jr JS, Donegan WL, Watson FR, Cunningham C.  The gross rates of growth of human mammary carcinoma.  Cancer. 1972;30(2):594–599.  doi:10.1002/1097-0142(197208)30:2¡594::aid-cncr2820300241¿3.0.co;2-2.
  * [54] Friberg S, Mattson S.  On the growth rates of human malignant tumors: implications for medical decision making.  Journal of surgical oncology. 1997;65:284–297.
  * [55] Awwad H.  Radiation Oncology: Radiobiological and Physiological Perspectives.  Springer Netherlands; 2013.
  * [56] Zabicki K, Colbert JA, Dominguez FJ, Gadd MA, Hughes KS, Jones JL, et al.  Breast Cancer Diagnosis in Women $\leq$ 40 versus 50 to 60 Years: Increasing Size and Stage Disparity Compared With Older Women Over Time.  Annals of Surgical Oncology. 2006;13(8):1072–1077.  doi:10.1245/aso.2006.03.055.
  * [57] Lee SH, Kim YS, Han W, Ryu HS, Chang JM, Cho N, et al.  Tumor growth rate of invasive breast cancers during wait times for surgery assessed by ultrasonography.  Medicine. 2016;95(37):e4874.  doi:10.1097/md.0000000000004874.
  * [58] de l’Aulnoit AH, Rogoz B, Pinçon C, de l’Aulnoit DH.  Metastasis-free interval in breast cancer patients: thirty-year trends and time dependency of prognostic factors. A retrospective analysis based on a single institution experience.  The Breast. 2018;37:80–88.  doi:10.1016/j.breast.2017.10.008.
  * [59] Bolin S, Nilsson E, Sjödahl R. Carcinoma of the colon and rectum–growth rate.; 1983.
  * [60] Tada M, Misaki F, Kawai K.  Growth rates of colorectal carcinoma and adenoma by roentgenologic follow-up observations.  Gastroenterologia Japonica. 1984;19:550–555.
  * [61] Choi SJ, Kim HS, Ahn SJ, Jeong YM, Choi HY.  Evaluation of the growth pattern of carcinoma of colon and rectum by MDCT.  Acta Radiologica. 2013;54(5):487–492.  doi:10.1177/0284185113475923.
  * [62] Finlay IG, Meek D, Bruntont F, McArdle CS.  Growth rate of hepatic metastases in colorectal carcinoma.  British Journal of Surgery. 1988;75(7):641–644.  doi:10.1002/bjs.1800750707.
  * [63] Tanaka K, Shimada H, Miura M, Fujii Y, Yamaguchi S, Endo I, et al.  Metastatic Tumor Doubling Time: Most Important Prehepatectomy Predictor of Survival and Nonrecurrence of Hepatic Colorectal Cancer Metastasis.  World Journal of Surgery. 2004;28(3):263–270.  doi:10.1007/s00268-003-7088-3.
  * [64] Tomimaru Y, Noura S, Ohue M, Okami J, Oda K, Higashiyama M, et al.  Metastatic Tumor Doubling Time Is an Independent Predictor of Intrapulmonary Recurrence after Pulmonary Resection of Solitary Pulmonary Metastasis from Colorectal Cancer.  Digestive Surgery. 2008;25(3):220–225.  doi:10.1159/000140693.
  * [65] Wilson MS, West CM, Wilson GD, Roberts SA, James RD, Schofield PF.  Intra-tumoral heterogeneity of tumour potential doubling times (Tpot) in colorectal cancer.  British journal of cancer. 1993;68:501–506.
  * [66] Kornprat P, Pollheimer MJ, Lindtner RA, Schlemmer A, Rehak P, Langner C.  Value of Tumor Size as a Prognostic Variable in Colorectal Cancer.  American Journal of Clinical Oncology. 2011;34(1):43–49.  doi:10.1097/coc.0b013e3181cae8dd.
  * [67] Ding Z, Wang Z, Huang S, Zhong S, Lin J.  Comparison of laparoscopic vs. open surgery for rectal cancer.  Molecular and Clinical Oncology. 2017;6(2):170–176.  doi:10.3892/mco.2016.1112.
  * [68] Waaijer A, Terhaard CHJ, Dehnad H, Hordijk GJ, van Leeuwen MS, Raaymakers CPJ, et al.  Waiting times for radiotherapy: consequences of volume increase for the TCP in oropharyngeal carcinoma.  Radiotherapy and Oncology. 2003;66(3):271–276.  doi:10.1016/s0167-8140(03)00036-7.
  * [69] Jensen AR, Nellemann HM, Overgaard J.  Tumor progression in waiting time for radiotherapy in head and neck cancer.  Radiotherapy and Oncology. 2007;84(1):5–10.  doi:10.1016/j.radonc.2007.04.001.
  * [70] Galante E, Gallus G, Chiesa F, Bono A, Bettoni I, Molinari R.  Growth rate of head and neck tumors.  European Journal of Cancer and Clinical Oncology. 1982;18(8):707–712.  doi:10.1016/0277-5379(82)90067-0.
  * [71] Umino S, Hayashi S, Ono S.  Doubling time of pulmonary metastases of adenoid cystic carcinoma.  International Journal of Oral and Maxillofacial Surgery. 1997;26:48.  doi:10.1016/s0901-5027(97)80987-3.
  * [72] Zackrisson B, Gustafsson H, Stenling R, Flygare P, Wilson GD.  Predictive value of potential doubling time in head and neck cancer patients treated by conventional radiotherapy.  International Journal of Radiation Oncology$\ast$Biology$\ast$Physics. 1997;38(4):677–683.  doi:10.1016/s0360-3016(97)00066-7.
  * [73] Muto M, Nakane M, Katada C, Sano Y, Ohtsu A, Esumi H, et al.  Squamous cell carcinoma in situ at oropharyngeal and hypopharyngeal mucosal sites.  Cancer. 2004;101(6):1375–1381.  doi:10.1002/cncr.20482.
  * [74] Markou K, Goudakos J, Triaridis S, Konstantinidis J, Vital V, Nikolaou A.  The role of tumor size and patient’s age as prognostic factors in laryngeal cancer.  Hippokratia. 2011;15(21607041):75–80.
  * [75] Kerr KM, Lamb D.  Actual growth rate and tumour cell proliferation in human pulmonary neoplasms.  British Journal Of Cancer. 1984;50:343.
  * [76] Arai T, Kuroishi T, Saito Y, Kurita Y, Naruke T, Kaneko M.  Tumor Doubling Time and Prognosis in Lung Cancer Patients: Evaluation from Chest Films and Clinical Follow-up Study.  Japanese Journal of Clinical Oncology. 1994;doi:10.1093/oxfordjournals.jjco.a039706.
  * [77] Detterbeck FC, Gibson CJ.  Turning Gray: The Natural History of Lung Cancer Over Time.  Journal of Thoracic Oncology. 2008;3(7):781–792.  doi:10.1097/jto.0b013e31817c9230.
  * [78] Henschke CI, Yankelevitz DF, Yip R, Reeves AP, Farooqi A, Xu D, et al.  Lung Cancers Diagnosed at Annual CT Screening: Volume Doubling Times.  Radiology. 2012;263(2):578–583.  doi:10.1148/radiol.12102489.
  * [79] Yoo H, Nam BH, Yang HS, Shin SH, Lee JS, Lee SH.  Growth rates of metastatic brain tumors in nonsmall cell lung cancer.  Cancer. 2008;113(5):1043–1047.  doi:10.1002/cncr.23676.
  * [80] Fowler JF.  Biological Factors Influencing Optimum Fractionation in Radiation Therapy.  Acta Oncologica. 2001;40(6):712–717.  doi:10.1080/02841860152619124.
  * [81] Bando T.  A new method of segmental resection for primary lung cancer: intermediate results.  European Journal of Cardio-Thoracic Surgery. 2002;21(5):894–899.  doi:10.1016/s1010-7940(02)00122-7.
  * [82] Strand TE.  Survival after resection for primary lung cancer: a population based study of 3211 resected patients.  Thorax. 2006;61(8):710–715.  doi:10.1136/thx.2005.056481.
  * [83] D’Amico AV, Hanks GE.  Linear regressive analysis using prostate-specific antigen doubling time for predicting tumor biology and clinical outcome in prostate cancer.  Cancer. 1993;72:2638–2643.
  * [84] Werahera PN, Glode LM, Rosa FGL, Lucia MS, Crawford ED, Easterday K, et al.  Proliferative Tumor Doubling Times of Prostatic Carcinoma.  Prostate Cancer. 2011;2011:1–7.  doi:10.1155/2011/301850.
  * [85] Zharinov GM, Bogomolov OA, Neklasova NN, Anisimov VN.  Pretreatment prostate specific antigen doubling time as prognostic factor in prostate cancer patients.  Oncoscience. 2017;4:7–13.  doi:10.18632/oncoscience.337.
  * [86] Berges RR, Vukanovic J, Epstein JI, CarMichel M, Cisek L, Johnson DE, et al.  Implication of cell kinetic changes during the progression of human prostatic cancer.  Clinical cancer research : an official journal of the American Association for Cancer Research. 1995;1:473–480.
  * [87] Haustermans KMG, Hofland I, Poppel HV, Oyen R, de Voorde WV, Begg AC, et al.  Cell kinetic measurements in prostate cancer.  International Journal of Radiation Oncology$\ast$Biology$\ast$Physics. 1997;37(5):1067–1070.  doi:10.1016/s0360-3016(96)00579-2.
  * [88] Renshaw AA, Richie JP, Loughlin KR, Jiroutek M, Chung A, D’Amico AV.  Maximum diameter of prostatic carcinoma is a simple, inexpensive, and independent predictor of prostate-specific antigen failure in radical prostatectomy specimens. Validation in a cohort of 434 patients.  American journal of clinical pathology. 1999;111:641–644.
  * [89] Johnson SB, Hamstra DA, Jackson WC, Zhou J, Foster B, Foster C, et al.  Larger Maximum Tumor Diameter at Radical Prostatectomy Is Associated With Increased Biochemical Failure, Metastasis, and Death From Prostate Cancer After Salvage Radiation for Prostate Cancer.  International Journal of Radiation Oncology$\ast$Biology$\ast$Physics. 2013;87(2):275–281.  doi:10.1016/j.ijrobp.2013.05.043.
  * [90] Serres S, Soto MS, Hamilton A, McAteer MA, Carbonell WS, Robson MD, et al.  Molecular MRI enables early and sensitive detection of brain metastases.  Proceedings of the National Academy of Sciences. 2012;109(17):6674–6679.  doi:10.1073/pnas.1117412109.
  * [91] Fujiwara S, Yao K, Nagahama T, Uchita K, Kanemitsu T, Tsurumi K, et al.  Can we accurately diagnose minute gastric cancers ($\leq 5$mm)? Chromoendoscopy (CE) vs magnifying endoscopy with narrow band imaging (M-NBI).  Gastric Cancer. 2015;18(3):590–596.
  * [92] Wang L.  Early Diagnosis of Breast Cancer.  Sensors. 2017;17(7):1572.  doi:10.3390/s17071572.
  * [93] Chignola R, Foroni RI.  Estimating the Growth Kinetics of Experimental Tumors From as Few as Two Determinations of Tumor Size: Implications for Clinical Oncology.  IEEE Transactions on Biomedical Engineering. 2005;52(5):808–815.  doi:10.1109/tbme.2005.845219.
  * [94] Fillon M.  Better Guidelines Needed for Cancer Survivorship Management.  CA: A Cancer Journal for Clinicians. 2018;68(6):392–393.  doi:10.3322/caac.21435.
  * [95] Lee SP, Sun JR, Qian H, McBride WH, Withers HR.  Characterization of Metastatic Tumor Formation by the Colony Size Distribution.  arXiv pre-print. 2006;.
  * [96] Schmid HP, McNeal JE, Stamey TA.  Observations on the doubling time of prostate cancer. The use of serial prostate-specific antigen in patients with untreated disease as a measure of increasing cancer volume.  Cancer. 1993;71(6):2031–2040.  doi:10.1002/1097-0142(19930315)71:6¡2031::aid-cncr2820710618¿3.0.co;2-q.
  * [97] Andre F, Slimane K, Bachelot T, Dunant A, Namer M, Barrelier A, et al.  Breast Cancer With Synchronous Metastases: Trends in Survival During a 14-Year Period.  Journal of Clinical Oncology. 2004;22(16):3302–3308.  doi:10.1200/jco.2004.08.095.
  * [98] Boutros C, Mazouni C, Lerebours F, Stevens D, Lei X, Gonzalez-Angulo AM, et al.  A preoperative nomogram to predict the risk of synchronous distant metastases at diagnosis of primary breast cancer.  British Journal of Cancer. 2015;112(6):992–997.  doi:10.1038/bjc.2015.34.
  * [99] Yilmaz U, Marks LB.  Estimating changes in the rate of synchronous and metachronous metastases over time: Analysis of SEER data.  Advances in Radiation Oncology. 2018;3(1):70–75.  doi:10.1016/j.adro.2017.09.007.
  * [100] Kim H, Choi DH, Park W, Huh SJ, Nam SJ, Lee JE, et al.  Prognostic factors for survivals from first relapse in breast cancer patients: analysis of deceased patients.  Radiation Oncology Journal. 2013;31(4):222.  doi:10.3857/roj.2013.31.4.222.
  * [101] Fitzpatrick DJ, Lai CS, Parkyn RF, Walters D, Humeniuk V, Walsh DCA.  Time to Breast Cancer Relapse Predicted By Primary Tumour Characteristics, Not Lymph Node Involvement.  World Journal of Surgery. 2013;38(7):1668–1675.  doi:10.1007/s00268-013-2397-7.
  * [102] Nowikiewicz T, Wiśniewska M, Wiśniewski M, Biedka M, Głowacka I, Kozak D, et al.  Overall survival and disease-free survival in breast cancer patients treated at the Oncology Centre in Bydgoszcz – analysis of more than six years of follow-up.  Współczesna Onkologia. 2015;4:284–289.  doi:10.5114/wo.2015.54387.
  * [103] Kemeny MM, Adak S, Gray B, Macdonald JS, Smith T, Lipsitz S, et al.  Combined-Modality Treatment for Resectable Metastatic Colorectal Carcinoma to the Liver: Surgical Resection of Hepatic Metastases in Combination With Continuous Infusion of Chemotherapy - An Intergroup Study.  Journal of Clinical Oncology. 2002;20(6):1499–1505.  doi:10.1200/JCO.2002.20.6.1499.
  * [104] Park JH, Kim TY, Lee KH, Han SW, Oh DY, Im SA, et al.  The beneficial effect of palliative resection in metastatic colorectal cancer.  British Journal Of Cancer. 2013;108:1425.
  * [105] Lykoudis PM, O’Reilly D, Nastos K, Fusai G.  Systematic review of surgical management of synchronous colorectal liver metastases.  Br J Surg. 2014;101(6):605–612.
  * [106] Elferink MAG, de Jong KP, Klaase JM, Siemerink EJ, de Wilt JHW.  Metachronous metastases from colorectal cancer: a population-based study in North-East Netherlands.  International Journal of Colorectal Disease. 2015;30(2):205–212.  doi:10.1007/s00384-014-2085-6.
  * [107] Holch JW, Demmer M, Lamersdorf C, Michl M, Schulz C, von Einem JC, et al.  Pattern and Dynamics of Distant Metastases in Metastatic Colorectal Cancer.  Visceral Medicine. 2017;33(1):70–75.  doi:10.1159/000454687.
  * [108] Hohenberger P, Schlag PM, Gerneth T, Herfarth C.  Pre- and postoperative carcinoembryonic antigen determinations in hepatic resection for colorectal metastases. Predictive value and implications for adjuvant treatment based on multivariate analysis.  Annals of Surgery. 1994;219:135–143.
  * [109] Nordlinger B, Van Cutsem E, Gruenberger T, Glimelius B, Poston G, Rougier P, et al.  Combination of surgery and chemotherapy and the role of targeted agents in the treatment of patients with colorectal liver metastases: recommendations from an expert panel.  Annals of Oncology. 2009;20(6):985–992.  doi:10.1093/annonc/mdn735.
  * [110] Sturesson C, Valdimarsson VT, Blomstrand E, Eriksson S, Nilsson JH, Syk I, et al.  Liver-first strategy for synchronous colorectal liver metastases - an intention-to-treat analysis.  HPB. 2017;19(1):52–58.  doi:https://doi.org/10.1016/j.hpb.2016.10.005.
  * [111] Ferlito A, Shaha AR, Silver CE, Rinaldo A, Mondin V.  Incidence and Sites of Distant Metastases from Head and Neck Cancer.  ORL. 2001;63(4):202–207.  doi:10.1159/000055740.
  * [112] Jain KS, Sikora AG, Baxi SS, Morris LGT.  Synchronous cancers in patients with head and neck cancer.  Cancer. 2013;119(10):1832–1837.  doi:10.1002/cncr.27988.
  * [113] Liu SA, Wong YK, Lin JC, Poon CK, Tung KC, Tsai WC.  Impact of recurrence interval on survival of oral cavity squamous cell carcinoma patients after local relapse.  Otolaryngology-Head and Neck Surgery. 2007;136(1):112–118.  doi:10.1016/j.otohns.2006.07.002.
  * [114] Ebrahimi A, Clark JR, Ahmadi N, Palme CE, Morgan GJ, Veness MJ.  Prognostic significance of disease-free interval in head and neck cutaneous squamous cell carcinoma with nodal metastases.  Head & Neck. 2012;35(8):1138–1143.  doi:10.1002/hed.23096.
  * [115] Wiegand S, Zimmermann A, Wilhelm T, Werner JA.  Survival After Distant Metastasis in Head and Neck Cancer.  Anticancer research. 2015;35:5499–5502.
  * [116] Tönnies M, Pfannschmidt J, Bauer TT, Kollmeier J, Tönnies S, Kaiser D.  Metastasectomy for Synchronous Solitary Non-Small Cell Lung Cancer Metastases.  The Annals of Thoracic Surgery. 2014;98(1):249–256.  doi:10.1016/j.athoracsur.2014.03.028.
  * [117] al Kattan K, Sepsas E, Fountain SW, Townsend ER.  Disease recurrence after resection for stage I lung cancer.  European journal of cardio-thoracic surgery : official journal of the European Association for Cardio-thoracic Surgery. 1997;12:380–384.
  * [118] Hung JJ, Jeng WJ, Hsu WH, Wu KJ, Chou TY, Hsieh CC, et al.  Prognostic factors of postrecurrence survival in completely resected stage I non-small cell lung cancer with distant metastasis.  Thorax. 2010;65(3):241–245.  doi:10.1136/thx.2008.110825.
  * [119] Farsi AA, Swaminath A, Ellis P.  Patterns of Relapse in Small Cell Lung Cancer (SCLC): A Retrospective Analysis of Outcomes from a Single Canadian Center.  Journal of Thoracic Oncology. 2017;12(1):S727–S728.  doi:10.1016/j.jtho.2016.11.962.
  * [120] Koo KC, Yoo H, Kim KH, Park SU, Han KS, Rha KH, et al.  Prognostic Impact of Synchronous Second Primary Malignancies on the Overall Survival of Patients with Metastatic Prostate Cancer.  Journal of Urology. 2015;193(4):1239–1244.  doi:10.1016/j.juro.2014.10.088.
  * [121] Fontenot PA Jr, Nehra A, Parker W, Wyre H, Mirza M, Duchene DA, et al.  Metastatic prostate cancer in the modern era of PSA screening.  International braz j urol. 2017;43(3):416–421.  doi:10.1590/s1677-5538.ibju.2016.0340.
  * [122] Almeida PL, Pereira BJ.  Local treatment of metastatic prostate cancer: what is the evidence so far?  Prostate Cancer. 2018;2018:1–7.  doi:10.1155/2018/2654572.
  * [123] Boorjian SA, Thompson RH, Tollefson MK, Rangel LJ, Bergstralh EJ, Blute ML, et al.  Long-Term Risk of Clinical Progression After Biochemical Recurrence Following Radical Prostatectomy: The Impact of Time from Surgery to Recurrence.  European Urology. 2011;59(6):893–899.  doi:10.1016/j.eururo.2011.02.026.
  * [124] Toussi A, Stewart-Merrill SB, Boorjian SA, Psutka SP, Thompson RH, Frank I, et al.  Standardizing the Definition of Biochemical Recurrence after Radical Prostatectomy—What Prostate Specific Antigen Cut Point Best Predicts a Durable Increase and Subsequent Systemic Progression?  Journal of Urology. 2016;195(6):1754–1759.  doi:10.1016/j.juro.2015.12.075.
  * [125] Obenauf AC, Massagué J.  Surviving at a Distance: Organ-Specific Metastasis.  Trends in Cancer. 2015;1(1):76–91.  doi:10.1016/j.trecan.2015.07.009.
  * [126] Durrett R.  Branching process models of cancer. vol. 1.1 of Stochastics in biological systems.  1st ed. Springer International Publishing; 2015.
  * [127] Haeno H, Gonen M, Davis MB, Herman JM, Iacobuzio-Donahue CA, Michor F.  Computational Modeling of Pancreatic Cancer Reveals Kinetics of Metastasis Suggesting Optimum Treatment Strategies.  Cell. 2012;148(1-2):362–375.  doi:10.1016/j.cell.2011.11.060.
  * [128] Yachida S, Jones S, Bozic I, Antal T, Leary R, Fu B, et al.  Distant metastasis occurs late during the genetic evolution of pancreatic cancer.  Nature. 2010;467:1114.
  * [129] Nicholson MD, Antal T.  Competing evolutionary paths in growing populations with applications to multidrug resistance.  PLOS Computational Biology. 2019;15(4):e1006866.  doi:10.1371/journal.pcbi.1006866.
  * [130] Gundem G, , Loo PV, Kremeyer B, Alexandrov LB, Tubio JMC, et al.  The evolutionary history of lethal metastatic prostate cancer.  Nature. 2015;520(7547):353–357.  doi:10.1038/nature14347.
  * [131] Bozic I, Antal T, Ohtsuki H, Carter H, Kim D, Chen S, et al.  Accumulation of driver and passenger mutations during tumor progression.  Proceedings of the National Academy of Sciences. 2010;107(43):18545–18550.  doi:10.1073/pnas.1010978107.
  * [132] Williams T.  The Basic Birth-Death Model for Microbial Infections.  Journal of the Royal Statistical Society Series B (Methodological). 1965;27(2):338–360.
  * [133] Waugh WAO.  Uses of the sojourn time series for Markovian birth process.  In: Proceedings of the Sixth Berkeley Symposium on Mathematical Statistics and Probability, Volume 3: Probability Theory. Berkeley, CA: University of California Press; 1972. p. 501–514.