This model employs the multi-stage carcinogenesis concept. As tissues might require anywhere between 2 to 8 driver mutations (denoted by S) for malignant transformation (37), we first identify different cell types that can lead to a malignant transformation based on number of driver mutations. Let us visualize the dynamics of 5 types of cells in a tissue ( Figure 1 ). “Type 0”, “Type 1”, “Type K”, “Type S-1” and “Type S” represent normal healthy cells with no mutation; premalignant cells with one cancer-related mutation; premalignant cells with K cancer-related mutation, premalignant cells with S-1 cancer-related mutations and cancer cells with S cancer-related mutations, respectively. Emergence of cancer cells must be preceded by that of premalignant cells with mutations from Type 1 cell to Type S-1 cell. Type K cell may or may not be present depending on number of mutations required for carcinogenesis. We assume that a normal healthy tissue consists of Type 0, Type 1, Type K and Type S-1 cells undergoing cellular turnover with a small probability of a mutation. Moran process is employed to consider the tissue turnover dynamics, where the total number of Type 0, Type 1, Type K and Type S-1 cells is kept constant as N (39). The turnover rate of a whole tissue is defined by d. Type S cells are considered as uncontrolled, highly proliferating cancer cells. The branching process is employed to consider the process of Type S proliferation (40).

Initially, N Type 0 cells occupy the tissue. There is a rare chance of a mutation every time a cell divides, and a daughter cell may change into a Type 1 cell with a mutation rate, μ ₁. Mutation rate, μ, refers to the sum total of the genomic or epigenetic factors affecting change from one cell type to another (41). When a cell dies, a cell to be divided in a tissue is selected depending on the cell fitness, r. The fitness of Type 0, Type 1, Type K, and Type S-1 are denoted by r ₀, r ₁, r_K and r_{S -1}, respectively. Cell fitness, r refers to the transcriptional and metabolic potential of a cell type to “out-compete” other cell types (42). A cell could divide to give rise to the same cell type or mutate to another cell type. When a Type 1 cell divides with a mutation, a daughter cell may change into a Type K cell with mutation rate μ_K (if more mutations are needed) or a Type S-1 cell with a mutation rate, μ_{S -1} (if additional mutation steps are not needed). Intermediate cell type, Type K, becomes Type S-1 after sequential accumulation of mutations. Finally, a Type S-1 cell is capable of mutating to become a malignant cell – Type S cell based on the mutation rate from Type S-1 to Type S cell, μ_S . Once a Type S cell appears, the cells proliferate indefinitely based on the growth rate of Type S cells, r_S , disrupting tissue dynamics and homeostasis. Type S cancer cells are “super-competitors” with outstanding metabolic prowess and assumed to increase exponentially with a net growth rate of r_S - d_S > 0, where d_S is a death rate.

We propose that the most important premalignant cells are the Type 1 cell that have acquired the first driver mutation and the Type S-1 cell that needs just one more driver mutation to become a cancer cell. These cells look phenotypically normal and are not regarded as important clinically but their genetic features are indispensable in cancer formation. As a result, the cell fitness of these 2 cell types must be taken into account in all computational analysis. For scenarios where additional mutations and more cell types are needed, we approximate intermediate mutational steps between Type 1 and Type S-1 by adjusting the values of μ_S-1 so that a low mutation rate, μ_S-1 , represents additional mutational steps. So, our computational analysis will be executed to account for the most important mutational events that affect cell fitness and all the mutation rates that can affect the number of steps required for carcinogenesis. In other words, we skip the state of Type K cell and a mutation rate stands in for the number steps.

The net growth of Type 0, Type 1, and Type S-1 cells is zero (equal frequency of cell division and death), while that of Type S cells is positive. Type 0, Type 1, and Type S-1 cells consist of a healthy tissue based on the Moran process, so r₀ , r₁ and r_S-1 are parameters to determine status of dividing cell and which daughter cell are obtainable at the time of a cell division. Alternatively, r_S is the growth rate, which determines the average number of increases in Type S cells during a unit time. When the number of Type S cells reaches 10⁹ at the first time, all the Type S cells are discarded to represent surgical resection, whereas the number of Type 0, Type 1 and Type S-1 cells in a tissue is preserved so that the time until the emergence of the recurrent tumor is influenced by the frequency of residual Type 0, Type 1 and especially, Type S-1 cells. Since the conversion from the number of cells to the tumor volume is frequently done using the following relationship as 10⁹ cells in a 1 cm³ tumor, the time of surgery in this model is conducted when the size of the tumor becomes 1 cm³. We describe it as “time of cancer detection”. After the first treatment, the simulation continues until the next Type S cell appears from the tissue and the number reaches 10⁹ again, representing the recurrence of the tumor after surgery.

To integrate the Moran process and branching process, we adopted stochastic simulations based on Gillespie’s algorithm (43) as follows: We firstly considered three events: (i) cell turnover in a healthy tissue as per Moran process, (ii) birth of a Type S cell as per Branching process, and (iii) death of a Type S cell. The rates of each event at time t is given by (i) dN (ii) r_SX_S (t), and (iii) d_SX_S (t), respectively. Here r_S , d_S , and X_S (t) are a proliferation rate, a death rate, and the number of Type S cells at time t, respectively. Then an average time until one of the three events happens, ∆ T , is given by

Let us first consider the case where a cell turnover happens. The probability that a cell turnover happens in ∆ T is given by d N ∙∆ T . In our model, a cell turnover in a healthy tissue is governed by a cell death. When one of N cells is randomly selected as a cell to die, and another cell is chosen to divide within the same time step to complete cell turnover. In a healthy tissue, there are three types of cells, corresponding to the number of acquired mutations, Type 0, Type 1, and Type S-1.The number of each cell type is denoted by X ₀, X ₁ and X_{S -1}, respectively. In brief, there are several possibilities of tissue composition transitions in the tissue dynamics and we consider the six events that affect the cell type composition of a tissue: (i) a type 0 cell increases by one while a type 1 cell decreases by one (ii) a type 0 cell increases by one while a type S-1 cell decreases by one (iii) a type 1 cell increases by one while a type 0 cell decrease by one (iv) a type 1 cell increases by one while a type S-1 cell decreases by one (v) a type S-1 cell increases by one while a type 0 cell decreases by one; or (vi) a type S-1 cell increases by one while a type 1 cell decreases by one. In such a condition, a Type 0 cell can increase by one if either a Type 1 or Type S-1 cell dies and a Type 0 cell divides without a mutation. Then the probability for these events leading to an increase in Type 0 cells are given by (i) X 1 N ∙ r 0 X 0 ( 1 − μ 1 ) F , and (ii) X S − 1 N ∙ r 0 X 0 ( 1 − μ 1 ) F . Here, F = r 0 X 0 + r 1 X 1 +   r S − 1 X S − 1 is a scaling factor for selecting a dividing cell. The probability of a Type 1 or Type S-1 cell death is given by X 1 N and X S − 1 N , respectively. Taken together, the transition probability that the number of Type 0 cell increases by one and that of Type 1 decreases by one is given by

and the probability that the number of Type 0 cell increases by one and that of Type S-1 decreases by one is given by

A Type 1 cell can increase by one if either a Type 0 or Type S-1 cell dies, and either a Type 1 cell divides without mutation or a Type 0 cell divides with mutation to become a Type 1 cell. Then the probabilities for these events leading to an increase in Type 1 cells are given by (iii) X 0 N ∙ r 1 X 1 ( 1 − μ S − 1 ) + r 0 X 0 μ 1 F , and (iv) X S − 1 N ∙ r 1 X 1 ( 1 − μ S − 1 ) + r 0 X 0 μ 1 F . Taken together, the transition probability that the number of Type 1 cell increases by one and that of Type 0 decreases by one is given by

and the probability that the number of Type 1 cell increases by one and that of Type S-1 decreases by one is given by

Similarly, a Type S-1 cell can increase by one if either a Type 0 or Type 1 cell dies, and either a Type S-1 cell divides without mutation or a Type 1 cell divides with mutation. The probabilities for the events leading to an increase in Type S-1 cells are given by (v) X 0 N ∙ r S − 1 X S − 1 ( 1 − μ S ) + r 1 X 1 μ S − 1 F , and (vi) X 1 N ∙ r S − 1 X S − 1 ( 1 − μ S ) + r 1 X 1 μ S − 1 F . Taken together, the transition probability that the number of Type S-1 cell increases by one and that of Type 0 decreases by one is given by

and the probability that the number of Type S-1 cell increases by one and that of Type 1 decreases by one is given by

In addition, a Type S cell can increase by one if a Type S-1 cell divides with mutation. The probability is given by r S − 1 X S − 1 μ S F . Since a Type S cell is not a component of a tissue, once a Type S appears by mutation, another round of selection for a dividing cell is performed according to the transition probabilities described above. This is because malignant Type S cell disrupts 2D lattice structure and the Moran process is no longer applicable to it.

Next, let us consider the case where Type S cell divides or dies. The probabilities of Type S cell division or death is given by r S X S ( t ) ∙ ∆ T and d S X S ( t ) ∙ ∆ T , respectively.

In summary, the time of one step in our simulation is calculated using Eq. (1) and in one time step, one of the following three processes occurs: (i) a cell turnover in a tissue, (ii) the birth of a Type S cell, or (iii) the death of a Type S cell. Initially, all the cells are Type 0. Once the number of Type S cells reaches 10⁹, computational surgical resection sets the number of Type S cells to be 0, keeping the cell type composition in a tissue remained and computational carcinogenic process restarts again. After that, the time until the number of Type S cells reaches 10⁹ is measured as recurrence time.

Two-dimensional lattice structure ( I × J ) is introduced to a tissue dynamics in our computational model. The transition probabilities are basically the same with or without spatial structure. The difference is the choice of a dividing cell. If a cell at position 〈 i , j 〉 dies, 4 adjacent cells – 〈 i , j − 1 〉 , 〈 i , j + 1 〉 , 〈 i − 1 , j 〉 and 〈 i + 1 , j 〉 can divide to replace it. The transition probabilities are calculated according to the cell type at those positions. We assume wall boundary condition to represent an asymmetric tissue structure.

As for the calculation of the Type S growth, we assume that when the number of cells is small, the stochastic effect should be considered. When the number of Type S cells exceed twice as large as the size of the normal tissue, 2N, growth can be regarded as a deterministic process. Then, the time duration from when the number of Type S cells is 2N to 10⁹, ∆t _s, is given by

The data used in our analysis were from TCGA Pan-Cancer Clinical Data Resource (34, 35) and are available in the cBio Cancer Genomics Portal (44, 45). We adopt the clinical data of locoregional recurrence from 8,957 patients with 27 different non-sarcoma, non-hematological cancer types. From these datasets, the inclusion criterium for our study was “disease-free” survival – patients with no detectable malignant disease after surgery or total remission. We excluded data of “progression-free” survival in order to eliminate patients who survived with detectable disease possibly as a result of treatment-resistant clones; and also excluded data containing metastatic progression. We also included data from other independent publications for extra validation. Sarcomas and hematological cancers were excluded due to their non-conformity to a 2-dimensional lattice structure.

Survival time analysis of clinical data is calculated using the Kaplan–Meier method from disease-free intervals mentioned in Clinical Data section. In this study, disease-free interval is defined as the survival time without cancer recurrence of each patient, which corresponds to the time to recurrence of each simulation trial.

The whole process of our model is conducted on C++. Simulation codes have been deposited in a GitHub repository (https://github.com/sharaf501/Heano-Lab-Codes). The survival time analysis and other statistical analysis is conducted on Prism (version 9.4.1). Mantel-Cox (log-rank) test is used to compare difference between survival curves. A p value less than 0.05 is considered to be statistically significant.

Firstly, we conducted stochastic triplicate simulations for the cancer initiation up to the time of cancer detection. We were curious to know what effect the presence or absence of the spatial structure would have on the model. We traced the time course of 4 cell populations – Type 0, Type 1, Type S-1 and Type S cells using a combination of various parameter sets. Lower mutation rate from Type 1 to Type S-1, μ_{S -1}, was additionally examined to account for additional premalignant cell types between Type 1 and Type S-1. In the model without spatial structure, we observed 3 patterns of cancer initiation based on frequency of non-malignant cell population at cancer detection ( Figure 2A ). Interestingly, all the patterns show a progressive decline of Type 0 cells until the entire tissue is dominated by Type 1, Type S-1 or both Type 1 and Type S-1 cells. By combining various parameter sets in our simulation, we extrapolated the varying distribution of the cancer initiation patterns ( Figure 2B and Supplementary Figure 1 ). Lower fitness of Type 1 cells, r ₁ generally favored Type S-1 cells dominance when fitness of Type S-1, r_{S -1}, is high. Higher r ₁ values favored Type 1 dominance while equal fitness of Type 1 and Type S-1 cells yielded Type 1 dominance or Type 1/S-1 co-dominance. Mutation rates generally affected time to cancer detection and appearance of dominance. We also extended the mutation rates from Type 1 to Type S-1 cell type to denote other additional mutation steps and found a consistent increase in cancer detection times but patterns generally remained the same. In some cases where low fitness of both Type 1 and Type S-1 were coupled with lower mutation rates, Type S malignant cells failed to appear at extended times and simulations were terminated. Most curiously, mutation rate from Type 0 to Type 1, μ ₁, did not affect the pattern of cancer initiation or time to cancer detection ( Supplementary Figure 1 ).

Next, we examined the time to recurrence after surgical resection and the proportion of Type S-1 cells at the time of surgery in varying parameter sets. We reasoned that since Type S-1 cells needs only one more step for malignant transformation; therefore, its proportion was thought to be critical for cancer recurrence. To do this, we ran 1,000 simulations for each parameter set and calculated the mean recurrence time ( Figures 3A–K ). We also ran similar simulations at higher cell number in a tissue, N, between 100 to 1,000 times to assess the effect of tissue size on the parameter dependency ( Figures 3B–L ). We found that higher fitness of Type 1 cells, r ₁ increased the mean recurrence time ( Figures 3A, B ), while mutation rate from Type 0 to Type 1, μ ₁, had no effect on mean recurrence time ( Figures 3G, H ). Other parameters however, showed a negative correlation to the mean recurrence time – higher parameter values resulted in shorter mean recurrence time. Higher tissue cell number yielded an overall shortening of mean recurrence time but parameter dependency remained the same. We also observed a reduction in the proportion of Type S-1 cells at the time of surgery when r ₁, r_S and μ_S increases ( Figures 3A–L ), while r_{S -1} and μ_S increases in the proportion of Type S-1 cells ( Figures 3C–J ). Mutation rate from Type 0 to Type 1, μ ₁, had little effect on the proportion of Type S-1 cells ( Figures 3G, H ).

We then incorporated the spatial structure framework into the model to investigate the effect of tissue positional influence in cancer initiation patterns and time of cancer detection. After triplicate simulations using various parameters, we identified seven distinct patterns of cancer initiation ( Figure 4A ) based on the composition of non-malignant cell population at cancer detection. Figure 4B and Supplementary Figure 2 showed the distribution of the patterns in a wide parameter range. Low r ₁ values showed Type 0 dominance at low r_{S -1} levels with failure to detect cancer cells at very low μ_{S -1} levels. With higher r ₁ values, Type 1 cell types begin to dominate. When combined with high μ_{S -1}, we saw Type 1/S-1 co-dominance (green region in Figure 2B ). When r ₁ and r_{S -1} are equal to fitness of Type 0 (r ₀), we saw Type 0/1 co-dominance or Type 0/1/S-1 co-dominance depending on μ_{S -1} values. We noticed a peculiar pattern of Type 0/S-1 co-dominance (black region in Figure 2B ) when r_{S -1}, μ_{S -1}, and μ_S were high with relatively lower r ₁ value. Type S-1 dominance (red region in Figure 2B ) was regarded as the most undesirable scenario due to the abundance of Type S-1 cells, indicating shorter recurrence time. We saw this pattern when r_{S -1}, μ ₁ and μ _S-1 were high, r ₁ was equal to 1.0 and μ _S was relatively small. Some parameter sets with low fitness failed to yield Type S cells at extended time points during the simulations. Here, the incorporation of the spatial structure to our simulation framework had remarkable alterations to the cancer initiation patterns and cancer detection time.

Subsequently, we examined the mean recurrence time after surgical resection and the proportion of Type S-1 lesions at the time of surgery in a vast parameter range with the influence of the spatial structure setting ( Figure 5 ). Similarly, we ran 100 to 500 simulations to obtain mean recurrence time ( Figures 5A–K ) and to check the effect of larger cell numbers ( Figures 5B–L ). Generally, we saw that the integration of the spatial structure to our simulation framework had noteworthy changes to the parameter dependency to recurrence time. When the size of the normal tissue was small, the effect of fitness advantage on the proportion of Type S-1 cells in a tissue became larger. Our simulation results showed that an increase in the cell fitness shortened the mean time to recurrence ( Figures 5C–F ). However, an increase in r ₁ was found to reduce the recurrence time but begin to increase slightly at much higher levels regardless of the tissue size. We also found a consistent reduction in the mean recurrence time as mutation rates μ ₁, μ_{S -1} and μ_S increased ( Figures 5G–L ). We also observed a reduction in the proportion of Type S-1 cells at the time of surgery with a spatial structure. Especially, when either r_S or μ_S was small, and any of r_{S- 1}, μ ₁, or μ_{S- 1}was large, the proportion of Type S-1 increased ( Figure 5 ).

By using our computational model with spatial structure, multiple runs of stochastic simulations were performed with multiple parameter sets and in silico Kaplan–Meier curves were made. The data points about the time when 0% to 100% of patients experienced recurrence with an interval of 4% (time when 0%, 4%, 8%, …, 100% of patients experienced recurrence) were employed to compare between the in silico and published clinical data (44, 45) of 27 cancer types. In this analysis, we adopted random sampling for parameters to obtain in silico recurrence data and determined the best parameter set for each cancer type that minimized the mean of squared logarithmic residuals (log-MSR) between outputs in silico and in public. The accepted parameter set ( Table 1 ) was used to extrapolate recurrence time which were then fitted to clinical data and disease-free survival curves were depicted ( Figure 6 ).

Mantel-Cox test was used to compare between the curves of simulated and clinical data revealing minimal statistical nonconformity. According to the estimated parameters ( Table 1 ), we firstly deduced a tissue-specific turnover per month from d_S . Kidney chromophobe had the fewest cellular turnover cycles per month while bladder urothelial carcinoma and colorectal adenocarcinoma had the highest turnover cycles. Moreover, colorectal adenocarcinoma, kidney chromophobe, renal clear cell carcinoma, thyroid carcinoma, adenoid cystic carcinoma and acral melanoma showed higher fitness of all their premalignant cells than normal cells. Of note, the proliferation rate of the Type S malignant cells, r_S , was estimated to be high in cholangiocarcinoma, liver hepatocellular carcinoma, mesothelioma and upper tract urothelial cancer while being relatively low in breast invasive carcinoma, kidney chromophobe and skin cutaneous melanoma. Kidney chromophobe had the lowest mutation rate from the final premalignant cell stage to malignant cells while cervical squamous cell carcinoma and prostate adenocarcinoma had the highest mutation rate. Figure 7 showed the negative correlation between mutational steps required for carcinogenesis (37) and overall mutation rates (μ_I ) obtained from our studies by multiplying the mutation rates for all steps.