Cooperation among tumor cell subpopulations leads to intratumor heterogeneity

A recent compilation shows that there are about forty types of cancer based on the site of origin, such as lung, kidney, colon cancers and so on (see Table 1)[2]. Due to the coexistence of different cell types in some organs, cancers can be further classified as carcinoma (epithelial cells), sarcoma (connective tissue), myeloma (plasma cells), leukemia (bone marrow), and lymphoma (cells of the lymph system), blastoma (precursor cells) etc. If we take kidney cancer as an example, there are renal cell carcinoma, transitional cell carcinoma, nephroblastomas, and renal sarcoma. The renal cell carcinoma (RCC) is the most common kidney cancer, contributing to 90 percent of the cases, and it can be classified into a few subgroups (clear cell RCC, Papillary RCC, Chromophobe RCC, duct RCC and others) according to the phenotype of cancer cells[12, 13]. Indeed, it is quite astonishing that there are more than 100 different types of cancer known at present (see the full cancer list in $https://www.cancer.gov/types\#k$)[1, 14], which explains the difficulty in determining broad principles that drive the origin and evolution of these diseases. Yes, we are facing not a single disease but hundreds of them. Although all of these diseases share some similar hallmarks as summarized in the landmark reviews[1, 15], there appears to be no universal treatment (nor will there be in the near future) for these diseases due to the lack of understanding of the underlying mechanisms of cancer evolution.

In addition to the many faces of cancer that appear in different organs, another challenge for cancer treatment comes from tumor evolution. As proposed by Nowell in 1976, cancer is now widely viewed as a clonal evolutionary process[6, 11]. Due to the imperfect genome replication process, environment exposures and heredity, a normal cell can transform into a cancer cell by acquiring genetic mutations, which occur on the time scale of a few decades[16, 11]. During cancer progression, tumor cells constantly face the selective pressures derived from their complex microenvironment, such as competition from the surrounding healthy cells for nutrients and space for growth[17]. In addition, hypoxia, pH change, immune surveillance and other potential factors all threaten the survival and growth of tumor cells[18, 19, 20]. Therefore, tumor cells generate new traits through continuous evolution, which is just one of the major reasons for the failure of chemotherapy, radiotherapy and other widely applied methods used to treat cancer patients[6, 11].

We have gained considerable knowledge from studies of tumor evolution through animal models[21, 22]. The advent of many modern techniques, such as next generation sequencing have helped us get a deeper understanding for the genetic basis of cancer[23]. However, it is still very difficult to study tumor evolution in humans due to the inaccessibility of tumor biopsies from patients at different time points. Instead, the tumor evolutionary history obtained in most studies are derived from patient samples at a single time point based on assumptions, which are frequently violated as the tumor evolves[24]. Hence, a clear and complete picture about tumor evolution has not emerged. Many models have been proposed to describe the tumor evolution process, such as the linear sequential model, branched, neutral and punctuated evolution of tumors (see Fig. 1)[24, 25]. One common feature in all these models is that different subpopulations could appear and coexist in a tumor, which is the cause of pervasive cancer heterogeneity.

Cancer is driven by accumulation of genetic mutations especially the “driver mutations”, which convey selective advantage to cancer cells. With the advent of sequencing technology, genome-wide association studies (GWAS) is affordable, thus giving us a powerful tool to search for cancer driver genes. The list of cancer driver genes is continuously growing and the number has reached around 300 recently[26, 27]. A few cancer genes such as TP53 appear in many cancers while each cancer type usually has its own specific driver mutations. Therefore, different cancer patients have distinct tumor evolutionary processes, which leads to the interpatient heterogeneity (see Fig. 2)[7, 8]. Personalized cancer medicine has been proposed, and is necessary due to this type of heterogeneity[28].

As the cancer cells escape the primary site and seed other sites of the body, they finish the transition into the last and fatal stage, referred to as cancer metastasis - responsible for $90\%$ of cancer patient death[29]. Whether new driver mutations are required for the cancer metastasis is unresolved[30, 31, 32, 33]. However, distinct new mutations can appear both at the primary and metastatic sites after the spatial isolation is established among cells at different sites[34]. It leads to the next level of heterogeneity, intrapatient heterogeneity which is one of the reasons for cancer recurrence after treatment (see Fig. 2).

In addition to the two types of cancer heterogeneities described above, the intratumor heterogeneity (ITH, see Fig. 2), which refers to the coexistence of different subpopulations in a single tumor has been found in many cancers[9, 35, 36, 37, 38, 39, 40, 41]. ITH plays a crucial role in almost all stages of cancer such as tumor progression, metastasis, drug resistance and recurrence of cancer[42, 43, 44]. Therefore, understanding the ITH is a major step in investigating cancer heterogeneity because it can help design better therapy to avoid drug resistance and cancer recurrence.

The evolutionary models mentioned above all point to the possibility that different subpopulations (with distinct mutations) can coexist in a single tumor. In these models, the subpopulations interact with each other mainly through cell competition irrespective of the selective action under linear and branched models or the drift effect in neutral and punctuated models. However, increasing experimental evidence found that the cooperation among distinct cell subpopulations in a tumor is essential for tumor maintenance [45], enhanced tumor growth [46], and even cancer metastasis [44, 47, 48]. Surprisingly, a minor subpopulation is sufficient to support the whole tumor growth and determine the clinical course [49, 50, 51]. Hence, the different types of cell-cell interactions have to be considered to better understand the mechanisms of ITH, which is often neglected in many theoretical models. The lack of understanding of ITH greatly impairs the progress of developing more effective therapies for cancer patients.

In the past few decades, many innovative mathematical models have been proposed to study the cancer evolutionary process, including cancer initiation, progression, and metastases[52]. A simple multistage model that cancer is driven by a number of driver mutation events was proposed to explain the age-dependent cancer incidence rate more than sixty years ago[53]. The Moran process accounts for the cell division, apoptosis, and mutation process in tumor evolution for a fixed cell population size. It has been widely used to investigate the accumulation of mutations and cancer initiation[54, 55, 56, 57, 58]. Other models which further include the spatial structure of tumors and its related microenvironment such as the oxygen and nutrient concentrations are considered to investigate the growth dynamics of tumor by using partial differential equations or an agent-based model[59, 60, 61, 62, 63].

Evolutionary models have also been used to study ITH[64, 65, 66, 67], and many insightful results have been obtained from these studies, such as the extent of ITH and factors influencing ITH. There are two crucial assumptions in these theoretical studies. The genetic mutations only confer fitness advantage to the cell itself (cell-autonomous effect), and tumor cells interact with each other through competition. However, there is increasing evidence that tumors cannot overcome the microenvironmental constraints only through autonomous increase of the cell growth rate[68, 69, 70, 71]. In contrast, enhancement of survival and proliferation rates through non-cell-autonomous effects by factors such as metalloproteinases and cytokines are critical for tumor progression[72]. Therefore, these factors secreted by certain tumor/normal cells (‘producer’) can bring cooperation among different cell subpopulations instead of competition. Such a cooperation can in principle accelerate tumor progression. Instead of one population accumulating all the mutations required for cancer development[15], a few partially transformed subpopulations with each containing one or two mutations can realize this procedure through their cooperation[10, 73].

The evolutionary game theory (EGT) provides a novel and unique avenue for investigating ITH accounting for the interactions among subpopulaitons of tumor or between cancer and normal cells[74, 75, 76, 77, 78, 79]. The EGT is a subfield of game theory (GT), which provides mathematical models for studying the strategic interaction among individuals[80, 81]. An individual (‘player’) receives a payoff during the game depending on both the player and also the behavior (‘strategy’) of others. For EGT in cancer research, the players are cancer or normal cells, and the payoff is their fitness while the strategies are phenotypes adopted by players[77]. The advantage of EGT is that it can describe the time-dependent evolution of the relative abundance of each cell type, determine the equilibrium conditions and the stability of phenotype coexist[52]. Here, we briefly describe the EGT. The evolutionary games for two cells of different types $A$, and $B$ are frequently represented by the pay-off matrix,

where $W_{IJ}(I,J\equiv A,B)$ represents the fitness of the cell type $I$ interacting with the cell type $J$. The cell-cell interaction can be direct or indirect and its effect can also be competition or cooperation due to the influences of space, nutrients, information, growth factors and other microenvironmental factors[17, 77, 82]. Therefore, the fitness function $W_{IJ}$ of the cells can also be expressed by very different mathematical formula[77]. As the two different cell populations, instead of just two cells, are well-mixed with each other[82], the average fitness ($w_{A}$, $w_{B}$) of the two cell types may be described as,

$$w_{A}=f_{A}W_{AA}+(1-f_{A})W_{AB}\,,$$

(1)

$$w_{B}=f_{B}W_{BA}+(1-f_{B})W_{BB}\,,$$

(2)

where $f_{A}$, and $f_{B}$ are the fractions of the cell type $A$, and $B$ in the population, respectively. From the fitness functions of the cells, the time-(in)dependent properties of the cell population as mentioned above can be calculated.

In a few latest studies[83, 84], the EGT has been adopted to explain ITH in pancreatic cancer starting from the simplest case with only two different types of cancer cells. One cell type can produce a growth factor while the other does not. The growth factor (often called “public goods”) can be consumed by both cell types and promote their proliferation. Qualitative conclusions are obtained for the coexistence of two cell subpopulations which are consistent with experimental observations.

We consider two types of cancer cells with one of them producing a public good (see Fig. 3). In the insightful in vitro experiment[83], the public good is the insulin-like growth factor II (IGF-II). First, we focus on this experiment and then we apply the same idea to describe the ITH observed in an in vivo experiment on glioblastoma multiforme[91]. In the first experiment, both cell types are taken from cancer cells of mice with insulinomas (a type of pancreatic cancer). One cell type can produce IGF-II (the producer (+/+)) while the non-producer (-/-) cells have IGF-II gene deleted, which means they can no longer produce IGF-II but can still consume it from the environment.

One critical element in using the concept of the public goods game is how to define the payoff function. Often a rather complex, frequency dependent function is assigned to both players[86, 83]. Here, we use the growth rate as the fitness functions for both the cell types as they can be measured in experiments directly. A non-linear growth function[83] is found for the growth rate ($w_{-}$) of the non-producer cells as a function of IGF-II concentration ($c$), which is well-described by the Hill-like function,

$$w_{-}=a_{1}+\lambda_{1}c^{\alpha}/(a_{2}^{\alpha}+c^{\alpha})\,.$$

(3)

where $a_{1}=2.0$, $\lambda_{1}=18.9$, $\alpha=0.7$, and $a_{2}=3.2$[76]. Although it has not been measured, we expect that the producer cells also have a similar functional form as Eq. (3) for their fitness ($w_{+}$), which leads to,

$$w_{+}=g(c)-p_{0}\,,$$

(4)

where $g(c)$ is the same Hill-like function as in Eq. (3), and $p_{0}$ is the cost paid by producer cells due to the IGF-II production. The IGF-II concentration ($c$) has two sources once two cell types are mixed with each other. One is from the production of the producer cells, and the other is supplied exogenously (in the culture medium). The IGF-II produced by the producer will depend on the fraction $f_{+}$ in the population. Therefore, the IGF-II concentration ($c_{+}$) available for the producer can be written as,

$$c_{+}=af_{+}+c_{0}\,,$$

(5)

where $a$ is the coefficient for the allocation of IGF-II produced by +/+ cells, and $c_{0}$ represents the exogenous supply of IGF-II. Similarly, the available IGF-II ($c_{-}$) for the non-producer is given by,

$$c_{-}=bf_{+}+c_{0}\,,$$

(6)

The parameter $b$ is the coefficient of allocation of IGF-II produced by +/+ cells. We use different parameters $a,b$ to indicate that the producer and non-producer can get different fraction of IGF-II that are produced by the former cells.

The evolution of the fraction $f_{+}$ ($f_{-}$) of producers (non-producers) can be derived from the replicator equations,

$$\frac{\partial f_{+}}{\partial t}=(w_{+}-\langle w\rangle)f_{+}\,,$$

(7)

and,

$$\frac{\partial f_{-}}{\partial t}=(w_{-}-\langle w\rangle)f_{-}\,,$$

(8)

where the average fitness $\langle w\rangle$ is,

$$\langle w\rangle=w_{+}f_{+}+w_{-}f_{-}\,,$$

(9)

with $f_{+}+f_{-}=1$ being the normalization condition. Let the number of producers and non-producers be $N_{+}$, and $N_{-}$, respectively. The whole population size, $N$, is the sum of $N_{+}$ and $N_{-}$. The system size, $N$, is a time dependent quantity, which is often neglected in other models for simplicity [83, 92]. Here, we describe the time-dependent changes in the system size, $N$, through

$$\frac{\partial N}{\partial t}=w_{+}N_{+}+w_{-}N_{-}=\langle w\rangle N\,.$$

(10)

Then, we can study the conditions for the coexistence of producers and non-producers from the equations above.

The first important result derived using our model is that the emergence of a stable ITH requires an unequal allocation of public goods between producer and non-producer cells. The allocation of the public goods produced by the producers is determined by the ratio, $b/a$, in Eqs. (5) and (6). A unity for the value of $b/a$ means that the two cell types share the public goods equally. Therefore, the two cell types will have the same fitness function, expect a constant shift $p_{0}$ due to the cost paid by the producer (see Eqs. (3) and (4) and the Fig. 4A). There is only one stable state under this condition (see the filled blue dot in Figs. 4A and D), a homogeneous state consisted of only non-producer would be expected as long as the initial fraction $f_{+}(0)<1$. On the other hand, given $b/a=0$ which means producers do not share any public good with the non-producer. Then, the fitness of the non-producer does not depend on the fraction of the producer any more. There is a heterogeneous state but it is unstable under this condition (see the open circle in Fig. 4A and the evolution of $f_{+}$ in Fig. 4E). However, we find a stable heterogeneous state if $0<b/a<1$ (see the filled blue dots in Fig. 4C and F), which indicates that the producer maintains a higher fraction of the public goods produced. There is a negative feedback mechanism close to this stable coexisting state, which can be observed from Fig. 4C.

To validate our conclusion that it is the unequal allocation of public goods leading to the establishment of a stable coexist state, we calculated the the internal equilibrium fractions ($0<f_{+}^{i}<1$) of +/+ cells at different concentrations of serum ($c_{0}$) which has been detected in experiments. The +/+ cell fraction is measured after 5 days co-culture of +/+ and -/- cells under eleven different initial fractions of cells (with $f_{+}(t=0)$ almost evenly distributed between 0 and 1) and different amounts of serum[83]. There are three free parameters in our model, a, b, and $p_{0}$. From the equal fitness ($\approx$ 14.4) of producer and non-producer cells and the +/+ cell fraction approaching 1 at the stable internal state under $c_{0}=0$ observed in the experiments[83], the parameter b $\approx$ 8 is obtained from Eqs. (3) and (6). Similarly, we can derive the values of $a$ and $p_{0}$ based on the experimental results of internal equilibrium states. We find that our theoretical model captures the experimental observations very well (see the upper panel in Fig. 5). No quantitative explanations have been reached for similar experiments from other models at present. Interestingly, the ratio of $b/a$ obtained through the comparison between our theoretical results and experiments is 0.1 which is just located at the interval we proposed above for the formation of a stable ITH (see the low panel of Fig. 5 for the illustration of one stable ITH state).

Another important result is that the cost $p_{0}$ paid by producers has a strong influence on cell cooperation. As the fitness of producers is reduced by $p_{0}$ due to the production of the public goods, we investigate whether this parameter has any influence on the cooperation of different subpopulations. From the diagram of the internal equilibrium fraction $f_{+}^{i}$ of the producer as a function of the exogenous resource level at varied values of $p_{0}$ (see Fig. 6), we found a critical concentration for exogenous resources above which there is only one stable homogeneous state consisted of non-producers. And this critical concentration increases as the value of $p_{0}$ decreases. A second critical value (see the star symbols in Fig. 6A) is observed for exogenous resources as $p_{0}$ takes small values. There is no stable heterogeneous states after the exogenous resource level is lower than this critical value.

To better understand the influence of $p_{0}$ on cell cooperation, we calculated a phase diagram in terms of the initial fraction $f_{+}(0)$ and percentage of serum (see Figs. 6B and 6C). There are three distinct phases represented in different colors (pink, blue, and purple). Two stable homogeneous states consist of either producers (pink region) or non-producers (blue region). The third phase corresponds to a stable heterogeneous state consisted of both subpopulations (purple region). The area of the three regions is strongly influenced by the value of $p_{0}$ (see Figs. 6B and 6C). A heterogeneous phase is established at modest levels of exogenous resources and higher $f_{+}(0)$. We also find that a higher cost paid by the producers actually helps the establishment of a heterogeneous system (see the change of the purple region in Figs. 6B and 6C).

Cooperation among different cell types has been found in several types of cancer such as breast cancer[45], prostate cancer[93] and Glioblastoma multiforme (GBM)[91]. Here, we take GBM, which is the most aggressive and fatal brain cancer [94] as an example, to show how the same model based on public goods explains equally well the mechanism of ITH in GBM. Two types of cancer cells are frequently found in GBM with different expressions of epidermal growth factor receptor (EGFR) [95]. Apparently, the coexistence of the two cell types in GBM promotes the cancer progression and leads to a worse prognosis of the patients [96, 97]. Several studies found that the GBM cells with rearrangement of EGFR gene ($\Delta$ cells) can secrete factors (Interleukin-6 and Leukemia inhibitory factor) to enhance cell proliferation and inhibit its apoptosis [91, 98].

A recent study investigated the interactions between GBM cells with amplified levels of EGFR (referred to as WT cells) and $\Delta$ cells systematically[91]. A fixed amount of cancer cells ($2\times 10^{5}$ cells) with different initial fractions of the two cell types are injected into athymic nude mice (4 to 5 weeks old). Then, the size of the established tumors was measured at different times (see the different symbols in Fig. 7). It is difficult for the WT cells alone to induce a new tumor in the mice (see the pink upside down-triangles in Fig. 7) while a tumor can be quickly generated by the $\Delta$ cell alone (see the the blue squares in Fig. 7). Meanwhile, a faster tumor growth rate is found as the initial fraction of $\Delta$ cells increases until it reaches 90$\%$ (see the inset of Fig. 7). However, the maximum growth rate is found at the mixture of 10$\%$ WT and 90$\%$ $\Delta$ cells but not 100$\%$ $\Delta$ cells which seems to indicate a cooperative relation between these two cell types.

We applied the public goods model to explain the non-trivial experimental observations summarized above. The WT and $\Delta$ cells are considered to be the non-producer and producer, respectively. The IL-6 and/or LIF secreted by $\Delta$ cells is the public goods. We used a simple fitness function ($w_{-}$) for the WT cell as above which is given by,

$$w_{-}=bf_{+}/(1+bf_{+})\,.$$

(11)

The fitness function ($w_{+}$) for $\Delta$ cells is,

$$w_{+}=af_{+}/(1+af_{+})-p_{0}\,.$$

(12)

The different parameters $a$ and $b$ indicate that the public goods (produced by the $\Delta$ cells) can be shared unequally between the two cell types. There are only three free parameters in our model, which can be calculated in the following way. The average fitness of the population is given by $\langle w\rangle=w_{+}$ as it contains only producer cells with $f_{+}=1$ (see Eq. (9)). Therefore, the population size, $N$, can be calculated from Eq. (10) which increases exponentially, $N=N_{0}e^{w_{+}t}$. From the growth curve (see the blue square in the upper panel of Fig. 6) of the $\Delta$ cells measured in experiments, we obtain,

$$\frac{a}{1+a}-p_{0}\approx 0.335\,,$$

(13)

where 0.335 is obtained from the exponential fit. It is found that the two types of cells grow at the same rate given $f_{+}=f_{-}$[91], which leads to,

$$\frac{0.5a}{1+0.5a}-p_{0}=\frac{0.5b}{1+0.5b}\,.$$

(14)

In addition, the evolution of $N(t)$ can be described by $N(t)=N_{0}e^{w_{-}t}$ given $f_{+}=f_{-}$, and $w_{+}=w_{-}$ as the average fitness $\langle w\rangle=0.5(w_{+}+w_{-})=w_{-}$. Thus, we obtain the following expression,

$$\frac{0.5b}{1+0.5b}\approx 0.321\,,$$

(15)

from the tumor growth curve (orange dots) in the upper panel of Fig. 7. The constant 0.321 is obtained from the exponential fit (see the orange dotted line in Fig. 7). Therefore, the three parameters in our model can be derived using Eqs. (13)-(15), which leads to, $a=68.4$, $b=0.946$, and $p_{0}=0.651$.

To test our model, we predicted the tumor growth behavior at 10$\%$, and 90$\%$ of $\Delta$ cells, and found that our predictions agree very well with experimental observations (see the lower panel of Fig. 7). We also calculated the tumor growth rate as a function of the fraction of $\Delta$ cells[76] and found that there is a maximum tumor growth rate at an intermediate value ($\approx$0.77) of the $\Delta$ cells fraction which is consistent with experimental observations and can be test in future experiments.

Our theory also provides some hints for the frequent observations of both WT cells and $\Delta$ cells in GBM patient and its poor prognosis. From the evolution of the fraction of $\Delta$ cells under different initial conditions (Fig. 8A), we found that the heterogeneous state is stable even $f_{+}(0)$ varies from 0.1 to 0.9. Therefore, a quick recurrence of GBM is possible if some $\Delta$ cells still exist in the tumor. Cooperation between different cancer types greatly enhanced their survival and progression. Hence, new drugs and methods should be developed to eliminate such a relation. Our theory predicts that a rich culture medium can encourage cell competition instead of cooperation, which might provide potential treatment strategies for GBM and other similar cancer types (see Fig. 8B).