Cancer recurrence times from a branching process model

TL;DR

This paper models cancer metastasis formation using a branching process, showing that early surgical resection is crucial to prevent recurrence and that delays significantly increase recurrence risk.

Contribution

It introduces a simple, parameterized model of metastasis formation that links primary tumor size, metastasis initiation, and recurrence timing across multiple cancer types.

Findings

Early resections can prevent recurrence in some cases.

Small delays in surgery increase recurrence probability.

Model predictions align with clinical literature data.

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.

Tables

Table 1. Table 1: Typical ranges of volume doubling times for the primary tumor ( D ​ T p ​ t 𝐷 subscript 𝑇 𝑝 𝑡 DT_{pt} ) and metastasis ( D ​ T m 𝐷 subscript 𝑇 𝑚 DT_{m} ), tumor potential doubling time ( T p ​ o ​ t subscript 𝑇 𝑝 𝑜 𝑡 T_{pot} ) and tumor diameter at resection ( d p ​ t subscript 𝑑 𝑝 𝑡 d_{pt} ) for breast, colorectal, headneck, lung and prostate cancer.

Cancer type	Parameter	Typical range	Estimate	References
Breast	$D{T}_{pt}$ (days)	$103-353$	$210$	[47, 48, 49, 50, 51, 52]
	$D{T}_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	$D{T}_{pt}$ (days)	$130-438$	$175$	[59, 60, 61]
	$D{T}_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	$D{T}_{pt}$ (days)	$15-256$	$84$	[68, 69]
	$D{T}_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	$D{T}_{pt}$ (days)	$22-269$	$168$	[75, 76, 54, 77, 78]
	$D{T}_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	$D{T}_{pt}$ (days)	$36-1080$	$392$	[83, 84, 85]
	$D{T}_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]

Table 2. 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 𝑞 q , the primary tumor size at resection N 𝑁 N and the minimal detectable size M 𝑀 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$

Table 3. Table 3: Resection sizes of the primary tumor which yield a 1 % percent 1 1\% and 99 % percent 99 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 𝑁 N , and tumor diameter, d 𝑑 d .

		Breast	Colorectal	Headneck	Lung	Prostate
$P\left(S_T\ge 1\right)>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$
$P\left(S_T\ge 1\right)>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$

Table 4. Table 4: Typical ranges of P ​ ( S T ≥ 1 ) P subscript 𝑆 𝑇 1 \,\mathrm{P}(S_{T}\geq 1) and E ​ [ τ − T ∣ U T ] E delimited-[] 𝜏 conditional 𝑇 subscript 𝑈 𝑇 \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\left(S_T\ge 1\right)$ ($\%$)	$5-10$	$6.13$	[97, 98, 99]
	$E\left[\tau -T\mid U_T\right]$ (days)	$590-1022$	$725$	[100, 101, 102]
Colorectal	$P\left(S_T\ge 1\right)$ ($\%$)	$15-25$	$20.17$	[103, 104, 105, 106, 107, 99]
	$E\left[\tau -T\mid U_T\right]$ (days)	$353-760$	$356$	[108, 109, 106, 110, 107]
Headneck	$P\left(S_T\ge 1\right)$ ($\%$)	$1-16.8$	$1.65$	[111, 112]
	$E\left[\tau -T\mid U_T\right]$ (days)	$219-623$	$435$	[113, 114, 115]
Lung	$P\left(S_T\ge 1\right)$ ($\%$)	$30-55.39$	$33.96$	[99, 116]
	$E\left[\tau -T\mid U_T\right]$ (days)	$210-602$	$249$	[117, 118, 119]
Prostate	$P\left(S_T\ge 1\right)$ ($\%$)	$10-34$	$13.53$	[120, 121, 122]
	$E\left[\tau -T\mid U_T\right]$ (days)	$730-1131$	$969$	[123, 124]

Equations

$$a_t=\mathrm{E}[K_t]=\nu (1-q){\int }_0^{t}n(s)ds$$

$$\mathrm{P}(K_t\ge 1)=1-e^{-a_t}$$

$$\mathrm{P}({\Theta }_i\le t\mid {\Omega }_{\infty }^{(i)})\sim G(t)\equiv e^{-(1-q)M{e}^{-\lambda t}}$$

$$b_t=\mathrm{E}[S_t]=\nu (1-q){\int }_0^{t}n(s)G(t-s)ds$$

$$\mathrm{P}(\tau \le t)=\mathrm{P}(S_t\ge 1)=1-e^{-b_t}$$

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

$$\mathrm{P}(\bar{\tau }\le t)=1-\exp \left(-\nu (1-q){\int }_{-\infty }^{t}n(t-s)e^{-(1-q)e^{-\lambda s}}ds\right)$$

$$\bar{\tau }\sim {\mathrm{Gumb}}_{\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)$$

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

$$\mathrm{P}(\tau \le t\mid K_T\ge 1)=\frac{\mathrm{P}(\tau \le t)}{\mathrm{P}(K_T\ge 1)}=\frac{1-e^{-b_t}}{1-e^{-a_T}}$$

$$\mathrm{P}(\tau \le t\mid K_T\ge 1)\to \mathrm{P}(\tau \le t)$$

$$S_t\sim \mathrm{Pois}\left(b_t\right),M_t\sim \mathrm{Pois}\left(c_t\right)$$

$$c_t=a_t-b_t$$

$$\mathrm{P}(S_T=0)=\mathrm{P}(\tau >T)=e^{-b_T}$$

$$P(K_T\ge 1\mid S_T=0)=\mathrm{P}(M_T\ge 1)=1-e^{-c_T}$$

$$\delta =\log 2/D{T}_{pt},\lambda =\log 2/D{T}_m$$

$$\alpha =1/T_{pot}$$

$$\nu =\frac{\delta e^{-{\gamma }_E}}{1-q}{e}^{-\delta \mathrm{E}[{\sigma }_1]}$$

$$(D{T}_{pt},D{T}_m,T_{pot},d_{pt},d_m,e^{\delta \mathrm{E}[{\sigma }_1]})$$

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

$$U_T:=\{K_T\ge 1,S_T=0\}$$

$$\begin{array}{rl}\mathrm{P}(U_T)& =\mathrm{P}(M_T\ge 1,S_T=0)=\mathrm{P}(M_T\ge 1)\mathrm{P}(S_T=0)\\ & =e^{-b_T}-e^{-a_T}=\mathrm{P}(K_T\ge 1)-\mathrm{P}(S_T\ge 1)\end{array}$$

$$\mathrm{P}(\tau \le t\mid U_T)=\frac{\mathrm{P}(T<\tau \le t\mid K_T\ge 1)}{\mathrm{P}(T<\tau \mid K_T\ge 1)}=\frac{e^{-b_T}-e^{-b_t}}{e^{-b_T}-e^{-a_T}}$$

$$\mathrm{P}(\tau \le t\mid S_T=0)=1-e^{-(b_t-b_T)}$$

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

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

$$\mathrm{P}(T_M\le t\mid {\Omega }_{\infty }\cap \{Z_0=1\})\sim e^{-(1-q)M{e}^{-\lambda t}}\equiv G(t)\text{ as }M\to \infty$$

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

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

$$Y\sim \mathrm{Gumb}(a,b)⟹\mathrm{E}[Y]=a+b{\gamma }_E,\mathrm{Var}(Y)=\frac{{\pi }^2}{6}{b}^2$$

Full text

1. 1 Introduction

2. 2 Results

3. 2.1 Model setup

4. 2.2 Time to reach detectable size

5. 2.3 Exponential population growth

6. 2.4 Relapse time with resection

7. 2.5 Metastasis classification

8. 3 Discussion

9. 3.1 Parameter estimation

10. 3.2 Model predictions

11. 3.3 Conclusions

12. 4 Materials and methods

13. 4.1 Single type process

14. 4.2 Scaled relapse time distribution

15. 4.3 Explicit results for exponential primary growth

16. 4.3.1 Scaled relapse time

17. 4.3.2 Small initiation limit

18. 4.3.3 Mean relapse time

19. 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

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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)

[TABLE]

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

[TABLE]

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.

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

[TABLE]

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

[TABLE]

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$.

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

[TABLE]

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

[TABLE]

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.

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

[TABLE]

where

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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.

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

[TABLE]

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

[TABLE]

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

[TABLE]

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].

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.

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.

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

[TABLE]

which has probability (see equations 10 and 9)

[TABLE]

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.

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

[TABLE]

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.

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]).

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

[TABLE]

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.

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)$

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

For both of these distributions we have

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

Hence, as $M$ tends to infinity we obtain

[TABLE]

where

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

This can be equivalently written as

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

where

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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)Mexp(-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.