Quantifying metastatic inefficiency: rare genotypes versus rare dynamics

We provide details of the calculation of various statistical quantities, starting with the probability of establishment of a micrometastasis of size N. By defining ${{\hat{P}}_{n}}(s)$ to be the Laplace transform of P_n (t), the master equations and boundary conditions take the form of the following $N+1$ algebraic equations:

$$\begin{eqnarray}&&s{{\hat{P}}_{0}}={{\hat{P}}_{1}},\end{eqnarray} \tag{ 10 }$$

$$\begin{eqnarray}&&s{{\hat{P}}_{1}}=2{{\hat{P}}_{2}}-(1+\theta ){{\hat{P}}_{1}}+1,\end{eqnarray} \tag{ 11 }$$

$$\begin{eqnarray}&&s{{\hat{P}}_{2}}=3{{\hat{P}}_{3}}-2(1+\theta ){{\hat{P}}_{2}}+\theta {{\hat{P}}_{1}},\end{eqnarray} \tag{ 12 }$$

$$\begin{eqnarray*}&&\vdots \end{eqnarray*}$$

$$\begin{eqnarray}&&s{{\hat{P}}_{n}}=(n+1){{\hat{P}}_{n+1}}-n(1+\theta ){{\hat{P}}_{n}}+(n-1)\theta {{\hat{P}}_{n-1}},\end{eqnarray} \tag{ 13 }$$

$$\begin{eqnarray*}&&\vdots \end{eqnarray*}$$

$$\begin{eqnarray}&&s{{\hat{P}}_{N-1}}=-(N-1)(1+\theta ){{\hat{P}}_{N-1}}+(N-2)\theta {{\hat{P}}_{N-2}},\end{eqnarray} \tag{ 14 }$$

$$\begin{eqnarray}&&s{{\hat{P}}_{N}}=(N-1)\theta {{\hat{P}}_{N-1}}.\end{eqnarray} \tag{ 15 }$$

Equations (11)–(14) form a closed set of $N-1$ equations and can be solved explicitly.

We now define the set of functions $\{{{H}_{n}}(s)\}$ via the recursion relation

$$\begin{eqnarray}&&{{H}_{n}}(s)=\frac{\theta }{s/n+\theta +1-{{H}_{n-1}}(s)},\end{eqnarray} \tag{ 16 }$$

with ${{H}_{0}}=0$. Then solving the algebraic equations (11)–(13) iteratively, one finds

$$\begin{eqnarray}&&\theta n{{\hat{P}}_{n}}={{H}_{1}}\cdots {{H}_{n}}+(n+1){{\hat{P}}_{n+1}}{{H}_{n}}.\end{eqnarray} \tag{ 17 }$$

Setting $n=N-2$ in (17) and eliminating ${{\hat{P}}_{N-2}}$ using (14) and then (5) we obtain a simple expression for ${{\hat{Q}}_{N}}(s)$:

$$\begin{eqnarray}&&{{\hat{Q}}_{N}}(s)={{H}_{1}}\cdots {{H}_{N-1}}.\end{eqnarray} \tag{ 18 }$$

It is subsequently convenient here to work with the set of functions $\{{{J}_{n}}(s)\}$ defined by

$$\begin{eqnarray}&&{{J}_{n}}(s)={{\theta }^{n}}/({{H}_{1}}\cdots {{H}_{n}}).\end{eqnarray} \tag{ 19 }$$

Then ${{H}_{n}}=\theta {{J}_{n-1}}/{{J}_{n}}$ and using (16) we get the recursion relation:

$$\begin{eqnarray}&&{{J}_{n}}(s)=\left( \frac{s}{n}+1+\theta \right){{J}_{n-1}}(s)-\theta {{J}_{n-2}}(s),\end{eqnarray} \tag{ 20 }$$

with ${{J}_{0}}\equiv 1$ and ${{J}_{1}}=s+1+\theta$. The function of interest is therefore

$$\begin{eqnarray}&&{{\hat{Q}}_{N}}(s)=\frac{{{\theta }^{N-1}}}{{{J}_{N-1}}(s)}.\end{eqnarray} \tag{ 21 }$$

Now, using the definition of the Laplace transform, the integrals of Q_N (t) with weights 1, t, and t² are respectively given by ${{\hat{Q}}_{N}}(0)$, $-\hat{Q}_{N}^{^{\prime} }(0)$, and $\hat{Q}_{N}^{^{\prime\prime} }(0)$. Therefore the probability of successful establishment (equation (6)) is

$$\begin{eqnarray}&&{{P}_{N}}=\mathop{\int }\limits_{0}^{\infty }{\rm d}t\ {{Q}_{N}}(t)={{\hat{Q}}_{N}}(0)=\frac{{{\theta }^{N-1}}}{{{J}_{N-1}}(0)}.\end{eqnarray} \tag{ 22 }$$

From equation (20) we have

$$\begin{eqnarray}&&{{J}_{n}}(0)=(1+\theta ){{J}_{n-1}}(0)-\theta {{J}_{n-2}}(0),\end{eqnarray} \tag{ 23 }$$

with ${{J}_{0}}(0)=1$ and ${{J}_{1}}(0)=1+\theta$. This recursion relation may be iterated to find the simple result

$$\begin{eqnarray}&&{{J}_{n}}(0)=1+\theta +{{\theta }^{2}}+...{{\theta }^{n}}\end{eqnarray} \tag{ 24 }$$

thus

$$\begin{eqnarray}&&{{J}_{N-1}}(0)=1+\theta +{{\theta }^{2}}+...{{\theta }^{N-1}}=\frac{(1-{{\theta }^{N}})}{(1-\theta )}.\end{eqnarray} \tag{ 25 }$$

Substituting in (22) we have our first main result, viz. the probability of establishment is given by:

$$\begin{eqnarray}&&{{P}_{N}}=\frac{(1-\theta ){{\theta }^{N-1}}}{(1-{{\theta }^{N}})}.\end{eqnarray} \tag{ 26 }$$

The temporal statistics can be obtained in a similar manner. From (21) we have that

$$\begin{eqnarray}&&\langle {{t}_{N}}\rangle =\frac{\mathop{\int }\limits_{0}^{\infty }{\rm d}t\ t\ {{Q}_{N}}(t)}{\mathop{\int }\limits_{0}^{\infty }{\rm d}t\ {{Q}_{N}}(t)}=-\frac{\hat{Q}_{N}^{^{\prime} }(0)}{{{{\hat{Q}}}_{N}}(0)}=\frac{J_{N-1}^{^{\prime} }(0)}{{{J}_{N-1}}(0)}\ .\end{eqnarray} \tag{ 27 }$$

Differentiating (20) with respect to s and setting s = 0 one finds

$$\begin{eqnarray}&&J_{n}^{^{\prime} }(0)=(1+\theta )J_{n-1}^{^{\prime} }(0)-\theta J_{n-2}^{^{\prime} }(0)+{{J}_{n}}(0)/n.\end{eqnarray} \tag{ 28 }$$

Solving by iteration for $J_{n}^{^{\prime} }(0)$ and using equation (27) gives an exact solution for the average time of colonization:

$$\begin{eqnarray}&&\langle {{t}_{N}}\rangle =\frac{1}{(1-\theta )}\mathop{\sum }\limits_{m=1}^{N}\frac{(1-{{\theta }^{m}})(1-{{\theta }^{N-m}})}{m(1-{{\theta }^{N}})}.\end{eqnarray} \tag{ 29 }$$

Following the same procedure, in solving by iteration for ${{J}_{n}}(0)^{\prime\prime}$, one finds after a lengthy calculation an exact solution for the variance in the time of colonization:

$$\begin{eqnarray}&&\begin{array}{l} \operatorname{var}({{t}_{N}})={{\langle {{t}_{N}}\rangle }^{2}}-\frac{2}{{{\left( 1-\theta \right)}^{2}}}\mathop{\sum }\limits_{m=2}^{N-1}\frac{(1-{{\theta }^{N-m}})}{m(1-{{\theta }^{N}})} \\ \qquad \qquad \,\times \mathop{\sum }\limits_{m^{\prime} =1}^{m-1}\frac{(1-{{\theta }^{m^{\prime} }})(1-{{\theta }^{m-m^{\prime} }})}{m^{\prime} }. \\ \end{array}\end{eqnarray} \tag{ 30 }$$

We are interested in micrometastases of moderate size, so it is useful to extract the asymptotic form of the temporal statistics for $N\gg 1$. Equation (26) gives the probability of establishment of a tumour of size N for any fitness θ. For $\theta <1$, but not too close to unity (i.e. $1-\theta \gg 1/N$), this equation takes the simple exponential form:

$$\begin{eqnarray}&&{{P}_{N}}\sim a{{\theta }^{N}},\end{eqnarray} \tag{ 31 }$$

with prefactor $a=(1-\theta )/\theta$. Quite intuitively, the probability of establishment decreases exponentially for increasing N. For example, a micro-tumour of 20 cells with fitness of $\theta =0.5$ has a probability of forming ${{P}_{20}}\sim {{(0.5)}^{20}}\sim {{10}^{-6}}$. Hence, one in a million progenitor cells will achieve this size for this fitness level.

The expressions for the average and variance in colonization times are too complicated to easily infer the dependence on the fitness θ and critical size N. By recasting the summations as integrals and performing an asymptotic analysis for $N\gg 1$ it is possible to extract the leading behaviour, which turns out to be quite simple. In these calculations, it is convenient to use the integral representation

$$\begin{eqnarray}&&\frac{1}{m}=\mathop{\int }\limits_{0}^{\infty }{\rm d}u\;{{e}^{-mu}}\end{eqnarray} \tag{ 32 }$$

in order to explicitly perform the summations, which have the form of a geometric series.

For the average colonization time we find, for $\theta <1$ and not too close to unity,

$$\begin{eqnarray}&&\langle {{t}_{N}}\rangle =\frac{{\rm log} N}{(1-\theta )}+\frac{(\gamma +{\rm log} (1-\theta ))}{(1-\theta )}+O(1/N),\end{eqnarray} \tag{ 33 }$$

where $\gamma =0.57721...$ is Eulerʼs constant [21]. Because of the logarithmic dependence of the average colonization time on N, we can rewrite this expression as a simple exponential growth law:

$$\begin{eqnarray}&&N=A(\theta ){\rm exp} \left( {{g}_{{\rm eff}}}\langle {{t}_{N}}\rangle \right),\end{eqnarray} \tag{ 34 }$$

where $A(\theta )={{e}^{-\gamma }}/(1-\theta )$ and the effective growth rate is ${{g}_{{\rm eff}}}=(1-\theta )$. Thus, we find that the average growth dynamics is exponential with a rate anti-correlated to the fitness. In other words, the less fit the cells are, the more aggressive the colonization dynamics of (rare) successful tumours will appear to be. On first reading, this result is counter-intuitive. One might imagine that the rare dynamics leading to colonies of low-fitness cells would be slow growing, relying on repeated chance avoidance of extinction. Despite the huge combinatorial weight of such meandering trajectories, the statistically relevant successful trajectories are those which grow rapidly, relying on repeated birth events and a relative absence of death events. This extreme bias becomes ever more important as $\theta$ decreases, hence the more rapid successful trajectories for cells of lower fitness.

Similar counter-intuitive results of this type are known in disparate fields, such as the 'bold play' strategy in gambling (successive 'all or nothing' gambling on high odds games is an optimal strategy for rare but large wins) [22], and transition rate theory in chemistry to quantify spontaneous dissociation of chemical bonds (stronger bonds break rarely but explosively). The Freidlin–Wentzell theory of large deviations provides a general mathematical framework to further understand such rare event dynamics, including the time-reversal duality discussed below [23]. We are obliged to ask whether this average time is a meaningful statistical measure. Perhaps the fluctuations in successful realizations are so great that the average time is a statistical oddity, with little bearing on the true nature of the dynamics. The variance of the distribution is informative in this regard.

Performing an asymptotic analysis of equation (30), and after a lengthy calculation, we find that the variance in colonization times is given by

$$\begin{eqnarray}&&\operatorname{var}\left( {{t}_{N}} \right)=\frac{\zeta (2)}{{{\left( 1-\theta \right)}^{2}}}+O(1/N),\end{eqnarray} \tag{ 35 }$$

where the Riemann Zeta function $\zeta (2)={{\pi }^{2}}/6$ [21]. Here we see that for large N the variance is independent of the critical tumour size and decreases with θ. This means that successful colonization trajectories are effectively more deterministic (i.e. less noisy) for lower fitness cells than for higher fitness cells. This is consistent with the preceding discussion, in that the statistically dominant growth trajectories of low-fitness cells are strongly biased to repetitive birth events, and will thus have a low degree of statistical variation.

In the limit of extremely low fitness, $\theta \to 0$, it is relatively straightforward to calculate the entire distribution of colonization times for $N\gg 1$, through a general solution of the generating equation (20). One finds that the colonization times are distributed according to the Gumbel distribution (well-known in extreme value statistics), centred at ${{t}^{*}}={\rm log} N$ (which, according to equation (33), is also the average colonization time in this limit). Explicitly, one finds for the distribution function of colonization times:

$$\begin{eqnarray}&&P({{t}_{N}})={\rm exp} \left( -\left( {{t}_{N}}-{{t}^{*}} \right)-{{e}^{-\left( {{t}_{N}}-{{t}^{*}} \right)}} \right).\end{eqnarray} \tag{ 36 }$$

This result suggests that the colonization time is drawn as an extreme value in samples of random variables with exponential distribution [24].

Despite the existence of exact solutions and asymptotic analysis, given the counter-intuitive nature of the results it is reassuring to check them against a numerical integration of the original master equation (1) with the boundary conditions equations (2)–(4). Integration was performed with a 2nd order Runge–Kutta scheme, and equations (6)–(9) were used to evaluate the statistics of interest. In figure 1 we compare the numerical solution to the exact results for the probability of establishment (26), the average time of colonization (29), and the standard deviation of the time of colonization (30). In each case we have perfect agreement as expected. We also compare a numerical evaluation of the exact results with the asymptotic analysis (for large N) in figure 2 for the average time of establishment (33) and the standard deviation of the time of colonization (35), and agreement is found, improving as N increases.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Numerical generation of stochastic realizations of the birth-death process are useful to get better intuition for the actual realizations of successful colonization events. This was achieved by implementing the Gillespie method [25] to generate many realizations of the birth/death stochastic process corresponding to the master equation (1). Results for the mean and standard deviation of establishment times are supported by these numerical results as shown in figure 3. Furthermore, we show some representative realizations of colonization for cells of two different fitnesses (figure 4) demonstrating explicitly that successful events arising from lower-fitness progenitors are indeed typically more explosive and deterministic in their dynamics than those for higher-fitness cells.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Given that our calculations for temporal statistics are valid for any value of $\theta$, it is interesting to compare our results for 'rare dynamics' (i.e. rare successful trajectories arising from cells with $\theta <1$), with results for 'rare genotypes' (i.e. successful trajectories arising from rare cells with $\theta >1$). So, we consider $\theta >1$ in the exact expressions for the growth statistics given in equations (26), (29) and (30). Taking the limit $N\gg 1$ we find that

$$\begin{eqnarray}&&{{P}_{N}}=1-\frac{1}{\theta }+O({{\theta }^{-N}}).\end{eqnarray} \tag{ 37 }$$

Thus for colonization from high-fitness cells, the fraction of successful colonies is of order unity as one would expect. Analysis of equations (29) and (30) shows the exact mapping:

$$\begin{eqnarray}&&\langle {{t}_{N}}\rangle {{|}_{\theta >1}}(\theta )=(1/\theta )\langle {{t}_{N}}\rangle {{|}_{\theta <1}}(1/\theta ),\end{eqnarray} \tag{ 38 }$$

$$\begin{eqnarray}&&\operatorname{var}({{t}_{N}}){{|}_{\theta >1}}(\theta )=(1/{{\theta }^{2}})\operatorname{var}({{t}_{N}}){{|}_{\theta <1}}(1/\theta ).\end{eqnarray} \tag{ 39 }$$

For $N\gg 1$, using the above relations with the asymptotic results for the temporal statistics found in the low-fitness regime, we have for the high-fitness colonization-time statistics:

$$\begin{eqnarray}&&\langle {{t}_{N}}\rangle =\frac{{\rm log} N}{(\theta -1)}+\frac{(\gamma +{\rm log} (1-1/\theta ))}{(\theta -1)}+O(1/N),\end{eqnarray} \tag{ 40 }$$

$$\begin{eqnarray}&&\operatorname{var}({{t}_{N}})=\frac{\zeta (2)}{{{\left( \theta -1 \right)}^{2}}}+O(1/N).\end{eqnarray} \tag{ 41 }$$

Interestingly, comparing the leading terms in these expressions with those for the low-fitness colonies, equations (33) and (35), we can see that the expressions are identical if ${{g}_{{\rm eff}}}=(1-\theta )$ for $\theta <1$ colonies is replaced by ${{g}_{{\rm eff}}}=(\theta -1)$ for $\theta >1$. The two solutions are 'dual' to one another for large N: the growth of a rare successful low-fitness colony, of fitness $\theta <1$, is statistically indistinguishable from the growth of a high-fitness colony of fitness $(2-\theta )>1$. This duality is illustrated in figure 5 where the exact expressions for mean and standard deviation of colonization times are plotted as a function of $\theta$ for $\theta \in (0,2)$.

Zoom In Zoom Out Reset image size

Relating this mathematical result back to the biological context of metastatic colonization leads to the remarkable conclusion: the dynamic progression of micro-metastases at early stages for rare dynamics and rare genotypes are indistinguishable. If one were to observe, in vivo, the early exponential growth of a micrometastasis, our results show that one would not be able to infer whether that event was seeded by a pre-adapted or non-pre-adapted cell. A distinction does occur for colony growth beyond the critical size N. Such colonies seeded by low-fitness cells are highly unlikely to sustain the rapid subsequent growth, and would enter a quasi-stable state in which adaptation and microenvironment reconfiguration will presumably enable further, and possibly slow, growth. By contrast, in colonies seeded by fit cells the critical size N is of little relevance, and exponential growth would presumably continue until access to nutrients becomes a limiting factor.

A duality of time-inversion is also implied in the exponential growth dynamics of non-pre-adapted cells given in (34). An initial tumour of size N comprising cells with fitness $\theta <1$ would most likely decay exponentially, with an effective decay rate of $1-\theta$. The dynamics would be the time-reversed dynamics of our rare event solution, a result reminiscent of reversible solutions in the Freidlin–Wentzell theory of rare event dynamics. Yet another duality exists, in fact, if one considers the dynamics of an initial tumour of size N comprising high-fitness cells ($\theta >1$), and asks for the temporal statistics of the rare event that such a tumour vanishes. The probability of such an event will be exponentially small, and will typically occur following an exponential decay, for the (time-reversed) reasons discussed at length above. This result implies that spontaneous remission of small tumours comprising fit cells, although rare, will occur rapidly and deterministically, an effect similar to the rare event study here, that would be interesting to pursue in future.

It is crucial to recognise that the fitness of a cell attempting colonization depends on both the cell phenotype and the nature of the secondary tissue. Clearly we expect significant heterogeneity in the millions of colonization initiation events happening during the progression of cancer, due both to the heterogeneity of cells leaving the primary tumour, and the heterogeneity of secondary tissue environments in which colonization is attempted. During later stages of metastasis genetic instability, adaptation and microevolution will induce cell heterogeneity within single metastatic tumours. Here though we are only concerned with the very early steps of metastatic colonization. Heterogeneity can then be described by introducing a distribution of fitnesses in the population of progenitor cells (i.e. those single DTCs attempting to proliferate to form a new colony). We denote this distribution ${{P}_{{\rm het}}}(\theta )$. Since this distribution involves properties of both DTCs and secondary tissue sites, it describes the attempted events of colonization throughout the body over the time course of metastatic disease.

For larger values of fitness, ${{P}_{{\rm het}}}(\theta )$ will be a rapidly decreasing function of $\theta$ because progressively fitter (or better adapted) cells should be progressively more rare in order to account for metastatic inefficiency. In this regard, it is helpful to introduce $\epsilon$ as the fraction of pre-adapted cells in the population:

$$\begin{eqnarray}&&\mathop{\int }\limits_{1}^{\infty }{\rm d}\theta \ {{P}_{{\rm het}}}(\theta )=\epsilon \ll 1,\end{eqnarray} \tag{ 42 }$$

based on the normalization condition

$$\begin{eqnarray}&&\mathop{\int }\limits_{0}^{\infty }{\rm d}\theta \ {{P}_{{\rm het}}}(\theta )=1.\end{eqnarray} \tag{ 43 }$$

The overall probability distribution for successful colonization events in the body is then given by the product of ${{P}_{N}}(\theta )$ and ${{P}_{{\rm het}}}(\theta )$. We will denote this distribution by ${{P}_{{\rm succ}}}$:

$$\begin{eqnarray}&&{{P}_{{\rm succ}}}(\theta )={{P}_{N}}(\theta ){{P}_{{\rm het}}}(\theta )=\frac{(1-\theta ){{\theta }^{N-1}}}{(1-{{\theta }^{N}})}{{P}_{{\rm het}}}(\theta )\ ,\end{eqnarray} \tag{ 44 }$$

where we have used equation (26). It is convenient to rewrite this expression as

$$\begin{eqnarray}&&{{P}_{{\rm succ}}}(\theta )=\frac{{{\theta }^{N-1}}}{\Sigma (\theta )}{{P}_{{\rm het}}}(\theta )\ ,\end{eqnarray} \tag{ 45 }$$

where $\Sigma (\theta )=1+\theta +...+{{\theta }^{N-1}}$.

The probability ${{P}_{N}}$ is an increasing function of θ, and the tail of ${{P}_{{\rm het}}}$ decays rapidly for large θ. Hence ${{P}_{{\rm succ}}}$ should have a maximum for some intermediate value $\theta ={{\theta }_{0}}$, corresponding to the fitness of the most abundant cohort of successful metastatic cells. By comparison, cells with $\theta \gg {{\theta }_{0}}$ are too rare in the initial population to have statistical significance in creating micro metastases, and likewise those with $\theta \ll {{\theta }_{0}}$ are too unfit. Therefore θ₀ represents the optimal balance between the rarity of initial phenotypes and the rarity of their corresponding dynamics to maximally contribute to nascent metastases. As such, if ${{\theta }_{0}}<1$, the colonization process is dominated by rare dynamics (figure 6(a)), whereas if ${{\theta }_{0}}>1$ colonization is dominated by rare genotypes (figure 6(c)). We can use the condition ${{\theta }_{0}}\equiv 1$ to define a phase boundary of the parameter space. This can be written as the condition:

$$\begin{eqnarray}&&{{\left. \frac{{\rm d}{{P}_{{\rm succ}}}}{{\rm d}\theta } \right|}_{\theta =1}}=0.\end{eqnarray} \tag{ 46 }$$

Combining equations (45) and (46), we find that the phase boundary condition provides a relationship between the threshold metastasis size N and the distribution of heterogeneity ${{P}_{{\rm het}}}$:

$$\begin{eqnarray}&&N=2{{\left| \frac{{\rm d}}{{\rm d}\theta }{\rm log} {{P}_{{\rm het}}}(\theta ) \right|}_{\theta =1}}+1\ .\end{eqnarray} \tag{ 47 }$$

Zoom In Zoom Out Reset image size

To proceed further, it is necessary to make an explicit choice for the form of ${{P}_{{\rm het}}}$. Unfortunately this distribution is not available experimentally, but one can make a pragmatic choice nonetheless. First, the details of the distribution for low fitness values are not very important since, from equation (47), one will be differentiating the distribution at the threshold fitness value of unity. Second, due to the experimental fact that metastasis is a very inefficient process [4], a reasonable assumption is that the tail of the distribution decreases exponentially. Such an exponential form could arise, for example, if accumulating discrete random events were required to endow cells with higher fitness values; an idea consistent with popular conceptions of metastatic potential of cells arising from successive mutations of key genes [26].

We therefore take ${{P}_{{\rm het}}}$ to have a generic exponential form:

$$\begin{eqnarray}&&{{P}_{{\rm het}}}(\theta )=A{\rm exp} (-b{{\theta }^{s}}).\end{eqnarray} \tag{ 48 }$$

The normalization condition (43) fixes $A=s{{b}^{1/s}}/\Gamma (1/s)$, where $\Gamma (z)$ is the Gamma function [21]. The second integral condition (42) relates b and s to as follows

$$\begin{eqnarray}&&\frac{{{b}^{1/s-1}}{{e}^{-b}}}{\Gamma (1/s)}=\epsilon .\end{eqnarray} \tag{ 49 }$$

Since high-fitness cells are rare, we have $\epsilon \ll 1$. This enables us to give an approximate solution to the transcendental equation (49):

$$\begin{eqnarray}&&b={\rm log} (1/\epsilon ).\end{eqnarray} \tag{ 50 }$$

Note this solution is exact if ${{P}_{{\rm het}}}$ is a pure exponential function (s = 1).

It is now straightforward to insert the explicit form of ${{P}_{{\rm het}}}(\theta )$ (equation 48) into the condition (47) to derive a relationship between N and . Taking $(s=1)$ for simplicity we get:

$$\begin{eqnarray}&&N=2b+1=2{\rm log} (1/\epsilon )+1.\end{eqnarray} \tag{ 51 }$$

This relation demarcates the two domains in the $(\epsilon ,N)$ parameter space as shown in figure 6(b). Specifically, rare dynamics dominates when the fraction of high-fitness cells is small and the critical size N is not too large. The impact on values of N defining the boundary is only weakly (logarithmically) dependent on .

Statistical arguments to follow will allow us to estimate relevant sizes for N. First though, in this subsection, we provide here a simple argument for the lower limit on the size of a nascent metastatic colony which will allow a protective microenvironment for non-surface cells within the cluster. This argument is mainly illustrative, relying purely on geometrical considerations. The true efficacy of a microenvironment depends on many factors such as the biochemical and signalling milieux of the tumour, as well as the morphology of its constituent cells. In addition, for the cluster to be stable, the rate of division of core cells must be high enough to at least replace surface cells which are lost to cell death. This balance will depend on the relative fitnesses of core and surface cells, which is a biological aspect we do not account for in the geometrical arguments below.

Assuming that cells in the metastasis are spherical and of identical sizes, the minimal number of cells required to shield one single core cell is 12, since this is the 'kissing number' for spheres in three dimensions [33]. If indeed a surface of 12 cells provides a sufficiently beneficial microenvironment to the single core cell, one can assign this core cell a high fitness. The subsequent dynamics are, however, not particularly robust, since the progeny of a single high-fitness cell are at significant risk of stochastic extinction. A more robust core would require more than one cell. Using simple geometrical considerations, and assuming a large cluster, one finds that the number of surface cells N_s required to cover N_i internal core cells is given by the relation

$$\begin{eqnarray}&&{{N}_{s}}\sim {{\left( N_{i}^{1/3}+2{{e}^{1/3}} \right)}^{3}}-{{N}_{i}},\end{eqnarray} \tag{ 52 }$$

where e is the packing efficiency, a value that typically varies between $e\sim 0.52$ for randomly packed spheres and $e\sim 0.74$ for Gaussʼs optimal sphere packing. According to this formula a core comprising 2–10 cells (allowing a more deterministic growth dynamics) would require approximately 25–50 surface cells (figure 7). A more accurate estimate can be gained from numerically packing spheres according to face-centred-cubic close packing, which produces somewhat lower estimates that equation (52). In this case 20-40 surface cells are required to enclose 2–10 core cells. Thus, on geometric grounds, the minimum stability threshold for a proto-metastasis, with greater than one core cell, is around $N\sim 20$–50 cells.

Zoom In Zoom Out Reset image size

It is possible that the stability threshold could be smaller if surface cells are afforded some degree of protection also. Since cells are likely to adapt their morphology to that of their neighbours, via adhesive interactions, one can think of surface cells in the tumour forming a contiguous surface layer having only one facet in direct contact with the secondary tissue. This partial screening may be sufficient to raise the fitness of surface cells above unity, in which case all cells in the proto-metastasis (both core and surface) would be capable of deterministic growth. In this case, the minimum stability threshold might be as small as approximately 12 cells.

To proceed further, in discussing the relevance of our theoretical results, it is necessary to make closer contact to human cancer, and infer what we can from the dearth of empirical data on the dynamics of metastatic colonization in humans. Indeed, little is known about the properties of DTCs in human cancer. Even measures such as the typical rate of shedding of cancer cells from the primary tumour are poorly known. As such, physiological parameters describing the number and diversity of DTCs will be subject to high uncertainties. Fortunately, for our purposes, it will transpire that only the logarithm of these parameters will be relevant to calibrating our model, so that approximate orders of magnitude in the physiological parameters will suffice.

For concreteness, consider individual DTC colonization attempts occurring over a period of one year in a hypothetical cancer patient who is at risk of metastatic disease, but not in late stage disease from previous successful metastatic colonization. The total number of such attempts, which we denote by M, will be very large, ranging from tens of thousands to hundreds of millions. The typical concentrations of CTCs found in patients with metastasis disease can be estimated at about one to ten cells per millilitre of blood [11, 27], which corresponds to about 10⁴–10⁵ total CTCs in the human body at any given time. Although a fraction of CTCs will die in the bloodstream, we expect that a significant fraction are rapidly arrested and extravasated from the vascular system as shown by Luzzi et al [4] (note, this is an extrapolation from a mouse model), and most of them are not in the blood stream for more than 24 hours [11, 28]. Thus, over a one year period we can expect approximately 10⁶–10⁷ different DTCs to be attempting colonization. This estimate appears low when compared to other experimental sources. For example, Butler et al [29] showed that approximately one million cells per day are shed per one gram of tumour in a rat mammary carcinoma model. If this figure were translated to humans, a one gram tumour would yield over two orders of magnitude greater number of DTCs than quoted above. A lower bound for this number comes from studies of DTCs in bone marrow of prostate tumour patients. In such cases it is estimated that 10⁴–10⁵ DTCs may be in evidence, most likely in a dormant state [30]. Given this great uncertainty, for our purposes we take M in the range 10⁴–10⁸.

Assuming that rare dynamics is the dominant mechanism for colonization, we can crudely estimate a range of values for N as follows: let m be the number of colonization attempts in the relevant range ${{\theta }_{0}}\pm \Delta {{\theta }_{0}}$ where, ultimately, successful events are likely to emerge. It is not possible to estimate systematically what fraction of the M attempts will contribute to m, so given the enormous range of uncertainty in M, we will assume that the range of m is equally large and uncertain. We denote by K the number of m attempts resulting in a metastatic colony. The great inefficiency of metastasis implies that successful colonization events are exceedingly rare [4] and as such K will be orders of magnitude smaller than m. From the definitions given above we have the direct relationship

$$\begin{eqnarray}&&m{{P}_{N}}({{\theta }_{0}})\simeq K.\end{eqnarray} \tag{ 53 }$$

Assuming that θ₀ is not too close to unity, equation (31) gives us ${{P}_{N}}({{\theta }_{0}})\sim a\theta _{0}^{N}$, which together with (53) allows us to express the stability threshold N in terms of the other quantities, viz.

$$\begin{eqnarray}&&N\sim \frac{{\rm log} \left( m/K \right)}{{\rm log} \left( 1/{{\theta }_{0}} \right)}\sim \frac{{\rm log} \left( m \right)}{{\rm log} \left( 1/{{\theta }_{0}} \right)},\end{eqnarray} \tag{ 54 }$$

where we have used in the second step the fact that $m\gg K$. On varying m within the range 10⁴–10⁸ and θ₀ within the range 0.25–0.85 one finds that N takes a modest range of values, between 10 and 100. The fact that the range of N is relatively well-defined given the enormous uncertainty in the number of attempts and the fitness of DTCs is due to the logarithmic dependence of N on these parameters.

There is a remarkable concurrence of this size range with both geometrical reasoning and biological observations. As discussed in the previous subsection geometrical estimates on the size of a tumour which would enable protection of the core from a hostile microenvironment yields a range for N of about 20–50 cells. Observations of induced metastatic colonization in zebrafish reported that clusters as small as 15–30 cells were able to induce angiogenesis for subsequent growth [31]. There have also been reports of stable metastatic cell clusters in mice of size 30–60 cells [32]. We believe the congruence of biological observations with geometrical considerations and our rare event analysis in identifying 10–100 cells as a threshold size for stability argues for closer experimental examination of such minute clusters, which we term 'proto-metastases'.

We have used equation (54) to determine the possible range of N given crude estimation of physiological ranges of θ₀ and m, finding this range to be 10–100. We can take a different strategy and use (54) to provide a range of values of θ₀ assuming a value for N motivated by biological data, e.g. N = 30 as discussed above. Thus, taking N = 30 in (54), and the physiological range $m\sim {{10}^{4}}$–10⁸, we find that θ₀ takes values in the rather confined range 0.55–0.75. The absolute probability of colonization for a given event will be trivially given by K/m, and thus will lie in the large range of 10⁻⁸–10⁻⁴, which simply reflects the large uncertainty in the number of DTCs. Using the range of θ₀ and N = 30 we can estimate the average time for a given successful colonization event using the leading term of equation (33). This equation gives the average time in units of the mean time (the inverse of the Poisson rate) of cell death, which we crudely estimate for non-dormant DTCs to be of order one day. One then finds $\langle {{t}_{30}}\rangle$ in the range 7.5–13.5 days, i.e. between one and two weeks. This time-scale of initial growth of a successful proto-metastasis is rapid compared to the time interval separating successful events, and relatedly, the typical time-scale of years over which human metastasis clinically progresses.