Summary of the supplementary notes from https://www.nature.com/articles/s41588-018-0128-6 by Williams et al. 2018
To model tumor growth, let us start with the assumption that we have 1 tumor cell initially:
After $\epsilon$ time has passed, each tumor cell can:
- Do nothing
- Divide Once into 2 Daughter Cells
- Undergo Cell Death
- Some combination of dividing and cell death
When $\epsilon \to 0$,
Let $b \epsilon$ be the probability of dividing once, and $d \epsilon$ be the probability of cell death
Any combination of dividing and cell death in time $\epsilon$ is therefore $\leq (b \epsilon) (d \epsilon) = O(\epsilon^2) \approx 0$
Let $D_i (t)$ be the number of cells at time $t$ which have undergone $i$ divisions.
After a cell death, one is subtracted from $D_i(t)$
After cell division, one is subtracted from $D_i(t)$ and $\textbf{two new cells}$ are added to $D_{i+1}(t)$
On average after time $\epsilon \to 0$, $\epsilon b N(t)$ cells divide and $\epsilon d N (t)$ cells die
Therefore $lim_{\epsilon \to 0} D_i(t+\epsilon) \approx lim_{\epsilon \to 0} [D_i(t) − (\epsilon b) D_i (t) − (\epsilon d) D_i(t) + \epsilon (2b) D_{i−1} (t)]$
$lim_{\epsilon \to 0}[D_i(t+\epsilon)−D_i(t)] \approx lim_{\epsilon \to 0}[−\epsilon(b+d)D_i(t)+\epsilon(2b)D_{i−1}(t)]$
$\frac{d}{dt} D_i(t) = lim_{\epsilon \to 0} \frac{D_i(t+\epsilon) − D_i(t)}{\epsilon} = −(b+d)D_i(t) + (2b)D_{i−1}(t)$
$\frac{d}{dt} D_0(t) = lim_{\epsilon \to 0} \frac{D_0(t+\epsilon) − D_0(t)}{\epsilon} = −(b+d)D_0(t)$
When there is only one initial cell, $D_0(0) = 1$ and $N(0)=1$:
$D_0(t)=e^{−(b+d)t}$
$D_i(t)=\frac{(2bt)^i}{i!}e^{−(b+d)t}$
Proof:¶
Base Case:¶
$\frac{d}{dt} e^{−(b+d)t} = −(b+d) e^{−(b+d)t}$
$\implies D_0(t) = e^{−(b+d)t}$
Inductive Hypothesis:¶
Assume for some $i−1 \geq 0$, $D_{i−1}(t) = \frac{(2bt)^i}{(i−1)!} D_0(t)$ satisfies the differential equation
Recursive Case:¶
$\frac{d}{dt} D_i(t) = 2b D_{i−1}(t)−(b+d) D_i(t)$
By the IH, $\frac{d}{dt} D_i(t) = 2b \frac{(2bt)^{i−1}}{(i−1)!} e^{−(b+d)t} −(b+d)D_i(t)$
$= i\frac{(2b)^i t^{i−1}}{(i)!} e^{−(b+d)t} − (b + d)\frac{(2bt)^i}{i!} e^{−(b+d)t}$
By the Reverse Chain Rule
$\frac{d}{dt} D_i(t) = \frac{d}{dt}[\frac{(2bt)^i}{i!} e^{−(b+d)t}]$
QED¶
Normalizing by the Total Number of Cells¶
The total number of cells is given by $N(t)=e^{(b−d)t}$
Dividing the previous equation by $N(t)$ gives the $\textbf{density}$ (proportion of cells) which have divided $\Gamma$ times.
$P_i(t)=\frac{(2bt)^i}{i!} \frac{e^{−(b+d)t}}{e^{(b−d)t}} = \frac{(2bt)^i}{i!} e^{−2bt}$
$P_0(t)=e^{−2bt}$
Naive Forward Simulation¶
import numpy as np
import math
import matplotlib.pyplot as plt
n_cells_stop = 2 ** 20
n_divisions = np.zeros(100).astype(np.int32)
n_divisions[0] = 1
b = 0.000001 # cell division rate
d = 0.0000001 # cell death rate
while n_divisions.sum() < n_cells_stop:
event = np.random.choice(['cell_birth', 'cell_death'], p=[b/(b+d),d/(b+d)]) # Assuming an event happens (otherwise nothing changes)
cell_chosen = np.random.choice(len(n_divisions), p=n_divisions / n_divisions.sum())
n_divisions[cell_chosen] -= 1
if event == 'cell_birth':
n_divisions[cell_chosen+1] += 2
Theoretical Prediction¶
def theoretical_predicition(n_cells_stop, b, d):
t = np.log(n_cells_stop) / (b-d)
P = np.zeros(n_divisions.shape)
for i, _ in enumerate(P):
P[i] = np.exp(-2*b*t) * (np.power(2 * b * t, i) / math.factorial(i))
return P
plt.ylabel('Density')
plt.xlabel('Number of Divisions')
plt.bar(np.arange(len(n_divisions)), n_divisions/n_divisions.sum(), color='black')
plt.plot(theoretical_predicition(n_cells_stop, b, d), color='red', linewidth=3, alpha=0.5, label='Theoretical Prediction')
plt.legend()
<matplotlib.legend.Legend at 0x117a78e20>
