Quantitative Interpretation of a Genetic Model of Carcinogenesis Using Computer Simulations

The genetic model of tumorigenesis by Vogelstein et al. (V theory) and the molecular definition of cancer hallmarks by Hanahan and Weinberg (W theory) represent two of the most comprehensive and systemic understandings of cancer. Here, we develop a mathematical model that quantitatively interprets these seminal cancer theories, starting from a set of equations describing the short life cycle of an individual cell in uterine epithelium during tissue regeneration. The process of malignant transformation of an individual cell is followed and the tissue (or tumor) is described as a composite of individual cells in order to quantitatively account for intra-tumor heterogeneity. Our model describes normal tissue regeneration, malignant transformation, cancer incidence including dormant/transient tumors, and tumor evolution. Further, a novel mechanism for the initiation of metastasis resulting from substantial cell death is proposed. Finally, model simulations suggest two different mechanisms of metastatic inefficiency for aggressive and less aggressive cancer cells. Our work suggests that cellular de-differentiation is one major oncogenic pathway, a hypothesis based on a numerical description of a cell's differentiation status that can effectively and mathematically interpret some major concepts in V/W theories such as progressive transformation of normal cells, tumor evolution, and cancer hallmarks. Our model is a mathematical interpretation of cancer phenotypes that complements the well developed V/W theories based upon description of causal biological and molecular events. It is possible that further developments incorporating patient- and tissue-specific variables may build an even more comprehensive model to explain clinical observations and provide some novel insights for understanding cancer.


Introduction
In the following figures we present several trajectories for simulations of the stochastic model. For each experiment, 101 simulations were performed. The trajectories were ranked by cell number within the mass at the time the simulations were terminated, with two exceptions. Some of the masses disappeared (cell number = 0) before termination. These masses are ranked by the order in which they disappeared. Some of the mases grew above a preset limit before the time interval of interest elapsed. These masses are ranked by the order in which they surpassed the preset limit. The trajectories presented in this supplemental correspond to those at rank 1, 25, 51, 76, and 101. That is, they represent the minimum terminal result, maximum terminal result, the median, and the results at the 25th and 75th percentiles.

Results
Flowchart for the evolution of a cell by the stochastic model.
The following step-by-step algorithm can be used to advance the cell number over time. This is the algorithm for the advancement of one cell by one time step ∆t. This would have to be repeated for each cell in the system. Over the course of one time step (∆t), the system is not stochastic, and therefore the differential equation for α(t) can be solved explicitly. We can then derive expressions for the advancement of N (exp) (t) = t t0 α(s)ds where t 0 is the time of birth of the cell, N (t) = 2 N (exp) (t) , I(t) = t 0 |α(t)|dt, and g(t) = ceiling I(t) . Let p denote the number of mutations gained per generation. 1. t new = t old + ∆t 2. Assign β. 3. Evaluate cellular status: -If 'yes', then -If 'no', then new ≥ 1, then wait to update k with mutations until step #9.

Is
• If 'yes', the cell dies.
• If 'no', then return to #1. 7. The cell has proliferated into 2 f loor(N (exp) new ) daughter cells. Replicate 2 f loor(N (exp) new ) − 1 new cells. 8. Assign each daughter cell α, g, I, k, α p , and k p from the parent cell. 9. For each cell: 10. Return to #1 for each daughter cell, and run the algorithm for each separately.
Note, however, that as we proceed through the algorithm in practice, we do not need to explicitly keep track of the values of N (t), g(t), α p (t), or k p (t).
The clone lifetime of uterine epithelial cells during menopause.
In this experiment, we examine the effect of microenvironmental stimulation on normal clonal growth.
We begin with one cell, with initial conditions Experiments are performed with sustained β = 0, β ∼ N (0, 1), β ∼ N (0, 2), β ∼ N (1, 1), and β ∼ N (2, 1) in the absence of mutations. Results are shown in Figure S1. There does not appear to be much deviation from the median behavior. However, it is clear that β can have the effect of increasing mass size, but shortening its lifetime.
The clone lifetime of uterine epithelial cells in the presence of oral contraceptives.
In this experiment, we examine the effect of hormonal stimulation on normal clonal growth, a situation mimicking that occuring via oral contraceptives. We begin with one cell, with initial conditions Experiments are performed with β ∼ N (60, 6), β ∼ N (90, 9), and β ∼ N (120, 12) applied for 21 days, followed by sustained β = 0. All experiments are performed in the absence of mutations. Results are shown in Figure S2. There again does not appear to be much deviation from the median behavior. However, β can have a massive effect on overall mass size.
Tumor evolution: growth of a heterogeneous mass consisting of cells with various differentiation coefficients.
In this experiment, we examine the evolution of a heterogeneous mass composed of mature cells. Clone size is illustrated in Figure S3. The median k-values of these masses are shown in Figure S4, illustrating the evolutionary advantage of cells with lower k-values. We start with a mass of 1000 cells, with initial conditions α p (0) = 0, α(0) ∼ N (0, 1).
Experiments are performed with k = 0, k ∼ N (2, 1), k ∼ N (3, 1), and k ∼ N (4, 1). All cells are exposed daily to β ∼ N (0, 1). If k < 0 for some cell, we reassign it k = 0. Experiments are performed in the absence of mutations and without advancing generations (g(t)). We can see from Figure S3 that there is much deviation from the median result, either in the cell number after two years, or the rapidity with which the masses approach the simulation limit of 10 12 cells. For mass sizes above 10 7 , we approximate the size by setting β = 0 and analyzing the resultant system of differential equations for clone growth. We can see from Figure S4 that cells with lower k-values have an evolutionary advantage ( dα dt = −kα + β) and are selected for as the mass progresses. Note that the masses that reached the simulation limit size of 10 12 cells had median k = 0.

Survival of metastatic cells in the circulation and establishment of metastatic lesions.
In this experiment, we examine cells undergoing the process of metastasis. We represent the general steps of the dislocation process with microenvironmental effects (β) representative of the relative difficulty of each step. We start with 10,000 cells, with initial conditions Experiments are performed with k = 0, k = 1, and k = 2. Cells are exposed to β ∼ N (−6, 0.2) for 3 days representing dissemination at the primary site, β ∼ N (−10, 0.2) for 2 days representing intravasation, β ∼ N (−20, 0.2) for 1 days representing circulation, β ∼ N (−10, 0.2) for 2 days representing extravasation, β ∼ N (−6, 0.2) for 3 days representing survival at the primary site. This is followed by β = 0 for 1800 days. Experiments are performed in the absence of mutations and without advancing generations (g(t)). Results are shown in Figure S5. We can see from the figure that there does not seem to be much deviation from the median behavior. However, we can see that, although higher k-values permit cells to survive the metastasis to the ectopic location, it is the cells with lower k-values that thrive in the new environment, producing masses.

Discussion
One observation from these simulations is that, although the environmental stimulation β can have a large effect on mass size, it is a consistent effect. That is, similar cells have similar outcomes; there is not much variation from median behavior. However, as we can see from Figure S3, variations in k can lead to wildly different behaviors among similar cells in similar environments. This dynamic between the effects of variations in β and k in this mathematical model may lead to many interesting results. Therefore, one point of future research will be to provide a global sensitivity analysis of the mathematical model along with biological interpretations of the results.