Modeling breast cancer progression to bone: how driver mutation order and metabolism matter

Abstract

Background

Not all the mutations are equally important for the development of metastasis. What about their order? The survival of cancer cells from the primary tumour site to the secondary seeding sites depends on the occurrence of very few driver mutations promoting oncogenic cell behaviours. Usually these driver mutations are among the most effective clinically actionable target markers. The quantitative evaluation of the effects of a mutation across primary and secondary sites is an important challenging problem that can lead to better predictability of cancer progression trajectory.

Results

We introduce a quantitative model in the framework of Cellular Automata to investigate the effects of metabolic mutations and mutation order on cancer stemness and tumour cell migration from breast, blood to bone metastasised sites. Our approach models three types of mutations: driver, the order of which is relevant for the dynamics, metabolic which support cancer growth and are estimated from existing databases, and non–driver mutations. We integrate the model with bioinformatics analysis on a cancer mutation database that shows metabolism-modifying alterations constitute an important class of key cancer mutations.

Conclusions

Our work provides a quantitative basis of how the order of driver mutations and the number of mutations altering metabolic processis matter for different cancer clones through their progression in breast, blood and bone compartments. This work is innovative because of multi compartment analysis and could impact proliferation of therapy-resistant clonal populations and patient survival. Mathematical modelling of the order of mutations is presented in terms of operators in an accessible way to the broad community of researchers in cancer models so to inspire further developments of this useful (and underused in biomedical models) methodology. We believe our results and the theoretical framework could also suggest experiments to measure the overall personalised cancer mutational signature.

Background

The development of tumours and the formation of metastasis are caused by cells which due to mutations are incapable of correctly sensing the external signalling of the environment and the surrounding cells so to adjust their metabolic and mitotic cycle activities. The purpose of exchange of signals between the cell and the environment is limiting their proliferation which, in the case of development of the cancer disease, becomes uncontrolled. There are different ways the single cells and organisms can defend against the development of tumour, for example, DNA self-repair mechanisms, signals inducing apoptosis in case of mechanical and biochemical instabilities of a cell and a prompt immune response which is capable of recognizing cancer cells. Mutations are accumulated by cells during their life and passed as inheritance along different progenies. In absence of external causes of mutations such as the exposition to radiations or chemicals, errors in the DNA replication and crossing over during the mitotic cycle are the major sources of mutations. Nowak proposed a model of cancer originating from propagation and accumulation of mutations occurring since early age [1]. Nevertheless, cancer is an ageing related disease as its development before the age of 30 years is rare.

Mutations may generate disruptions of cell cycle proliferation checkpoints, generate cytoskeleton instabilities inducing cell death or metabolic pathway deregulations increasing the cell proliferation, variation in the antigens and receptors on the membrane which transduce external signals (for example avoiding the immune response). These mechanisms produce the spreading of the sub–populations of cells (clones) with accumulated mutations. The populations of cancer cells in the primary and secondary sites are characterized by a large number of genes with altered gene expression due to mutations. The tumour heterogeneity results in genotypic and phenotypic (for example drug resistance) variance in the neoplastic cell populations. The heterogeneity increases during time as a consequence of new mutations acquired by cells during subsequent proliferations. Comparing gene expressions among cells in the primary site with those in the metastatic sites, one can observe the heterogeneities between different sites are not identical in variance, as well as, in the type of mutations accumulated; furthermore, looking at overlapping subsets of mutations between different tumour sites, it is possible to derive the time a secondary site originated, and if it is due to an offspring from the primary site or a nearby secondary metastasis. If we consider cancer as an evolving disease in competition for oxygen and energetic resources with other healthy cells, mutational heterogeneity represents a strategy initiated by an early small number of cancer stem cells. Therefore, the increase of heterogeneity can be seen as a diffusion process with a drift in the multidimensional space of the cell state belonging to {0,1}^N where N is the dimension of the space and represents the maximum number of genes affected by the disease during all its evolution. We believe that in order to relate cancer evolution with patient’s survival we need to take into account the characteristics of cancer stem cells, the classes of mutations and for some classes, also the order of mutations.

The work is structured in the following way. In the next subsections, we discuss the role of cancer stemness, and we define the type of mutations modelled and their effects on cells. In the “Model limitations” section, we introduce the concept of order of driver mutations, and we present the corresponding mathematical formulation. After which, we describe the set of rules driving the model dynamics from which we derive the master equations in the physical time. We model the effects of metabolic mutations on the cell cycle in terms of waiting time distributions and compute the final form of the master equation depending on the transition rates. The definition of the functional form of the transition rates in terms of the cancer stemness follows. Further discussion on the order of mutations in terms of ladder operators and the mathematical derivation of the effective driver mutations is addressed in the last method subsection. In the “Results” section, we present how simulations are carried out and the analysis of data supporting both the metabolic and driver mutations followed by the discussion and comparison of the three cases of interest numerically simulated.

The role of Cancer Stemness

Stem cells are capable of both self-renewing and differentiating [2]; this means they preserve themselves during proliferation without undergoing extinction due to differentiation, and they are a source for more committed cells [3]. The process of cell differentiation is mainly caused by epigenetic changes, and it results in the appearance of new cell phenotypes. These changes in the cell state are induced by external signalling or by internal variations of the cell dynamics like methylation or segregation of factors during mitosis. Not all the signals and changes involved in the differentiations are persistent or permanent. The loss of the new acquired phenotype is called de-differentiation. Nevertheless, the restoration of the external niche preserving the stemness or the circulation of factors inducing the cell stem state might not suffice to re-establish the stem condition in differentiated cells or in cells proliferating in a stem-like favourable condition [4]. Therefore, differentiated cells tendentiously do not de-differentiate.

The renewal condition is met when a cell will always undergo asymmetric division or undifferentiated symmetric proliferation. Stem cells are considered renewal at the level of single cells, meaning after proliferating a stem cell produces at least a daughter equal to itself. On the other hand, the hypotheses advanced by Knoblich and Morrison in [5, 6] imply that stem cells might symmetrically proliferate differentiated cells when the pool of (neighbour) stem cells do not risk extinction or has reached the maximum number the niche can sustain in a stem state. This last statement makes the definition of stem cell very close to the concept of population of cells which is in a state of high fitness given that such state is the one with more undifferentiated characteristics when compared with its progenies. The stemness of a cell is a behaviour which is not only present or absent, but there are different degrees of stemness actuation. A cell can lose its stemness when it cannot produce a more differentiated progeny or it cannot indefinitely proliferate; vice versa, a cell can increase its stemness by de-differentiation.

In an elevated mutational rate condition, changes in the DNA interfere or compete with the cell differentiation process, and it may become the predominant source of stemness variations [3]. Cancer stem cells are capable of producing a heterogeneous population which presents different phenotypes. The mutations will change the internal state of the cell affecting the transcription and translation processes. Cancer cells which will become more unstable and less fit will disappear. Cancer cells with a high fitness phenotype will survive if they proliferate asymmetrically or they do not acquire too many degrading mutations. Cancer stemness is the behaviour of cancer cells to preserve their phenotype and generate a heterogeneous population, part of which can show more stem-like behaviours [7].

Classes of mutations

Cancer mutations have been distinguished in drivers and non–drivers (or passengers). The accumulation of evidences of clonal heterogeneity and the observation of the arise of drug resistance in clonal sub–populations, suggest that mutations usually classified as non–drivers may have an important role in the fitness of the cancer cell and in the evolution and physiopathology of cancer. Similarly, mutations that alter the metabolism may modify the fitness of the cancer cell. The evolution of cancer clones in the context of the driver and non–driver mutations has been visualized as a braided river, with the capacity to diverge and converge [8].

Recently 122 potential immune response drivers–genetic regions in which mutations correlate with the presence or absence of immune cells infiltrating the tumours, have been reported. Current research into how tumours hide focuses on the over-expression of so-called checkpoint inhibitors in cancer cells [9, 10].

There are multiple levels of heterogeneity: inter-patient, intra-patient and intra-tumour. Inter-patient heterogeneity exists at the level where even within cancer of the same type, patients exhibit differences in terms of both biology and prognosis. In breast cancer, based on gene expression, patients can be classified into at least four broad subtypes, namely basal (triple negative breast cancer (TNBC)), HER2+, and luminal A and B subtypes. Luminal A and B are estrogen positive cancers, with luminal A having the best prognosis; HER2+ has overexpression of the HER2 growth-enhancing gene. These are slow growing tumours that respond well to treatment. The Basal or TNBC subtype is triple negative for estrogen, progestin receptors and HER2. This is the most aggressive subtype and is unresponsive to treatment. Extending on this, intra-patient heterogeneity is defined by the differences in the primary and the metastatic tumour sites; these include morphological and genetic differences, and differences in terms of tumour aggressiveness and proliferation. Intratumour heterogeneity is then explained by differences between regions within the same tumour mass [11–17].

The probing of the mutational space represents an advantage that cancer cell clonal populations have adopted. Also, the temporal aspect of mutations is of great importance in the evolution of the disease [18, 19].

In this work, we group the mutations in three main categories: driver mutations, metabolic mutations and non–driver mutations.

The driver mutations are mutations of specific genes responsible for the deregulation of pathways involved in the production of cytomembrane proteins which induce new pro-oncogenic behaviours; they increase the capability of cancer cells to self-renewal. Therefore, the driver mutations increase the fitness of cancer cells in a tissue and their capability of creating a progeny carrying all the driver and non–driver mutations accumulated. Hence, according to the definition in [20], driver mutations increase the cancer cell stemness, even though it is not possible to explicitly say a cell with a higher number of driver mutations corresponds to less uncommitted cell phenotype as in the case of healthy stem cells. A relevant aspect in the evolution of the disease is that the most important driver mutationss are very few. Works from A. Trumpp [21] show a number between 4 and 5 drivers mutations may be responsible for high correlation with survival. Recent works have provided evidence, at least in one cancer type, the order of the driver mutations affects the cancer clone trajectories and the patient survival [18, 19].

The role of non–driver mutations in cancer development and evolution has been modelled in [22]. The authors consider two cases: 1) the neutral non–driver mutations which do not affect the survival of cells and 2) the deleterious non–driver mutations as opposing to the proliferative effects resulting from the pro-oncogenic driver mutations, and consequently, causing a possible melting down of the cancer cell sub–populations. The non–driver mutations are those mutations which do not promote the formation of cancer stem cells and do not provide new oncogenic behaviours for the cells. Non–driver mutations generally affect the genes which are not related to the pathways promoting proliferation, migration, resistance to immune response and the formation of colonies. Cancer cells acquiring non–driver mutations become more unstable both in the cytoskeletric structure and in their fitness.

The metabolic mutations include all the set of mutations relative to the cell energetic needs and the controls of cell proliferation [23]. We consider that cells with accumulated metabolic mutations have an accelerated proliferation cycle because there is a sustained availability of ATP, a reduced sensitivity to cell activity inhibiting signals, and many of the checking procedures which are present during the various proliferation phases are not activated. Consequently, a higher number of these mutations further speed up the cell proliferation by reducing the time necessary to activate all these processes and the costs in terms of resources and energy.

Metabolic mutations differentiate from driver and non-driver mutations in that they do not determine how often the mitotic cycle of a cell is initiated or how frequent an apoptotic signal is triggered in a population of cells, but they affect how fast the clock of a cell runs. If each cell acquires one mutation at a time, then cells with few mutations are younger than cells with a large number of mutations. Metabolic mutations increase the mixing of groups of cells with different ages. In a heterogeneous population of cancer cells, at any physical time, there are cells that have different numbers of total mutations and have completed their mitotic cycle a different number of times based on their internal clock independently from their fitness. A high number of metabolic mutations is advantageous for those cells with a high fitness, but it has a negative ageing effect on populations with negative selective pressure.

Having considered the evolution of cancer cells in the primary site into circulating tumour cells, we have considered driver mutations inducing apoptotic or proliferative behaviours via transduction signalling. Instead, mutations responsible for activating or deactivating cellular processes prevalently involved with the internal dynamics of the cells have been relegated to non–driver and metabolic mutations. This choice takes into account the known fact that mutations affecting the cellular dynamics strongly related with the cytoplasm and nuclear domains, especially at the initial stage of cancer development, will dictate the fate of a cell and the occurrence of consecutive mutations which are less directly linked with extracellular domain interactions and the migratory process.

Results

Modelling the order of driver mutations

Although the effects of mutations on cancer heterogeneity and phenotypic plasticity of cancer evolution has led to interesting pure theoretical models (see for instance [24]), and it has also been addressed by combining theoretical modelling, bioinformatics and experimental data (see for instance [25–28] among others), there is a paucity of models that take into account the effects of the order of mutations. In [29], mutations are accumulated in a predefined sequence, and all the cells, sooner or later, will follow the same trajectory if they are not eliminated, or do not die along the way. In other cases, the gene expression of a sample of cancer cells are analysed at different times [30, 31], or in various sites [32] in order to identify the order in which the mutations occurred [33–35].

The fate of a cell is strongly related to the environment it is surrounded by. In our case, we consider that the order of mutations is given only by the environment, and it is not related to the dynamic, or evolution of the cancer cells. Indeed, the mutation order does not represent a tight preferential order of occurrence of mutations. On the contrary, the mutation order does not induce any constraint or bias on the occurrence of mutations which is instead random. The driver mutations affect the fitness of a tumoural mutated cell in a given compartment.

Considering mutation processes for which the weight of further mutations depends on their order of occurrence would require following the complete mutation history of a cell; therefore, in order to limit the set of values a driver mutation can induce on the state and fitness of a cell, we considered only phenotypical driver mutations for which the selection pressure depends uniquely on the environment.

In this work, the order among driver mutations is given by the maximization of survival of cancer cells which travel through a sequence of compartments from primary tumour site to secondary cancer site. We define a path as the sequence of tissues a hypothetical cancer cell will visit during its migrations before forming a metastasis. Different cancer cells can traverse different organs or tissues defining different paths. All the paths begin at the primary tumour site, and each of them ends in a site where cancer cells can form a secondary colony. Nevertheless, due to multiple metastasization sites, the final point of the paths can differ. Generally, a path has a different probability of being traversed than others. This means that among the total amount of cells forming the primary tumour only a portion will travel through a specific path, and some of the paths will be more common than others.

Let us define the set of m compartments $\mathcal {C}=\{c_{1},c_{2},\ldots, c_{m}\}$ where, for convenience, c₁ represents the primary site and the last s compartments {c_m−s+1,…,c_m} are the secondary sites. Let us also introduce the transition probability T_ij from the compartment c_i to another compartment c_j for those cells which already have all possible useful mutations to pass through all the compartments. The transition matrix T is such that $\sum _{j} T_{i j}=1$, and all the elements on the diagonal are one for all the secondary sites and zero otherwise. For each compartment c_j, we define the cell density ρ_j(n) after n transitions. If σ(i) is a permutation of the m compartments $\mathcal {C}$ visited by the cells, then we say a path of length l is given by

$$\begin{array}{@{}rcl@{}} \mathcal{P}_{\mathcal{C}}(\sigma)=\{\sigma(1),\sigma(2),\ldots,\sigma(l)\} \end{array} $$

and, after n≤l transitions, the population which has followed the $\mathcal {P}_{\mathcal {C}}(\sigma)$ path is:

$$\begin{array}{@{}rcl@{}} \rho_{\sigma(n)}=\prod_{i=1}^{n-1}T_{\sigma(i),\sigma(i+1)}\rho_{\sigma(1)}. \end{array} $$

If at the beginning all the cells are in c₁ so that ρ_j(n=0)=δ_j,1, then we must impose that a path’ s initial point is σ(1)=1. Furthermore, a path ends when it reaches for the first time one of the absorption points meaning that, after l steps, c_σ(l) corresponds to one of the secondary sites (σ(l)>m−s). Many paths will have a final population density ρ(l)=0 because the transition T_ij=0 before reaching the end of the path. If $\mathcal {P}_{\mathcal {C}}$ is the set of paths obtained by all possible permutations of the compartments, then $\overline {\mathcal {P}_{\mathcal {C}}}=\max _{\rho (l)}(\mathcal {P}_{\mathcal {C}})$ are the most probable paths. For the sake of simplicity, from here on, we restrict ourselves to the case $\overline {\mathcal {P}_{\mathcal {C}}}$ contains only one element representing the path with higher probability of being traversed. Nevertheless, there are no extra difficulties in considering multiple paths one at a time, and then studying the jointly resultant dynamics. Hence, the order we consider is given by the order of the tissues (compartments) the cancer cell will visit during its migration from the primary tumour site towards the secondary metastatic sites. The reason for introducing this type of order is related to the fact that a cancer cell migrating into a different compartment is a rare event which is almost impossible if the cell does not develop the right behaviours due to the occurrence of proper mutations.

A similar explanation is valid in the case of the survival (fitness landscape) of a cell while it remains for a period of time in a compartment. In different tissues, the cells need different membrane proteins to separate from the other surrounding cells, to migrate, to access or exit a specific compartment and to be recognized by the immune system as self. All these behaviours are not necessary at the same time, and in all the compartments, but only a few of them at a time are required during the travel through a compartment or in the transitions between compartments. As a consequence of the concomitant facts that a tumour cell needs specific driver mutations and that mutations can not be counterbalanced, the originating site and preferential way for a cell to reach the secondary sites determine the right order of accumulation of driver mutations. Therefore, if we order the compartments in the sequential order the cancer cells will traverse them, then we can establish the proper order of driver mutations for a high fitness and a short permanence during the development of metastases.

After we have specified the path $\mathcal {P}$, let us introduce the index $k \in {\mathbb {N}}$ identifying the compartments in the order given by the path. In the case the cells need only one driver mutation during each migratory transition between two compartments the index k can also be associated to the required driver mutations necessary to continue the travel along the path. If d_k is the driver mutation needed to go from compartment k to compartment k+1, then the cell which follows the specific path defined as an increasing value of k at each migratory transition till compartment $\overline {k}=l$, will also generate the sequence D={d₁,d₂,…,d_l} which defines the right mutation order.

We can extend the previous case to situations where the number of driver mutations for each compartment is more than one. Hence, if $\mathcal {C}=\{c_{1},\ldots,c_{l}\}$ is the sequence of compartment in the order visited by the cancer cells, and, for any k, S_k is the set of driver mutations necessary in the compartment c_k to improve the fitness landscape or to increase the probability of migrating in a different compartment, then we can define two sequences of sets as follow:

$$\begin{array}{*{20}l} S^{\prime}_{k} &=&\bigcup_{j=1}^{k} S_{j}, \qquad\quad 1\leq k\leq l \end{array} $$

(1)

$$\begin{array}{*{20}l} S^{\prime\prime}_{k}&=&S^{\prime}_{k}\setminus S_{k-1}, \qquad\quad 1\leq k\leq l \end{array} $$

(2)

$$\begin{array}{*{20}l} S^{\prime\prime}_{0}&=&\emptyset. \end{array} $$

(3)

The k-th element in the sequence $\mathcal {S^{\prime }}=\{S^{\prime }_{1},\ldots,S^{\prime }_{l}\}$ in Eq. 1 is the cumulated mutations of cells following the right order of mutations which have just entered in the c_k+1 compartment, while the k-th element in the sequence $\mathcal {S^{\prime \prime }}=\{S^{\prime \prime }_{1},\ldots,S^{\prime \prime }_{l}\}$ in Eq. 2 describes all the mutations for a cell following the right order of mutations that has just entered in the c_k compartment necessary in order to pass in the k+1-th compartment. If S_i∩S_j=∅ for any 1≤i,j≤l, then $S^{\prime \prime }_{k}=S_{k}$, but in general mutations in a set $\mathcal {S}=\{S_{1},\ldots,S_{l}\}$ can be necessary in other compartments. The right order of mutations, hence, also represents the most fit and the fastest path in the mutation space. Any deviation from the right order of mutations reduces to an increased cancer cells apoptosis rate and a longer period of time for cells in a given compartment.

The order of mutations determines many aspects of the life of a cancer stem cell. In general, it is safe to say that the mutational selection is a filter dynamically influenced by the order of appearance of driver mutations. During the loss of stemness and acquisition of different phenotypes, a cancer cell experiences different environments, and in each of them, it has a different value of fitness depending on the presence of specific mutations. These particular mutations may have already occurred in the history of the cancer clone so to assure the survival. We believe that the fitness of a cancer cell is a complex function that is modified by environmental factors such as immune system, drugs and therapies. The introduction of the order of driver mutations does not completely constraint the evolution of the system; indeed, in terms of single cells, the presence of this order does not at all affect the appearance of a mutation and is absolutely not perceived by the cells meaning that cells cannot use the right order of mutations to make any present choice or forecast any future state of the system. Therefore, single cells will continue to acquire mutations randomly. The cells and their progenies will explore the mutation space following different trajectories not limited by a unique mutation order. Eventually, the trajectory will feel the effects of the mutation order in terms of dynamics.

The interactions between a cell and the surrounding environment is what determines the dynamics of the cell itself [2]. As approximation, in average, a compartment is uniform and ready to interact with those cells which present the right membrane proteins. On the other hand, the cells, most of the time, are neutral to the environment interactions similar to a neutral molecule in an electric field.

We consider that driver mutations become effective only if they are preceded by other specific driver mutations. In other words, we introduce a unique sequential order of the driver mutations. If a driver mutation is acquired in a different order from the one predefined, it remains ineffective until all the driver mutations, which should have been occurred before according to the sequential order, are activated.

Data analysis supporting metabolism-related mutations

The Cosmic (http://cancer.sanger.ac.uk/) is the world’s largest and most comprehensive database resource for exploring the impact of somatic mutations. Other valuable databases include The Gene Expression Omnibus (GEO, www.ncpi.nlm.nih.gov/gds) among others. We have developed a R program that performs statistical analysis on the Cosmic database, identifies and evaluates the functional impact of each mutation by using a combination of pathways and gene ontology approaches. The functional annotation and enrichment analysis allows to classify large lists of genes into metabolism related groups. We have then extracted useful statistical estimators that allows to quantify the role of the metabolism-related mutations in the context of driver and non–driver mutations and the heterogeneity of proliferation rates among the cells. We have analysed the mutations with different threshold of FATHMM-MKL value which measures the impact of a mutation. The FATHMM-MKL algorithm predicts the functional, molecular and phenotypic consequences of protein missense variants using hidden Markov models [36].

This analysis shows the importance of the metabolic mutations in the context of the mutation impact calculated with the FATHMM score (Table 1).

Table 1 Go Ontologies (from the top biological processes, molecular functions, cellular components) of the mutations with the highest FATHMM score >0.9; the last section shows the relative enrichment of the mutations with the highest FATHMM with respect to a background of low FATHMM score below 0.7

Numerical results: cases study

To highlight the role of effective driver mutations and the role of the order in which the driver mutations occur, we look at the differences with the unordered driver mutations dynamics in which cancer cells evolve and migrate to different tissue depending only on the number of driver mutations m_d and not on the specificity of the corresponding genes. We simulate the two distinct cases and compare the recovering capabilities of cancer cell populations after a drug targeting sub–populations of cells within a specific range of the state σ={d,m_n,m_m,k} has been administered in one compartment. The effect of the drug is simulated by killing all the cells which match a given signature at the time of administration t_drug. Therefore, let us define a drug target signature σ_drug={d_drug,m_{n drug},m_{m drug},k_drug}, the amount of driver mutations m_{d drug} and the vector of driver mutations d_drug belonging to the same space of the cell driver mutations $\{0,1\}^{\overline {\mathrm {m_{d}}}}$, so that all the cells immediately sent to apoptosis are those for which hold true both the relations:

$$\begin{array}{*{20}l} \sigma_{\text{drug}}\ \circ\ \sigma=\sigma_{\text{drug}}\ \circ\ \sigma_{\text{drug}}, \end{array} $$

(4a)

$$\begin{array}{*{20}l} {m}_{\mathrm{d} \, \text{drug}} ||d||_{1}={{m}_{\mathrm{d} \, \text{drug}}}^{2}{\kern17pt} \end{array} $$

(4b)

where ∘ is the Hadamard product and ||·||₁ is the taxi cab norm of a vector. Then, we observe the recovery of the populations after the perturbation induced by the drug affecting the target cells is applied at the same time in the ordered and unordered simulations. In the simulations, we have considered the biological significant case of breast cancer cells metastasising in the bone. The process encompasses three different tissues. The mammary ducts are the primary site where cells develop malignant mutations and become tumourigenic. The second tissue is the circulatory system. Cancer cells which undergo epithelial to mesenchymal transition will initiate to migrate and, chemoattracted to regions with higher concentration of oxygen and nutrients, eventually they will intravasate through the proximal blood vessels. There, the cancer cells, also named Circulating Tumour Cells (CTC), will travel in the circulatory system. The CTCs which are recognized as self by the immune system will survive longer in the circulatory system, but only the cells capable of extravasating while reach the secondary sites. The third tissue is the bone. Breast cancer cells’ preferred metastasisation sites are regions enriched in TGF– β, and bone tissue has one of the highest concentration of TGF– β. Bone tissue represents the statistically preferred secondary site for breast cancer cells. Nonetheless, not all cancer cells reaching the secondary site will be enough tumourigenic to generate metastases, but only those which undergo mesenchymal to epithelial transition show sufficient cancer stem-like behaviours with the potentiality to form a secondary colony. Cell survival and capability of migrating is due to specific membrane proteins involved in breast cancer metastasising in the bone. A small subset of four genes corresponding to membrane proteins which positively correlate with the evolution of the disease and the survival of the patients has been experimentally found in [29] and are EPCAM, CD47, CD44 and MET. Mutations of these genes correspond to driver mutations. The EPCAM is related to the cells’ connectivity of the mammary duct lumen and plays a role in the epithelial to mesenchymal transition [37]. The CD47 is responsible for the production of membrane proteins increasing the survival of cancer cells by helping them to evade macrophage phagocytosis. CD44 variant isoforms are highly expressed in carcinomas of epithelial origin and relate to tumour progression and metastatic potential of some cancers. The protein CD44 bound with hyaluronic acid and sugar rich coating molecules is found on the surface of endothelial cells and tumour cells. Different amounts of CD44 changes the arrangement of the coating sugars causing the exposition of the embedded adhesion mediated selectin and integrins favouring cell extravasation [38]. The MET gene regulates mesenchymal to epithelial transition, and its overexpression positively correlates with metastasis formation. The four driver mutations present a sequential order given by their respective functionality in the three tissues, hence the order of the genes and their respective mutations is {EPCAM, CD47, CD44, MET}. Using the definition in Eq. 2, we derive Table 2 which puts in relation each compartment with a subset of mutations necessary to reach or seed in the next compartment, hence, defining the right order of driver mutations. In the circulatory system, CTCs need two driver mutations to survive and reach the bone tissue. The order of mutations between CD47 and CD44 is not important and both are considered effective drivers. The last row of Table 2 shows the minimum and maximum number of effective driver mutations a cell can have in each compartment.

Table 2 Ordered subsets of driver mutations

We have simulated three distinct cases of drug targets. The first two cases are less realistic, but are proposed to show the effect of drugs acting on unordered driver mutation systems, hence targeting cells with non specificity of the mutated gene. While, the last case is a realistic example of a drug affecting cells in a specific compartment with specific driver mutations. The drug target signature for the cases studied are:

All the components of the drug signature equal to zero are disregarded. Consequently, cells can have any value corresponding to each of the null components and still be a drug target.

In case 1 of Table 3, all the circulating tumour cells which have acquired a single driver mutation are killed by the drug at time t_drug as shown in Fig. 1. On the one hand, this means in the simulations with no ordered driver mutations, the cells targeted in the circulatory system could have acquired any of the four driver mutations listed in Table 2. Considering that when the order of mutations does not affect the dynamics, the cells which have one of any of the driver mutations can already intravasate, and, by entering in the circulatory system, they allow the restoration of the number of CTCs with one driver constituting a reservoir for cells which later acquire a second and a third driver mutation by asymmetrical division. Furthermore, the drug does not target sub–populations of cells having two or more driver mutations. Indeed, the dynamics of the sub–populations of cells having more than one driver when the drug is administered are close to the simulations with no drugs meaning they are mainly reservoirs of their own sub–populations due to symmetrical proliferation.

Fig. 1

Plot of the number of cell populations in function of the time. The 3 tissues are shown in the columns ordered from left to right following the cell traversing. All the combinations of ordered (w/ order) and unordered (w/o order) mutation dynamics together with drug (w/ drug) and without drug (w/o drug) administration are considered and shown in rows. Each curve represents the sub–population of cells with a specific number of driver mutations: m_d=0 is red, m_d=1 is green, m_d=2 is pink, m_d=3 is blue, and m_d=4 is yellow. The event of the drug administration is at time t_drug which is marked with a vertical line for easier comparison. The target of the drug, Eqs. 4a and 4b, and the compartment affected is shown at the bottom

Table 3 Numerical values for different types of treatments

On the other hand, in simulations with ordered driver mutations, a drug targeting any cell with one driver (case 1) is more effective when compared with the simulations with no order of driver mutations only in terms of total number of CTCs. The reason for less CTCs is that the minimum requirement to intravasate is having the mutation of the EPCAM gene which is more restrictive when compared to the unordered driver dynamics. Indeed, the sub–population of cancer cells with two driver mutations is close to the corresponding sub–population in the simulations with no order. The sub–population of cells exiting the mummary duct and entering the circulatory system is represented in both ordered and unordered driver simulations by the sub–population with one driver mutation referring to the EPCAM gene in the former and to a generic driver in the latter. Case 1 highlights a set of situations where perturbing sub–populations of cells in a compartment with few driver mutations compared to the total number of mutations necessary to change compartment, hence which are upstream in respect to the local progenies, does not affect sub–populations with more driver mutations, and therefore, the role of asymmetric proliferation is less relevant due to the loss of cancer stemness.

In case 2 of Table 3, the drug targets cancer cells in the second compartment with two driver mutations. Cell sub–populations are plotted in Fig. 2. As in the previous case, the results present different dynamics and drug conditions in all three compartments. The effect of the treatment is equal in both simulations with order and without order of driver mutation respectively. Indeed, the sub–population of circulating tumour cells with two driver mutations are going to zero at time t_drug. The differences between the two simulations is given by the dynamics of recovering of the targeted sub–populations and by the competition with the non targeted ones.

Fig. 2

Plot of the number of cell populations in function of the time. The 3 tissues are shown in the columns ordered from left to right following the cell traversing. All the combinations of ordered (w/ order) and unordered (w/o order) mutation dynamics together with drug (w/ drug) and without drug (w/o drug) administration are considered and shown in rows. Each curve represents the sub–population of cells with a specific number of driver mutations: m_d=0 is red, m_d=1 is green, m_d=2 is pink, m_d=3 is blue, and m_d=4 is yellow. The event of the drug administration is at time t_drug which is marked with a vertical line for easier comparison. The target of the drug, Eqs.4a and 4b, and the compartment affected is shown at the bottom

Similarly to case 1, when no treatment is administered, the number of CTCs with one driver mutation is larger in the simulations with unordered dynamics due to a less restrictive filter on the transition rates of cancer cells to pass from the mammary duct into the circulatory system than in the simulations with ordered dynamics. On the contrary, sub–populations with two driver mutations is larger in the ordered dynamics than in the unordered one because there are larger reservoirs of cancer cells in the mammary duct which have a null transition rate as a consequence of the wrong order of driver mutations. These reservoirs are composed of heterogeneous populations of cancer cells which go from committed phenotypes to more stem-like phenotypes. When the necessary mutations are acquired (in this specific case the required mutation is on the EPCAM gene) by any of the cells in the reservoirs, their transition rate becomes different from zero, and they contribute to the sub–populations with higher numbers of driver mutations shown in the last row of Table 2. Cells with more stemness will have more chances to survive and divide asymmetrically, while committed cells will have less possibility to contribute to the evolution of the disease. Similarly, many CTCs are obliged to remain in circulatory system until they have all the required mutations necessary for the extravasation. Hence, in the simulations with order of driver mutations, there are CTCs with more than two driver mutations, while in the simulations without order, there is only a sub–population of cells with two driver mutations which are already ready to extravasate.

The evolution of sub–populations of cancer cells reaching the bone compartment are similar in both ordered and unordered types of simulations when the system is not perturbed at all. Instead, when a perturbation is applied at time t_drug, the simulations with unordered dynamics show no sub–population of cells capable of seeding and metastasise in the bone, while, in the simulations with order of driver mutation dynamics, the number of aggressive cancer cells capable of reaching the bone tissue is even larger than in the absence of perturbation. This major difference between the two dynamics stems from the fact that in simulations with ordered driver mutations there are CTCs having high value of cancer-stemness which are missing only few (in this case just the CD47 and CD44 genes) driver mutations, but are constrained to remain in the circulatory system. The lack of extravasating capabilities make these sub–populations of cancer cells apparently less aggressive, when instead they are a few mutational steps away from the aggressive phenotype. Furthermore, part of these cells acquired a large number of metabolic mutations which contribute to the speed of the cell cycle and their effective capability of forming secondary colonies.

Case 3 of Table 3 differentiates from the former cases for the specificity of the drug that targets only circulating tumour cells with the CD47 gene mutated, see Fig. 3. In the dynamics with ordered driver mutations, the sub–populations affected are those with two and three driver mutations, but the treatment leaves untouched cells with only one driver mutation referring to the EPCAM gene. In this case, we can see the corresponding curves do not go to zero as in the former cases because the respective sub–populations contain cells without the CD47 mutation. In the unordered dynamics, the drug also hit the sub–population of cells with one driver mutation, and it is more effective compared with the ordered dynamics. Precisely in terms of combinatorics, the drug in the unordered dynamics targets 1/4 of all the combinations of cells with one driver mutations, 1/2 of the combinations of cells with two driver mutations and 3/4 of the combinations of cells with three driver mutations. In the dynamics with ordered driver mutations, the drug targets none of the cells with one driver mutation, 1/3 of the combinations of cells with two drivers and 2/3 of the combinations of cells with three driver mutations.

Fig. 3

Plot of the number of cell populations in function of the time. The 3 tissues are shown in the columns ordered from left to right following the cell traversing. All the combinations of ordered (w/ order) and unordered (w/o order) mutation dynamics together with drug (w/ drug) and without drug (w/o drug) administration are considered and shown in rows. Each curve represents the sub–population of cells with a specific number of driver mutations: m_d=0 is red, m_d=1 is green, m_d=2 is pink, m_d=3 is blue, and m_d=4 is yellow. The event of the drug administration is at time t_drug which is marked with a vertical line for easier comparison. The target of the drug, Eqs. 4a and 4b, and the compartment affected is shown at the bottom

Discussions

The order of driver mutations introduces a strong constraint in the migration of cancer cells from one compartment to another compared to a cancer evolution driven by unordered dynamics, and the former limits the cell heterogeneity much more than the latter. Even though the ordered dynamics is slower as consequence of the limits imposed on the transition rates by driver mutations waiting to become effective driver mutations, the ordered and unordered dynamics are comparable up to opportune rescaling of parameters. Nevertheless, the two dynamics are characterized by completely different types of evolution of the disease after the same kind of drug treatment. Simulations based on the unordered driver mutations show a slower recovery of the cancer populations and a retardation in the appearance of highly tumorigenic sub–populations, while simulations with ordered driver mutations show a faster restoration of the targeted cells and the anticipation of the appearance of aggressive sub–populations.

The perturbation induced by the drug produces different results on the two dynamics by breaking the symmetries between the heterogeneous cancer cell sub–populations. In the examples shown on breast cancer metastasizing in the bone, the restoration induced by the ordered dynamics strongly affects the extravasating cancer cells because, although the cells with mutation in the CD47 gene are sent to apoptosis by the drug, many other cells which present much less effective driver mutations (in this specific case only one effective driver mutation) survive, and some of these surviving cells are mutationally a few steps away from becoming highly tumorigenic.

Hence, drugs targeting the larger sub–populations of cells just entering the compartment or the sub–populations surviving longer when not treated are less effective in the long time regime because they eliminate the cells with a high number of effective driver mutations. However, they leave cells with an even higher number of driver mutations with a non aggressive phenotype, but which are missing a few mutations to unlock the complete sequence of the right order of mutations and switch to an aggressive tumorigenic phenotype. Furthermore, cells with few effective driver mutations can posses high stem-like potentiality due to acquired driver mutations which are not effective. These cells, in the long time regime, represent a source for many other cancer cells speeding up the restoration of cancer sub–populations.

Conclusions

We investigated the role of heterogeneity in cancer cells by embedding a new internal structure driving the development of the disease in the dynamics of the mutation process. Such internal structure is given by the order of occurrence of mutations and the metabolic cell cycle acceleration. The structure introduced in the dynamics stems from a phenomenological derivation of the effects of mutations which results in an increase/decrease of the capability of cancer cells to survive, differentiate, proliferate and metastasise. On one side, we have considered driver mutations responsible for the increase of cancer cell stemness and their order relevant for the development of oncogenic phenotypes involved in the process of metastasization in a secondary site of the cell and its progenies. On the other side, metabolic mutations are related to the production and consumption of ATP and are involved in the acceleration of the mitotic rate.

The order of mutations and acceleration of the dynamics have been considered relevant in some specific types of cancers like melanomas, but their importance and effects have been experimentally difficult to measure. To understand their effects on cancer evolution, we introduced an analytical description of the order of driver mutations and the effects of metabolic mutations in terms of operators, and we have included them in the derivation of the master equations for the cancer cell populations.

This model has been applied to the case of breast cancer metastasising in the bone tissue for which we simulated the evolution of cancer in presence and absence of the order of driver mutations. To highlight the differences between the two types of dynamics, we have compared them before and after they are affected by an external perturbation represented by the action of an anti-tumoral drug. The numerical results show the order of mutations introduce a slower dynamics of cancer cells reaching the bone than in simulations with no order of driver mutations. Nevertheless, for realistic drug perturbations of the two dynamics, the recovery of aggressive cancer sub-populations is faster in presence of order of driver mutations than in the case of cancer cells non sensitive to the mutation order.

This model pinpoints how the order of mutations and metabolic mutations in cancer may be responsible for the short time of relapse after drug treatments. Furthermore, this work helps us to understand the difficulties in experimentally detecting the order of mutation in cancer, their role in the dynamics and their response due to the interactions of the system with drugs.

There is a growing number of algorithms and tools to detect driver mutations, see for instance [39]. The model proposed in this work is based on a large number of evidences on clonal evolutionary cancer trajectories depending on mutation order. We cite here for instance [40–43] among others. Note that in [42], the effect of drugs (chemotherapy) is found important for the selection on the order.

The details of our theoretical model have been implemented in a software developed in Mathematica based on an Cellular Automata framework which is highly adaptable and allows us to simulate different hypothesis and evolution models under different numerical parameters some of which as been shown in this work.

Model limitations

In our model, metabolic mutations and driver mutations are two disjoint sets of mutations. This decision simplifies the derivation of the master equation driving the dynamics of cell populations. In some specific cases, separate classes of mutations may represent an unrealistic constraint for the underlying biological process; nonetheless, such limitations can be easily solved by adding specific transition and branching terms for the occurrence of each mutation which at the same time, is metabolic and driver. We proposed a model which is Markovian even though it includes the order of driver mutations. Markovianity is maintained by considering the driver mutations effective or non effective due to their order. In cases where the effect of each driver mutation on the fitness of the cell is more complex than just a stepwise binary state driven by the order, the Markovian behaviour is lost, and the model complexity increases with the number of mutations requiring suitable numerical oriented approaches. The number of distinct cell populations is described by the cell state, and in our final analyses, we grouped them by the number of effective driver mutations to derive our results. The type and dimensionality of the mutational space represent a constraint on the variability of the cell state and the processes which could be described. A mutational space which discriminates only if specific mutations occurred is a simplification to the detriment of heterogeneity.

Methods

Cancer dynamics and evolution

In order to mimic the tumour evolution, each cancer cell will evolve following four types of actions depending on the internal state of the cell and the environment constraints. The actions are asymmetric proliferation, symmetric proliferation, apoptosis and migration. These actions can be expressed in reaction form as follow:

$$\begin{array}{*{20}l} C_{\{\boldsymbol{d},m_{\mathrm{m}},m_{\mathrm{n}},k\}}&\rightarrow C_{\{{\boldsymbol{d}}+ {\boldsymbol{u}}_{{\boldsymbol{i}}},m_{\mathrm{m}} + 1, m_{\mathrm{n}},k\}} + C_{\{{\boldsymbol{d}}, m_{\mathrm{m}} + 1, m_{\mathrm{n}},k\}}, \end{array} $$

(5a)

$$\begin{array}{*{20}l} C_{\{\boldsymbol{d},m_{\mathrm{m}},m_{\mathrm{n}},k\}}&\rightarrow C_{\{\boldsymbol{d}, m_{\mathrm{m}} + 1, m_{\mathrm{n}} + 1, k\}} + C_{\{\boldsymbol{d}, m_{\mathrm{m}} + 1, m_{\mathrm{n}},k\}}, \end{array} $$

(5b)

$$\begin{array}{*{20}l} C_{\{\boldsymbol{d},m_{\mathrm{m}},m_{\mathrm{n}},k\}}&\rightarrow C_{\{\boldsymbol{d}+ \boldsymbol{u}_{\boldsymbol{i}},m_{\mathrm{m}} + 1, m_{\mathrm{n}},k\}}\,+\, C_{\{\boldsymbol{d}+ \boldsymbol{u}_{\boldsymbol{i}},m_{\mathrm{m}} + 1, m_{\mathrm{n}},k\}}, \end{array} $$

(6a)

$$\begin{array}{*{20}l} C_{\{\boldsymbol{d},m_{\mathrm{m}},m_{\mathrm{n}},k\}}&\rightarrow C_{\{\boldsymbol{d}, m_{\mathrm{m}} + 1, m_{\mathrm{n}} + 1, k\}} + C_{\{\boldsymbol{d}, m_{\mathrm{m}} + 1, m_{\mathrm{n}} + 1, k\}}, \end{array} $$

(6b)

$$ C_{\{\boldsymbol{d},m_{\mathrm{m}},m_{\mathrm{n}},k\}}\rightarrow 0, $$

(7)

$$ C_{\{\boldsymbol{d},m_{\mathrm{m}},m_{\mathrm{n}},k\}} \rightarrow C_{\{\boldsymbol{d},m_{\mathrm{m}},m_{\mathrm{n}},k + 1\}}, $$

(8)

where $\phantom {\dot {i}\!}C_{\{\boldsymbol {d}, m_{\mathrm {m}}, m_{\mathrm {n}},k\}}$ is the cell having d specific driver mutations, an amount of m_m metabolic mutations, m_n non–driver mutations and being located inside the compartment k. The reactions Eqs. 5a, 5b, 6a, 6b, 7 and 8 indicate that the cells C on the left side of → undergo a process which will end with their annihilation and the creation of new cells specified by the right side of →. The states of the cells are expressed between the brackets {…}.

The vector d has a dimension $\overline {m_{\mathrm {d}}}$ and belongs to the space $\{0,1\}^{\overline {m_{\mathrm {d}}}}$. Each component d_j of the vector d takes into account the state of the single driver mutation $j \in [1, \overline {m_{\mathrm {d}}}]$, and d_j=1 if the mutation is present or d_j=0 otherwise. The set of vectors u_i have all the components equal to zero except the i-th component which is equal to one. The number of the acquired and unordered metabolic and non–driver mutations are represented by the integer numbers m_m and m_n, respectively. In the model proposed, we impose that each mutation can occur at most one time independently on the type and severity of the mutation. This persistent accumulation effect on the mutation space is more explicit for the driver mutations. In fact, the initial state of the driver mutations for any cells is given by the vector d=0 having all the components equal to zero, and after a driver mutation occurs, the corresponding component become equal to one. In the same way, a further mutation relative to the same gene is excluded, meaning the component will never switch back to zero, and that gene will never be chosen again among the possible mutating genes, hence reducing the mutation space. It is worth to remark the state of the system changes only when there is a mutation corresponding to a single gene, and its dynamic is sensible only to such mutations. Despite the previous statements, the model can be easily extended to include a set of genes for each driver mutation. These sets do not need to be disjoint, but they can share some common genes, and when a mutation for a specific gene occurs, its state (or counting) is changed in each set where the gene is included, see Additional file 1. Adding the above model extension will allow us to trace more precisely the diffusion in the metabolic mutation space caused by each driver mutation and their temporal order.

The reactions Eqs. 5a and 5b describe the proliferation of cells principally characterized by stem-like behaviours, hence following an asymmetric mytotic cycle. The two reactions represent the cases when a driver and a non–driver mutation occurs, respectively. The reactions Eqs. 6a and 6b show the symmetric proliferation when a driver and a non–driver mutation occurs, respectively. The reactions Eqs. 6a and 6b are related to more differentiated cells which show progenitor-like characteristics. The reaction Eq. 7 represents cells undergoing apoptosis. Generally this is the fate followed by cells with high instabilities or committed cells with low stem-like characteristics, which have a reduced proliferation activities even when stimulated. The last reaction, Eq. 8, takes into account the migration of cells which have developed mutations allowing them to move and pass from an environment to another. The migration process depends on the occurrence of proper driver mutations which need to be developed in the proper order given by the order of the environments, k→k+1.

We have assumed driver and non–driver mutations are two disjoint classes of mutations, but both can be included in the wider class of the metabolic mutations. Hence, in Eqs. 5a, 5b, 6a and 6b at each occurrence of a driver or non–driver mutation, there is an increase of one of the number of metabolic mutations. The reasons for this choice are that, on one side, the metabolic mutations mark the advancing age of the cell in terms of how many mutations steps it has acquired (see below for connections with stemness and commitment), and on the other side, the advantages and disadvantages produced by the mutations can be considered independent from the speed of the cell metabolic activity. In order to simplify the complex representation of the underlying biological system, we do not consider reactions which result in the increasing of only metabolic mutations.

We have six reactions when indeed there are only four distinct actions occurring in the system. The cause for these over number of reactions is related to the choice of representing the system as if there where multiple classes of mutations, and each mutation can be definitely tagged as belonging to one and only one of these classes. However, from a biological and experimental point of view, probably, not all the mutations involved can be marked with a distinct colour, but more like if they were a mix of shades. Indeed, a gene (and the molecules derived from its transcription and translation) can have a function in a compartment or at a given state and play another role somewhere else. For the sake of simplicity, we have grouped all the non–driver mutations in an integer variable with the purpose of counting them and neglecting the order of their occurrence. First, as we said above, non–driver mutations are not strictly necessary for the worsening and the progress of the disease, but are mutations bringing small disadvantages to the cells when compared with the large advantages of the driver mutations. Second, the order of occurrence of the non-driver mutations in respect to all the mutations is less important than for the driver mutations, because they are weakly linked to the driver mutations and to the environment fitness constraints as well as compartment transition rates.

The rates of the reactions Eqs. 5a, 5b, 6a, 6b, 7 and 8 are not constant, but depend on the cancer stemness. Therefore, each cell has specific probabilities of proliferating or undergoing apoptosis. If r∈[1,4]⊂N is the index addressing one of the four reactions defined above, then we can define the probability function P_r(s) that the reaction r occurs as follow:

$$\begin{array}{*{20}l} &P_{r}: s\in \mathit{S}\rightarrow [0,1], \end{array} $$

(9)

$$\begin{array}{*{20}l} &\!\!\!\!\!\!\!\!\sum_{r} P_{r}(s)=1 \forall s, \end{array} $$

(10)

where s is a real positive value in [0,1] and it represents the cancer stemness of the cell. The stemness can be generated by a random process.

The cancer stemness does not depend explicitly on the gene expression, but only on the mutation state. Therefore, on the complete mutational space, the cells proceed along trajectories of cancer stemness which increases or decreases due to acquisition of new driver and new non–driver mutations, while the stemmness does not depend directly on the number of metabolic mutations. Eventually, if cells survive, they will be attracted toward the asymptotic value given by s(max(m_e), max(m_n)). The amount of cells per compartment is limited by a carrying capacity which determines the maximum number of cells that can be stored in a given compartment. The maximum number of cells depends on the oxygen/energetic needs and it changes in function of the type of the metabolic activities. The cumulative energetic needs are computed as:

$$\begin{array}{@{}rcl@{}} E=\sum_{j} e(j) N(j) \end{array} $$

(11)

where j is the number of driver mutations accumulated, e is the energetic needs depending on the accumulated driver mutations and N is the total number of cells in a compartment with given j mutations.

Master equation of the mutation process

Let us consider the cell state σ={d,m_n,m_m,k} with d driver mutations, m_n passenger mutations, m_m metabolic mutations and in the compartment k. the index total m_m represents also the age of the trajectory in the natural time of the occurrence of the events m_m=d·d+m_n.

The rates of transition involved in each action depends on the cancer cell stemness which is linked to the internal state of the cell, σ={d,m_n,m_m,k}, and to the number of mutations, m_d and m_n. Also the probabilities of acquiring a driver or a passenger mutation depend on the number of mutations already accumulated by the cells.

The conditional probability density ρ(d,m_n,m_m) that a cell starting from the state σ₀={0,0,0,k} arrives after m_m steps to the state d,m_n,m_m,k is:

$${} {\begin{aligned} &\rho(\boldsymbol{d},m_{\mathrm{n}},m_{\mathrm{m}},k)=\\ &\left[{\frac{1}{\overline{m_{\mathrm{d}}}}\sum_{i=1}^{\overline{m_{\mathrm{d}}}} r_{\text{asym}} (m_{\mathrm{d}}\,-\,1, m_{\mathrm{n}},k) r_{\mathrm{d}}(m_{\mathrm{d}}-1,m_{\mathrm{n}},k) \rho(\boldsymbol{d}-\boldsymbol{u}_{\boldsymbol{i}},m_{\mathrm{n}}, m_{\mathrm{m}}-1,k)}\right.\\ &\left. + r_{\text{asym}}(m_{\mathrm{d}},m_{\mathrm{n}}-1,k) r_{\mathrm{n}} (m_{\mathrm{d}},m_{\mathrm{n}}-1,k) \, \rho(\boldsymbol{d},m_{\mathrm{n}}-1,m_{\mathrm{m}}-1,k){\vphantom{\sum_{i=1}^{\overline{m_{\mathrm{d}}}}}}\right]\\ &\!+2\left[\frac{1}{\overline{m_{\mathrm{d}}}}\sum_{i=1}^{\overline{m_{\mathrm{d}}}} r_{\text{sym}} (m_{\mathrm{d}}\,-\,1, m_{\mathrm{n}},k) r_{\mathrm{d}}(m_{\mathrm{d}}\,-\,1,m_{\mathrm{n}},k) \rho(\boldsymbol{d}\,-\,\boldsymbol{u}_{\boldsymbol{i}},m_{\mathrm{n}}, m_{\mathrm{m}}\,-\,1,k)\right.\\ &\left. + r_{\text{sym}}(m_{\mathrm{d}},m_{\mathrm{n}}-1,k) r_{\mathrm{n}} (m_{\mathrm{d}},m_{\mathrm{n}}-1,k) \, \rho(\boldsymbol{d},m_{\mathrm{n}}-1,m_{\mathrm{m}}-1,k){\vphantom{\sum_{i=1}^{\overline{m_{\mathrm{d}}}}}}\right]\\ &-r_{\text{sym}}(m_{\mathrm{d}},m_{\mathrm{n}},k) \, \rho(\boldsymbol{d}, m_{\mathrm{n}}, m_{\mathrm{m}},k)-r_{\text{apop}}(m_{\mathrm{d}}, m_{\mathrm{n}},k) \rho(\boldsymbol{d}, m_{\mathrm{n}}, m_{\mathrm{m}},k)\\ &+\!\left[r_{\text{pass}}(m_{\mathrm{d}}, m_{\mathrm{n}}, k\,-\,1)\rho(\boldsymbol{d}, m_{\mathrm{n}}, m_{\mathrm{m}}, k\,-\,1)\,-\, r_{\text{pass}}(m_{\mathrm{d}}, m_{\mathrm{n}},k)\rho(\boldsymbol{d}, m_{\mathrm{n}}, m_{\mathrm{m}},k)\!\right] \end{aligned}} $$

(12)

The first term enclosed in square brackets describes the increase of number of cells in the state σ={d,m_n,m_m,k} due to asymmetric proliferation. The mutation occurring during the asymmetric proliferation can be driver or passenger; hence r_d(σ)+r_n(σ)=1. The second term enclosed in square brackets takes into account the increase of cells due to symmetric proliferation, while the third term express the fact that both the daughter cells equally change their state and there is no self-renewal. The fourth term is the decreasing of cells due to apoptosis, and the last term is due to the change of compartment.

In order to compare the dynamics with the biological process, we need to switch from the natural time of the events to the physical time [44], meaning the time is measured with a macroscopic clock advancing with a regular periodic step. We do so by using the waiting time distribution ψ(i,t) that the i-th mutational event occurs exactly at time t.

To take into account the effect of acceleration due to metabolic mutations, we introduce an exponential waiting time which depends on the number of metabolic mutations m_m:

$${} \psi(m_{\mathrm{m}},t)=\alpha(r+m_{\mathrm{m}}) e^{-\alpha r t} e^{-\alpha m_{\mathrm{m}} t}, $$

(13)

where $\frac {1}{\alpha r}$ is the cell cycle mean time of cells with no mutations, and α is the increase of the cell cycle rate caused by each metabolic mutation.

After substituting Eq. 13 into Eq. 12 and applying some mathematical transformations to ρ(d,m_n,m_m,k) (see the Additional file 2 for further details), we obtain the final form of master equation for the probability density p(d,m_n,m_m,k,t) of finding a cell with σ={d,m_n,m_m,k} mutations at time t with any possible order of the occurrence of d:

$$ \begin{aligned} &\quad\,\,\,\,\,\frac{\mathrm{d}}{\mathrm{d}t} p(\boldsymbol{d}, m_{\mathrm{n}}, m_{\mathrm{m}}, k, t)= {\alpha (r+ m_{\mathrm{m}}-1)} \\ &\left\{\left[\frac{1}{\overline{m_{\mathrm{d}}}}\sum_{i=1}^{\overline{m_{\mathrm{d}}}} r_{\text{asym}}\, r_{\mathrm{d}} p(\boldsymbol{d}-\boldsymbol{u}_{\boldsymbol{i}}, m_{\mathrm{n}}, m_{\mathrm{m}}-1,k,t)\right.\right.\\ &\left.\left.\quad\,\,\,\,\,+ r_{\text{asym}} \, r_{\mathrm{n}}\, p(\boldsymbol{d}, m_{\mathrm{n}}, m_{\mathrm{m}}-1, k,t){\vphantom{\sum_{i=1}^{\overline{m_{\mathrm{d}}}}}}\right]\right.\\ &+2 \left[\frac{1}{\overline{m_{\mathrm{d}}}}\sum_{i=1}^{m_{\mathrm{m}}} r_{\text{sym}}\, r_{\mathrm{d}} p(\boldsymbol{d}-\boldsymbol{u}_{\boldsymbol{i}}, m_{\mathrm{n}}, m_{\mathrm{m}}-1,k,t)\right.\\ &\left.\left.\quad\,\,\,\,\,+ r_{\text{sym}} \, r_{\mathrm{n}}\, p(\boldsymbol{d}, m_{\mathrm{n}}, m_{\mathrm{m}}-1, k,t){\vphantom{\sum_{i=1}^{m_{\mathrm{m}}}}}\right]\right\} +{\alpha(r+m_{\mathrm{m}})}\\ &\quad\,\,\,\,\,\left\{ \,- r_{\text{sym}} \, p(\boldsymbol{d}, m_{\mathrm{n}}, m_{\mathrm{m}}, k,s)-r_{apop} \, p(\boldsymbol{d}, m_{\mathrm{n}}, m_{\mathrm{m}},k, t)\right.\\ &\left.+\left[r_{\text{pass}} \, p(\boldsymbol{d}, m_{\mathrm{n}}, m_{\mathrm{m}}, k-1, t)- r_{\text{pass}} \, p(\boldsymbol{d}, m_{\mathrm{n}}, m_{\mathrm{m}},k,t) \right]\right\}. \end{aligned} $$

(14)

Equations and range of parameters

To completely describe the evolution of cancer cells’ probability densities given in Eq. 14, we need to define the transition rates for each type of action which depend on the mutational state of the cells in such a way to mimic the effects of the cancer stemness and the order of driver mutations. Therefore, we introduce a sufficiently generic form valid for all 4 type of actions described in Eqs. 5a, 5b, 6a, 6b, 7 and 8 which includes an explicit dependency on the cancer stemness s and on the amount of effective driver mutations m_e:

$$\begin{array}{@{}rcl@{}} r_{\text{act}}=\frac{1}{N} \, A_{\text{act}} \, \chi_{\text{act}}(s) \, k_{\text{act}}(m_{\mathrm{e}}, k). \end{array} $$

(15)

The rates are given by the product of a constant rate amplitude A_act modulated by a support function χ_act(s) depending only on the stemness of a cell multiplied by a filter function k_act(m_e,k) which is a function of the number of effective driver mutations and the compartment k divide by the normalization factor N. Each of the terms in the rate are action type dependent. The amplitude A_act determines how much more recurrent is the specific action in a given compartment. The support function χ_act∈[0,1] has the purpose of determining which actions are active for a specific value of the stemness cs. The filter k_act∈{0,1} selects the subgroup of cell states which present specific properties to perform an action. For each cell state, the normalization factor is defined as:

$$\begin{array}{@{}rcl@{}} N=\sum_{\text{act}} r_{\text{act}}. \end{array} $$

Because the cancer stemness plays a role in the survival and proliferation capability of the cells and is less important in the transition from one tissue to another, we consider χ_pass=1. Vice versa the effective driver mutations are mainly relevant for the migration of cells from one compartment to another, hence k_asym=k_sym=k_apop=1. For the remaining modulating functions, if the actions are allowed only for one contiguous compact interval of the cancer stemness s, then a valid choice for the support is:

$$\chi_{\text{act}}(s)=\Theta(s-w/2+c_{\text{act}}) \, \Theta(w/2-s-+c_{\text{act}}), $$

where Θ is a step function, and w is the extension of the support centred around c_act.

The former condition can be relaxed and more phenomenological derived functional forms, which do not present discontinuities and are derivable, can be chosen as support functions. Nevertheless, the product of Heaviside functions makes the problem simpler and allow us to describe a large variety of relevant cases.

To obtain the effect of asymmetric division for more stem-like cells and apoptotic tendency for more committed cells, then:

$$\begin{array}{@{}rcl@{}} c_{\text{asym}}\,< c_{\text{sym}}\,< c_{\text{apop}}<\overline{s}. \end{array} $$

These non-holonomic constraints, do not imply the support functions belong to disjoint range of the cancer stemness domain. Indeed, they can overlap so to mimic more realistic biological cases where for the same value of cancer stemness corresponds to multiple possible actions. On the other hand, for disjoint support functions, the cancer stemness strictly dictates the cell behaviour.

The filter function k_pass(m_e,k) is a simple step function which is zero if the cell does not have all the necessary driver mutations S_k defined in Eq. 2, and it is 1 if m_e≥S_k. It is worth to remark that effective driver mutations can be derived from the vector of driver mutations d by applying ladder operators similar to those used in quantum mechanics for describing the energetic state of a harmonic oscillator the order of driver mutations (for more details see Additional file 1).

The cancer stemness is a positive real value function in $[0,\overline {s}]$ defined on the sub-space of driver and non–driver mutations. As discussed in “The role of Cancer Stemness” section, the cancer stemness is related to the population of cancer cells due to their fitness. While the driver mutations (and the respective gene expressions) are responsible for the production of proteins and receptors on the surface of the cell membrane, the non–driver mutations are responsible for the cell instabilities related to the cytoskeleton. Hence a possible choice for the cancer stemness is $s(m_{\mathrm {d}}, m_{\mathrm {n}})=\overline {s}\frac {m_{\mathrm {d}}}{\sqrt {1+ m_{\mathrm {n}}}}$. The surface produced by s is a strictly monotonically increasing function of the driver mutations m_d and a strictly monotonically decreasing function of the non–driver mutation m_n, while it does not depend on the metabolic mutations m_m. The variation of the cancer stemness for each cell is a random process driven by the acquisition of mutations. If we regard the m_d and m_n as continuous parameters is it possible to define isocurves of cancer stemness. Cells whose trajectories remain on the same isocurve while acquiring mutations will not change their (average) behavior, while cells crossing the isocurves will become more stem–like if going toward higher levels or more committed if going toward lower levels.

Simulations

The simulations are done in the framework of CA where each unit represents a cancer cell having a position σ={d,m_n,m_m,k} in the cell state space composed by the mutation space and the tissue space. A simulation consists in a set of cells beginning with specific initial conditions and evolving together with the advancing of time. At each time, each cell continues to occupy the same position until an event occurs after which the corresponding cell moves to a new position in the state space. The evolution of cells, the creation of new others and their annihilation are given by stochastic events following the rules in Eqs. 5a, 5b, 6a, 6b, 7 and 8. There are two types of constraints for the cells. One is due to the limited extension of the cell state space given by the ranges $\phantom {\dot {i}\!}\{0,1\}^{m_{\mathrm {d}}}$, $\{0,\dots, \overline {m_{\mathrm {n}}}\}$, $\{0,\dots,\overline {m_{\mathrm {m}}}\}$ and $\{1,\dots,\overline {k}\}$ which are already included in the master equation. The second constraint is a volume exclusion limit which impose a maximum capacity of cells having the same number of driver mutations. The update of the system during the simulations is based on the asynchronous self-clocked time scheme where each cell has an independent timer. When a cell undergoes a specific action all the resulting cells from the product of the reaction set their timers to a random period given by the waiting time distribution Eq. 13 which depends on the cell state. Hence, the waiting times for the new events are chosen randomly after an event occurs and are based on the actual new cell state. This interval of time remains constant independently from other events occurring meanwhile. The type of action performed at the end of a time period is randomly chosen with probabilities given in Eq. 15 which depend only on the internal state of the cell. Because the state of a cell does not change between events, we are in a case of non competition between the reactions; therefore the type of action performed by a cell at the end of the time period can be tossed and decided at the beginning of the time period. At the end of the time period, the cell checks all the constraints and if the limits are not reached, the reaction is performed,otherwise it is put back. Simulations with the same initial conditions and parameters are repeated multiple times. At each time, the number of cells in the same cell state are averaged over the ensemble of simulations. Eventually, the mean trajectories and the deviations from the mean are computed. In some cases the trajectories of the cells are projected into a sub–space of the entire cell state and then averaged so to reduce the number of degrees of freedom and highlight the main characteristics of the systems corresponding to the observables which are experimentally detected. Our software is available upon request from the first author.