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{(2 b t)^i}{i!}e^{-(b+d)t}$
Proof:¶
$\textbf{Base Case:}$
$\frac{d}{dt} e^{-(b+d) t} = -(b+d) e^{-(b+d)t}$
$\implies D_0(t) = e^{-(b+d) t}$
$\textbf{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
$\textbf{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}]$
$\textbf{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 0x7514f13fb500>
