\begin{equation} W\left(GLY\right) = 1 + a * \left(f_1 + 1\right) \end{equation}

\begin{equation} W\left(VOP\right) = 1 + a * f_1 + v * \left(f_2 + 1\right) - c \end{equation}

\begin{equation} W\left(DEF\right) = 1 + a * f_1 + v * f_2 \end{equation}

Where $f_1$ is the GLY cells' frequency and $f_2$ is the VOF cells' frequency at a given time. All fitness equations start from balance by the sum of 1. Finally, depending of the parameter's values, the model can lead to three different situations (as in other examples, the different types are one mutation away from WT): ### Fully glycolytic tumours: If the fitness benefit of a single unit of acidification is higher than the maximum benefit from the club good for aerobic cells, then GLY cells will always have a strictly higher fitness than aerobic cells, and be selected for. In this scenario, the population will converge towards all GLY, regardless of the initial proportions (as long as there is at least some GLY in the population). ```{r glvop1, message=FALSE} # Definition of the function for creating the corresponding dataframe. avc <- function (a, v, c) { data.frame(Genotype = c("WT", "GLY", "VOP", "DEF"), Fitness = c("1", paste0("1 + ",a," * (f_GLY + 1)"), paste0("1 + ",a," * f_GLY + ",v," * (f_VOP + 1) - ",c), paste0("1 + ",a," * f_GLY + ",v," * f_VOP") )) } # Specification of the different effects on fitness. afavc <- allFitnessEffects(genotFitness = avc(2.5, 2, 1), frequencyDependentFitness = TRUE, frequencyType = "rel") ## For real, you would probably want to run ## this multiple times with oncoSimulPop simulation <- oncoSimulIndiv(afavc, model = "McFL", onlyCancer = FALSE, finalTime = 15, mu = 1e-3, initSize = 4000, keepPhylog = TRUE, seed = NULL, errorHitMaxTries = FALSE, errorHitWallTime = FALSE) ``` ```{r glvop2, message=FALSE} # Representation of the plot of one simulation as an example (the others are # highly similar). plot(simulation, show = "genotypes", type = "line", ylab = "Number of individuals", main = "Fully glycolytic tumours", font.main=2, font.lab=2, cex.main=1.4, cex.lab=1.1, las = 1) ``` ### Fully angiogenic tumours: If the benefit to VOP from their extra unit of vascularisation is higher than the cost c to produce that unit, then VOP will always have a strictly higher fitness than DEF, selecting the proportion of VOP cells towards 1. In addition, if the maximum possible benefit of the club good to aerobic cells is higher than the benefit of an extra unit of acidification, then for sufficiently high number of VOP, GLY will have lower fitness than aerobic cells. When both conditions are satisfied, the population will converge towards all VOP. ```{r glvop3, message=FALSE} # Definition of the function for creating the corresponding dataframe. avc <- function (a, v, c) { data.frame(Genotype = c("WT", "GLY", "VOP", "DEF"), Fitness = c("1", paste0("1 + ",a," * (f_GLY + 1)"), paste0("1 + ",a," * f_GLY + ",v, " * (f_VOP + 1) - ",c), paste0("1 + ",a," * f_GLY + ",v, " * f_VOP") )) } # Specification of the different effects on fitness. afavc <- allFitnessEffects(genotFitness = avc(2.5, 7, 1), frequencyDependentFitness = TRUE, frequencyType = "rel") simulation <- oncoSimulIndiv(afavc, model = "McFL", onlyCancer = FALSE, finalTime = 15, mu = 1e-4, initSize = 4000, keepPhylog = TRUE, seed = NULL, errorHitMaxTries = FALSE, errorHitWallTime = FALSE) ``` ```{r glvop4, message=FALSE} ## We get a huge number of VOP very quickly ## (too quickly?) plot(simulation, show = "genotypes", type = "line", ylab = "Number of individuals", main = "Fully angiogenic tumours", font.main=2, font.lab=2, cex.main=1.4, cex.lab=1.1, las = 1) ``` ### Heterogeneous tumours: If the benefit from an extra unit of vascularisation in a fully aerobic group is lower than the cost c to produce that unit, then for a sufficiently low proportion of GLY and thus sufficiently large number of aerobic cells sharing the club good, DEF will have higher fitness than VOP. This will lead to a decrease in the proportion of VOP among aerobic cells and thus a decrease in the average fitness of aerobic cells. A lower fitness in aerobic cells will lead to an increase in the proportion of GLY until the aerobic groups (among which the club good is split) get sufficiently small and fitness starts to favour VOP over DEF, swinging the dynamics back. ```{r glvop5, message=FALSE} # Definition of the function for creating the corresponding dataframe. avc <- function (a, v, c) { data.frame(Genotype = c("WT", "GLY", "VOP", "DEF"), Fitness = c("1", paste0("1 + ",a," * (f_GLY + 1)"), paste0("1 + ",a," * f_GLY + ",v," * (f_VOP + 1) - ",c), paste0("1 + ",a," * f_GLY + ",v," * f_VOP") )) } # Specification of the different effects on fitness. afavc <- allFitnessEffects(genotFitness = avc(7.5, 2, 1), frequencyDependentFitness = TRUE, frequencyType = "rel") # Launching of the simulation (20 times). simulation <- oncoSimulIndiv(afavc, model = "McFL", onlyCancer = FALSE, finalTime = 25, mu = 1e-4, initSize = 4000, keepPhylog = TRUE, seed = NULL, errorHitMaxTries = FALSE, errorHitWallTime = FALSE) ``` ```{r glvop6, message=FALSE} # Representation of the plot of one simulation as an example (the others are # highly similar). plot(simulation, show = "genotypes", type = "line", ylab = "Number of individuals", main = "Heterogeneous tumours", font.main=2, font.lab=2, cex.main=1.4, cex.lab=1.1, las = 1) ``` ## Prostate cancer tumour–stroma interactions {#prostatestroma} This example is based on @basanta_investigating_2012. The authors apply evolutionary game theory to model the behavior and progression of a prostate tumour formed by three different cell populations: stromal cells, a dependant tumour phenotype capable of co-opting stromal cells to support its growth and an independent tumour phenotype that does not require microenvironmental support, be it stromal associated or not. To enable this, the model has four variables, which is the minimun necessary to describe the relationships in terms of costs and benefits between the different types of cells and, of course, to describe the progression of the cancer. The different cell types, hence, are as follows: 1. S: stromal cells. 2. D: microenvironmental-dependent tumour cells. 3. I: microenvironmental-independent tumour cells. And the parameters that describe the relationships are as follows: * $\alpha$: benefit derived from the cooperation between a S cell and a D cell. * $\beta$: cost of extracting resources from the microenvironment. * $\gamma$: cost of being microenvironmentally independent. * $\rho$: benefit derived by D from paracrine growth factors produced by I cells. Table \@ref(tab:tableprostate) shows the payoffs for each cell type when interacting with others. We consider no other phenotypes are relevant in the context of the game and disregard spatial considerations. . | S | D | I - | ------------------ | ---------- | -------------------- **S** | $0$ | $\alpha$ | $0$ **D** | $1+\alpha-\beta$ | $1-2\beta$ | $1-\beta+\rho$ **I** | $1-\gamma$ | $1-\gamma$ | $1-\gamma$ : (\#tab:tableprostate) Payoff table that represents the interactions between the three cell types considered by de model As in @basanta_investigating_2012, the I cells are relatively independent from the microenvironment and produce their own growth factors (e.g. testosterone) and thus are considered to have a comparatively constant fitness $(1-\gamma)$, where $\gamma$ represents the fitness cost for I cells to be independent. The D cells rely more on their microenvironment for survival and growth at a fitness cost $(\beta)$ that represents the scarcity of resources or space that I cells can procure themselves. A resource-poor microenvironment would then be characterised by a higher value of $\beta$. As I cells produce space and shareable growth factors, this model assumes that D cells derive a fitness advantage from their interactions with I cells represented by the variable $\rho$. On the other hand, D cells interacting with other D cells will have a harder time sharing existing microenvironmental resources with other equally dependant cells and thus are assumed to have double the cost $2\beta$ for relying on the microenvironment for survival and growth and thus have a fitness of $1–2\beta$. The S cells can interact with tumour cells. In a normal situation, this population are relatively growth quiescent with low rates of proliferation and death. For this reason the fitness benefit derived by stromal cells from the interactions with tumour cells is assumed to be zero. However, they are able to undergo rapid proliferation and produce growth factors if they are stimulated by factors produced by I cells, giving rise to a mutualistic relationship. This relationship is represented by the parameter $\alpha$. A low $\alpha$ represents tumours in which the stroma cannot be co-opted and vice versa. From these variables, the fitness of each cell population $(W\left(S\right), W\left(I\right), W\left(D\right))$ is as follows:

\begin{equation} W\left(S\right) = 1 + f_3\alpha (\#eq:fitnessS) \end{equation}

\begin{equation} W\left(I\right) = 2 - \gamma (\#eq:fitnessI) \end{equation}

\begin{equation} W\left(D\right) = 1 + \left(1 - f_2 - f_3\right)\left(1 - \beta + \alpha\right) +\\ f_2\left(1 - \beta + \rho\right) + f_3(1 - 2\beta) + \\ 1 - \beta + \alpha + f_2\left(\rho - \alpha\right) - f_3\left(\beta + \alpha\right) (\#eq:fitnessD) \end{equation}

where $f_2$ is the frequency of I cells and $f_3$ is the frequency of D cells at a given time. All fitness equations start from balance by the sum of 1. ### Simulations First, we define the fitness of the different genotypes (see Equations \@ref(eq:fitnessS), \@ref(eq:fitnessI) and \@ref(eq:fitnessD)) through the function _fitness_rel_ that builds a data frame. It is important to note that this program models a situation where, from a WT cell population, the rest of the cell population types are formed. However, this model has also stromal cells that are not formed from a WT, since they are not tumour cells although interacting with it. Hence, for this model, we can not represent scenarios with total biological accuracy, something that we must consider when interpreting the results. ```{r example5_fitness, message=FALSE} fitness_rel <- function(a, b, r, g, gt = c("WT", "S", "I", "D")) { data.frame( Genotype = gt, Fitness = c("1", paste0("1 + ", a, " * f_D"), paste0("1 + 1 - ", g), paste0("1 + (1 - f_I - f_D) * (1 - ", b, " + ", a, ") + f_I * (1 - ", b, " + ", r, ") + f_D * (1 - 2 * ", b, ") + 1 - ", b, " + ", a, " + f_I * (", r, " - ", a, ") - f_D * (", b, " + ", a, ")")) ) } ``` Then, we are going to model different scenarios that represent different biological situations. In this case, we are going to explain four possible situations. **Note:** for these simulations the values of paratemers are normalised in the range (0 : 1) so 1 represents the maximum value for any parameter being positive of negative to fitness depending on the parameter. #### First scenario In this simulation, we are modelling a situation where the environment is relatively resource-poor. In addition, we set a intermediate cooperation between D-D and D-I and a very low benefit from coexistence of D with I. * $\alpha$ (a) = 0.5: intermediate cooperation between D and D cells. * $\beta$ (b) = 0.7: relatively resource-poor microenvironment. * $\rho$ (p) = 0.1: low benefit of D cells. * $\gamma$ (g) = 0.8: high cost of independence of I cells. We can observe that high values of $\alpha$ and low values of $\rho$ are translated in a larger profit of D cells from his interaction with S cells than from his interaction with I cells. Also, because of the high cost of independence of I cells ($\gamma$), it is not surprise that this population ends up becoming extinct. Finally, the tumour is composed by two cellular types: D and S cells. ```{r example5scen1, message=FALSE} scen1 <- allFitnessEffects(genotFitness = fitness_rel(a = 0.5, b = 0.7, r = 0.1, g = 0.8), frequencyDependentFitness = TRUE, frequencyType = "rel") set.seed(1) simulScen1 <- oncoSimulIndiv(scen1, model = "McFL", onlyCancer = FALSE, finalTime = 70, mu = 1e-4, initSize = 5000, keepPhylog = TRUE, seed = NULL, errorHitMaxTries = FALSE, errorHitWallTime = FALSE) op <- par(mfrow = c(1, 2)) plot(simulScen1, show = "genotypes", type = "line", main = "First scenario", cex.main = 1.4, cex.lab = 1.1, las = 1) plot(simulScen1, show = "genotypes", main = "First scenario", cex.main = 1.4, cex.lab = 1.1, las = 1) par(op) ``` To understand the stability of the results, we should run multiple simulations. We will not pursue that here. Note that the results can be sensitive to the initial population size and the mutation rate. #### Second scenario In this case, we set $\alpha$ lower than in the first scenario and we enable the indepenence of I cells through a lower $\gamma$. * $\alpha$ (a) = 0.3: low cooperation between D cells. * $\beta$ (b) = 0.7: relatively resource-poor microenvironment. * $\rho$ (r) = 0.1: low benefit of D cells from coexisting with I cells. * $\gamma$ (g) = 0.7: lower cost of independence of I cells than in the first scenario. Because of we are easing the possibility of independence of I cells, instead of extinguishing as in the first scenario, they compose the bulk of the tumour along with D cells in spite of the low benefit of cooperation between them (low $\rho$). Besides, we can observe that the population of I cells is bigger than the population of D cells, being at the end of the simulation in balance. On the other hand, stromal cells drop at the beginning of the simulation. ```{r example5scen2, message=FALSE} scen2 <- allFitnessEffects(genotFitness = fitness_rel(a = 0.3, b = 0.7, r = 0.1, g = 0.7), frequencyDependentFitness = TRUE, frequencyType = "rel") set.seed(1) simulScen2 <- oncoSimulIndiv(scen2, model = "McFL", onlyCancer = FALSE, finalTime = 70, mu = 1e-4, initSize = 4000, keepPhylog = TRUE, seed = NULL, errorHitMaxTries = FALSE, errorHitWallTime = FALSE) op <- par(mfrow = c(1, 2)) plot(simulScen2, show = "genotypes", type = "line", main = "Second scenario", cex.main = 1.4, cex.lab = 1.1, las = 1) plot(simulScen2, show = "genotypes", main = "Second scenario", cex.main = 1.4, cex.lab = 1.1, las = 1) par(op) ``` #### Third scenario In this case, we have a extreme situation where the microenvironment is rich (high $\beta$) and the independence costs are very low ($\gamma$) in relation with the previous scenarios. * $\alpha$ (a) = 0.2: low cooperation between D cells. * $\beta$ (b) = 0.3: rich microenvironment, being beneficial for D cells. * $\rho$ (r) = 0.1: low benefit of D cells from coexisting with I cells. * $\gamma$ (g) = 0.3: independence costs very low, being beneficial for I cells. Although $\gamma$ and $\rho$ are very low, which could make us think that I cells will control the tumour, we can observe that the fact that the microenvironment is very rich (with a low value of $\beta$) allows to D cells lead the progression of the tumour over the rest of cell populations, including I cells. ```{r example5scen3, message=FALSE} scen3 <- allFitnessEffects(genotFitness = fitness_rel(a = 0.2, b = 0.3, r = 0.1, g = 0.3), frequencyDependentFitness = TRUE, frequencyType = "rel") set.seed(1) simulScen3 <- oncoSimulIndiv(scen3, model = "McFL", onlyCancer = FALSE, finalTime = 50, mu = 1e-4, initSize = 4000, keepPhylog = TRUE, seed = NULL, errorHitMaxTries = FALSE, errorHitWallTime = FALSE) op <- par(mfrow = c(1, 2)) plot(simulScen3, show = "genotypes", type = "line", main = "Third scenario", cex.main = 1.4, cex.lab = 1.1, las = 1) plot(simulScen3, show = "genotypes", main = "Third scenario", cex.main = 1.4, cex.lab = 1.1, las = 1) par(op) ``` #### Fourth scenario This is a variation of the third scenario to illustrate that, if we set a microenvironment more rich than in the previous scenario, we get a cooperation between D and I cells, although we can still observe the superiority of D cells over I cells. * $\alpha$ (a) = 0.2: low cooperation between D cells. * $\beta$ (b) = 0.4: lower richness than in the previous scenario. * $\rho$ (r) = 0.1: low benefit of D cells from coexisting with I cells. * $\gamma$ (g) = 0.3: independence costs very low, being beneficial for I cells. ```{r example5scen4, message=FALSE} scen4 <- allFitnessEffects(genotFitness = fitness_rel(a = 0.2, b = 0.4, r = 0.1, g = 0.3), frequencyDependentFitness = TRUE, frequencyType = "rel") ## Set a different seed to show the results better since ## with set.seed(1) the progression of I cells was not shown set.seed(2) simulScen4 <- oncoSimulIndiv(scen4, model = "McFL", onlyCancer = FALSE, finalTime = 40, mu = 1e-4, initSize = 4000, keepPhylog = TRUE, seed = NULL, errorHitMaxTries = FALSE, errorHitWallTime = FALSE) op <- par(mfrow = c(1, 2)) plot(simulScen4, show = "genotypes", type = "line", main = "Fourth scenario", cex.main = 1.4, cex.lab = 1.1, las = 1) plot(simulScen4, show = "genotypes", main = "Fourth scenario", cex.main = 1.4, cex.lab = 1.1, las = 1) par(op) ``` In this case, we can see that there is more variation in the size of population of I cells. There are cases where the I cell population cooperates with D cells and, in others, there is not cooperation. You can examine this running multiple simulations (or manually rerun the example above changing the seed). ## Evolutionary Dynamics of Tumor-Stroma Interactions in Multiple Myeloma {#edmyel} This example is based on @sartakhti2016. The authors provide a frequency-dependent model to study the growth of malignant plasma cells in multiple myeloma. Assuming that cancer cells and stromal cells cooperate by exchanging diffusible factors, the study is carried out in the framework of evolutionary game theory. We first need to define a payoff strategy for this kind of scenario. The following definitions are needed: * There are $n$ phenotypes in a population denoted by $\{ P_1, \ldots, P_n \}$. * Each phenotype can produce one diffusible factor $\{ G_1, \ldots, G_n \}$. * Each diffusible factor $j$ has a different effect $r_{i,j}$ on the other phenotypes $i$. * The cost for $P_i$ for growth factor $G_i$ is denoted as $c_i$. * $M$ is the number of cells within the diffusion range. * There are $M_j$ individuals of type $P_j$ among the other group members. Then, the payoff for strategy $P_j$ is: \begin{align} \pi_{P_j}(M_1,\ldots,M_n)=\frac{(M_j+1)\times c_j}{M}r_{j,j} + \sum_{i=1, i \neq j}^n \frac{M_i \times c_i}{M}r_{j,i} - c_j \;. \end{align} In multiple Myeloma we have three different types of cells that have autocrine and paracrine effects on the cells within their diffusion range: Malignant plasma cells (MM), Osteoblasts (OB) and Osteoclasts (OC). The specification of fitness is the following (see [[10]](#Sartakhti)): \begin{align} W_{OC} &= \frac{(f_{1}(M-1)+1)c_{1}}{M}r11+ \frac{((1-f_{3})(M-1)-f_{1}(M-1)-1)c_{2}}{M}r12 \\ &+\frac{(M-(1-f_3)(M-1))c_{3}}{M}r13-c_{1} \\ W_{OB} &= \frac{(f_{2}(M-1)+1)c_{2}}{M}r22+ \frac{((1-f_{1})(M-1)-f_{2}(M-1)-1)c_{3}}{M}r23 \\ &+\frac{(M-(1-f_{1})(M-1))c_{1}}{M}r21-c_{2} \\ W_{MM} &= \frac{(f_{3}(M-1)+1)c_{3}}{M}r33+ \frac{((1-f_{2})(M-1)-f_{3}(M-1)-1)c_{1}}{M}r31 \\ &+\frac{(M-(1-f_{2})(M-1))c_{2}}{M}r32-c_{3} \;, \end{align} where $f_1$, $f_2$ and $f_3$ denote the frequency of the phenotype OC, OB and MM in the population. The multiplication factors for diffusible factors produced by the cells are shown in the following table (taken from @sartakhti2016): ```{r smyeloma, eval=TRUE,echo=FALSE, fig.cap="Multiplication factor in myeloma interaction. Table 1 in Sartakhti et al., 2016, 'Evolutionary Dynamics of Tumor-Stroma Interactions in Multiple Myeloma' https://doi.org/10.1371/journal.pone.0168856 ."} knitr::include_graphics("Myeloma_interaction.png") ``` Several scenarios varying the values of the parameters are shown in @sartakhti2016. Here we reproduce some of them. ### Simulations First, we define the fitness of the different genotypes through the function _fitness_rel_. ```{r smyelo56j} f_cells <- function(c1, c2, c3, r11, r12, r13, r21, r22, r23, r31, r32, r33, M, awt = 1e-4, gt = c("WT", "OC", "OB", "MM")) { data.frame(Genotype = gt, Fitness = c( paste0("max(0.1, 1 - ", awt, " * (f_OC + f_OB+f_MM)*N)"), paste0("1", "+(((f_OC * (", M, "-1)+1)*", c1, ")/", M, ")*",r11, "+((((1-f_MM) * (", M, "-1)-f_OC*(", M, "-1)-1)*", c2, ")/", M, ")*", r12, "+(((", M, "-(1-f_MM)*(", M, "-1))*", c3, ")/", M, ")*", r13, "-", c1 ), paste0("1", "+(((f_OB*(", M, "-1)+1)*", c2, ")/", M, ")*", r22, "+((((1-f_OC)*(", M, "-1)-f_OB*(", M, "-1)-1)*", c3, ")/", M, ")*", r23, "+(((", M, "-(1-f_OC)*(", M, "-1))*", c1, ")/", M, ")*", r21, "-", c2 ), paste0("1", "+(((f_MM*(", M, "-1)+1)*", c3, ")/", M, ")*", r33, "+((((1-f_OB)*(", M, "-1)-f_MM*(", M, "-1)-1)*", c1, ")/", M, ")*", r31, "+(((", M, "-(1-f_OB)*(", M, "-1))*", c2, ")/", M, ")*", r32, "-", c3 ) ) ,stringsAsFactors = FALSE ) } ``` It is important to note that, in order to exactly reproduce the experiments of the paper, we need to create an initial population with three different types of cell, but we do not need the presence of a wild type. For this reason, we will increase the probability of mutation of the wild type, which will disappear in early stages of the simulation; this is a procedure we have used before in several cases too (e.g., \@ref(hawkdove)). This is something that we must consider when interpreting the results. ### Scenario 1 Here we model a common situation in multiple myeloma in which $c_1