Evolutionary Dynamics of Tumor-Stroma Interactions in Multiple Myeloma

Cancer cells and stromal cells cooperate by exchanging diffusible factors that sustain tumor growth, a form of frequency-dependent selection that can be studied in the framework of evolutionary game theory. In the case of multiple myeloma, three types of cells (malignant plasma cells, osteoblasts and osteoclasts) exchange growth factors with different effects, and tumor-stroma interactions have been analysed using a model of cooperation with pairwise interactions. Here we show that a model in which growth factors have autocrine and paracrine effects on multiple cells, a more realistic assumption for tumor-stroma interactions, leads to different results, with implications for disease progression and treatment. In particular, the model reveals that reducing the number of malignant plasma cells below a critical threshold can lead to their extinction and thus to restore a healthy balance between osteoclast and osteoblast, a result in line with current therapies against multiple myeloma.


Introduction
The production of growth factors is one of the most important determinants of cancer development [1]. Clones producing different growth factors can sustain each other's growth [2,3], and non-producer cells can rely on the growth factors diffusing from neighboring cells [4]. The production of growth factors by cancer cells, therefore, is a form of cooperation that can be studied in the framework of evolutionary game theory [5,6].
Cancer cells also produce diffusible factors that induce stromal cells to release other growth factors that support tumor proliferation [7]. Consider, for example, multiple myeloma, a type of cancer of plasma cells [8][9][10] in which different cell types-malignant plasma cells (MM) themselves, as well as osteoclasts (OC) and osteoblasts (OB)-contribute to bone resorption and bone formation by exchanging diffusible factors (Fig 1). In a healthy bone, osteoclasts demolish bone tissue and osteoblasts regenerate it, two processes that balance each other maintaining bone health. MM cells produce cytokines like interleukin-1 (IL-1), IL-3, tumor necrosis factor (TNF-α), receptor activator of nuclear factor-kB ligand (RANKL) and macrophage inflammatory protein (MIP-1α), which activate OC, and consequently increase resorption. MM cells also produce factors, such as IL-1 and Dickkopf-related protein 1 (DKK-1), which have inhibitory effects on OB differentiation. Multiple myeloma, therefore, alters the normal OB-OC equilibrium and induces, among other symptoms, bone fracture. Some factors like IL-3 have stimulatory effects on OC and inhibitory effects on OB differentiation. OC also secretes growth factors like MIP-1 and IL-6 that affect MM cell proliferation [11,12]. In short, each cell type has stimulatory or inhibitory effects on the other types and on itself, leading to frequency-dependent selection, which can be analyzed using evolutionary game theory.
Multiple myeloma has been analyzed before in the framework of evolutionary game theory by Dingli et al. [13] using a model with pairwise interactions, that is, assuming that the effect of the growth factors produced by each cell is limited to one other companion cell. The growth factors produced by cancer and stromal cells, however, have autocrine and paracrine effects on multiple cells, therefore it would be more appropriate to model the effect of growth factors as a multiplayer game with collective interactions, rather than pairwise. While models with pairwise interactions are often used in game theory, it is known that their results cannot always be extended to games with collective interactions, and can actually lead to misleading conclusions [14].

Fitness functions
Assume that there are n phenotypes in a population denoted by P 1 . . .P n . Each phenotype can produce one of n diffusible factors G 1 . . .G n , respectively. Each diffusible factor j has a different effect (r i,j ) on the other phenotypes i (for example, G 1 can confer a net benefit on P 1 and P 2 , and cause a cost to P 3 ). Therefore, if there are N j individuals of type P j among the other group members (where N 1 +N 1 +. . .+N n = N, and N is the number of cells within the diffusion range of the growth factors), the payoff for strategy P j is: Where c i (for i 2{1,. . .,n}) is the contribution of (that is, the cost for) P i for growth factor G i , and r i,j is the multiplication factor that specifies the effect of G j on individuals belonging to type P i . For example, in a game with three phenotypes, the above equation would yield: Where k, j, and l are the numbers of individuals that invest in the growth factors G 1 , G 2 and G 3 respectively.

Dynamics in infinite populations
In a well-mixed, infinite population we assume, as is standard, that groups are formed at random at each generation, after which fitness is calculated. Let us call x, y and z the frequencies of P 1 , P 2 and P 3 in the population, respectively (x+y+z = 1). N is the number of cells within the diffusion range of the growth factors, S is the number of individuals belonging to either type P 1 or P 2 , and m is the number of individuals of type P 1 . Fitness is calculated by weighting the payoffs obtained in the randomly formed groups, weighted by the probability that such groups occur. The probability that a group contains S-m-1 individuals of type P 2 , N-S individuals of type P 3 and m individuals of type P 1 (in addition to the focal individual) is given by B (S-1, Therefore, the payoff for P 1 is: Evolutionary Dynamics of Multiple Myeloma and the fitness function for P 1 is: Likewise, the fitness functions for P 2 and P 3 are given by: where The fitness functions can be reduced to: Let us denote with x(t), y(t) and z(t) the frequencies of phenotypes P 1 , P 2 and P 3 , respectively. We assume that the frequencies change according to the replicator dynamics [15]: where W ¼ xW 1 þ yW 2 þ zW 3 . The equilibria of the system can be found by setting Equation set Eqs 15-17 equal to zero [15].

Dynamics of diffusible factors production in multiple myeloma
We follow Dingli et al. [13] in their assumptions on the effect of each cell on another. Essentially OC and OB cells are in equilibrium in the absence of MM cells, while MM and OC cells have a stimulatory effect on each other, and MM cells inhibit OB cells and OB cells have little or no effect on MM cells [8-12, 16, 17]. Table 1 summarizes these effects, while Table 2 describes the parameters used. We can assume that r i,i = 0 that is, each cell type has no net effect on itself. This does not mean that there are no autocrine effects; it implies, instead, that production and consumption of the growth factors are linear measures of the density of the cell types (MM cells, for instance, produce IL6 that stimulate MM cell proliferation: the more MM cells, the more IL6 is available, but also used by more MM cells). By setting r i,i = 0 we imply that the balance changes only due to paracrine effects (in this case, the number of OC cells). r MM,OB = 0 because G OB (i.e., the diffusible factors produced by OB) has no effect on MM. r OB,OC = r OC,OB because OB and OC balance each other; for the same reasons r MM,OC = r OC,MM .
Given the values in Table 2 and Eqs 12-14, the fitness of the three types of cells are: Eqs 15-17 and 18-20 can be used to investigate the evolutionary dynamics of the three cell types over time.

Results
Clearly, the three vertices of the simplex (x OC = 1, x OB = 0, x MM = 0), (x OC = 0, x OB = 1, x MM = 0) and (x OC = 0, x OB = 0, x MM = 1) are fixed points of the game. Other fixed points are the three points on the edge of the simplex: and one in the interior of the simplex: Since the fixed points and Eqs 18-20 are functions of several parameters, which makes the analysis exceedingly complex, we investigate the dynamics and the stability in three more specific scenarios. In each scenario we investigate the effect of the multiplication factors and group size on the dynamics. In what follows we show examples of the results for specific values of the parameters. Increasing the value of a, b and d, or reducing the value of c 1 , c 2 , c 3 changes the trajectories and velocity of the dynamics but the direction and the position of the equilibria remains unaltered.

Scenario 1
In multiple myeloma the contributions of the three types of cells (OC, OB, and MM) are clearly different. As mentioned above, MM cells produce growth factors and cytokines at higher levels than OB and OC cells. In the presence of MM cells, OC cells are also stimulated to produce more growth factors. Hence, it seems reasonable to study a scenario in which c 2 <c 1 <c 3 . Whenever the net benefit of diffusible factors that are secreted by MM cells is greater than the benefit that OC cells can obtain through the diffusible factors produced by OB, there exists a polymorphic stable point between MM and OC. Fig 2 shows an example of this scenario. While OB and OC cells are in equilibrium in the absence of MM cells, the presence of MM cells destabilizes that equilibrium and makes the population evolve to a stable mixed equilibrium of MM and OC cells, regardless of the initial frequencies (Fig 3). As the population approaches the stable point, OB cells disappear. In this process, the risk of bone fracture dramatically increases, a typical occurrence in multiple myeloma, which the model therefore explains as a perturbation of the OB-OC equilibrium by MM cells. Even a tiny fraction of MM cells in the population is able to change the dynamics and lead to an increased risk of bone fracture. Increasing the contribution (that is, the amount of diffusible factors secreted) by OC and MM cells increases their own fitness, while for OB increasing contributions decreases fitness (Fig 4).
The size (N) of randomly composed group in our model is, essentially, the number of cells in the diffusion range of the growth factors of the focal individual. Group size has an important role in the evolutionary dynamic of the game. For example, if N = 2 (i.e. pairwise interaction) the dynamics described above changes and the game has two stable points on the OB-OC and OC-MM edges, and an interior saddle point (Fig 5). The presence of this interior saddle point reveals already a fundamental difference between the pairwise model [13] and our model with collective interactions (Fig 2): in our model a single MM cell is enough to lead the system from the healthy OB-OC equilibrium to the OC-MM equilibrium, whereas in the pairwise model a large number of MM cells are necessary. The effect of N on the position of the two fixed points of our system is described in the      If the costs for OC, OB, and MM cells are equal (i.e., c 1 = c 2 = c 3 ), then the following are additional fixed points: In this scenario we have three type of stable points: (i) a polymorphic stable point between OC and OB, (ii) a polymorphic stable point between OC and MM and (iii) two stable points on OC-OB and OC-MM edges of the simplex.    If a = 1 and N = 2 the dynamics is equivalent to the one described by Dingli et al. [13]. N however is important: while N has no effects on the position of the fixed points, it does affect their stability (a result that is not captured by models with pairwise interactions, in which N = 2): if a<b and a<2N/(N-1), all eigenvalues of the Jacobean matrix are negative at (x OC = 1/2, x OB = 1/2, x MM = 0), which is therefore a stable point of the system (and a healthy OC-OB balance); if b+d>a and b<2N/ (N-1), all eigenvalues of the Jacobean matrix are negative at (x OC = 1/2, x OB = 0, x MM = 1/2), which is therefore a stable point of the system (and an unhealthy one in which osteoblasts disappear).

Exploiting evolutionary dynamics in response to therapy
In the three scenarios we have analysed there are polymorphic equilibria either on the OC-MM border or on the OC-OB border, or a monomorphic stable point on the MM vertex. If the only stable point of the dynamics includes at least a fraction of MM cells, clearly killing MM cells only slows down the dynamics, without leading to an effective, long-lasting treatment. However, in scenarios that include the OC-OB equilibrium, changing the parameters of the system or the fraction of cell types may affect the dynamics and change the equilibrium of the system.
For example, consider the case in which there are two stable points, on the OC-OB and OC-MM borders (Fig 8b). It may be possible, by removing some of the MM cells, or by adding OB cells, to make the population evolve to the OC-OB mixed equilibrium, rather than to the MM-OC equilibrium (Fig 11). The same result can be achieved by changing the parameters of the game. Understanding that tumor-stroma interactions are a system with multiple equilibria, therefore, could help devise strategies to change the dynamics and make the tumor evolve to the healthy OC-OB equilibrium.

Differences with the model of pairwise interactions
Given the linear nature of the payoffs, it is possible, in principle, to construct a game with pairwise interactions that is equivalent to our multi-player collective-interaction game, similar to what can be done for linear public goods games [18]. The payoff of the transformed pairwise game equivalent to ours is given by Note that this is different from the pairwise game analyzed by Dingli et al. [13]; the core difference between their model and ours is that in our model interactions are collective (with N>2) rather than pairwise (N = 2). Dingli et al. [13], however, use the same assumptions we use about the effect of the growth factors produced by each type on each other cell type. That is, their model is equivalent to ours if we set N = 2, a = 1 and c 1 = c 2 = c 3 . These assumptions, however, seem unrealistic. In reality, not only are the costs different, but most importantly, clearly N>2 (that is, growth factors have a collective effect on a large number of cells, not just one cell engaged in a pairwise interaction). As we have seen, N has a strong effect on the stability of the fixed points. It is not surprising then, that our results (the number and nature of the equilibria) are not equivalent to the model with pairwise interactions [13] if N>2. Fig 12 shows some of these differences, and we discuss different cases here.
b>1: In this case the net benefit that OC cells obtain from MM cells is greater than what they get from OB cells. In the pairwise model of Dingli et al. [13] the system has only one stable point on the OC-MM edge and in this situation multiple myeloma eventually leads to bone fracture, due to the lack of OB cells (left panel of Fig 12a). According to our model, instead, the system can have different dynamics and stable points depending on the costs and the size of the group (that is, the diffusion range of the growth factors) (middle and right panels of  Fig 12a). The differences have implications for therapy: reducing the number of MM cells in our model may help to restore a healthy OB-OC balance (see middle panel of Fig 12a), an effect that is not captured by the pairwise model of Dingli et al. [13]). In our model reducing the number of MM cells results in a recovery of the normal balance OC-OB if for any p 1 and p 2 (where p 1 +p 2 = 1, p 1 ,p 2 2[0,1]): b<1, b+d<1: Dingli et al. [13] showed that in this case, in a pairwise model, the normal balance OC-OB is the only stable equilibrium point of the system (left panel of Fig 12b). This is the case in our model as well (middle panel of Fig 12b) if c 2 <c 1 <c 3 , but not if we change the values of the costs and N (right panel of Fig 12b). In our model, instead, an additional MM-OC stable equilibrium is possible. b<1, b+d>1: In this case the pairwise model [13] has two stable points on the OC-OB and OC-MM edges (left panel of Fig 12c), while our model has different dynamics and equilibria (middle and right panels of Fig 12c). In particular, in our model the OC-MM stable equilibrium can be absent, which suggests that reducing b, for instance by inhibiting MIP-1α or IL-1β [19,20] may be beneficial. This was already observed by Dingli et al. [13], and is in agreement with clinical observations [21,22]. In our model this result is more realistic under certain parameters (middle panel of Fig 12c), as the coexistence of OB and OC can be the only stable equilibrium (whereas in the pairwise model of Dingli et al. [13], there is an additional OC-MM equilibrium, and hence the population must be pushed outside its basin of attraction); on the other hand, for different parameters (right panel of Fig 12c), a pure MM population is the only stable outcome of the dynamics.

Discussion
As the interactions between cancer cells and stromal cells depend on the effect of diffusible factors with autocrine and paracrine effects, tumor-stroma interactions are frequency-dependent processes that can be analyzed in the framework of evolutionary game theory. Our analysis of tumor-stroma interactions in multiple myeloma departs from previous work [13] because of our assumption that cells are engaged in collective interactions rather than pairwise interactions [13]. Assuming collective interactions is more realistic, as growth factors produced by both cancer and stromal cells have autocrine and paracrine effects that are clearly not limited to just one other interacting cells; growth factors have collective effects on multiple cells, and each cell's fitness is influenced by the production of growth factors by all the cells within the diffusion range of the factors.
As we have seen, the dynamics of tumor-stroma cooperation revealed by our model with collective interactions is more complex that observed in a model with pairwise interactions [13], and it leads to fundamental differences in the results. Like in the pairwise model analyzed by Dingli et al. [13] we observe the co-existence of two strategies when one of the three types is introduced, even at low frequency, in the population. In the presence of MM cells, the nature of the OB-OC equilibrium can change from a stable coexistence (a healthy balance of OB and OC cells) to a saddle point leading to an OC-MM equilibrium, which explains in part the pathological condition of bones observed with multiple myeloma, as the lack of OB cells, increases the risk of bone fracture [23,24] a result also observed in the pairwise model [13].
Our results, however, show a fundamental difference. Dingli et al. [13] observe that when the net benefit that OC cells obtain from MM cells is greater than what they get from OB cells (b>1), which, as they point out, is generally the case [16,25], the only stable equilibrium is the co-existence of MM and OC cells. Only if b<1 OB and OC cells can re-establish a stable healthy equilibrium, free from MM cells. In the pairwise model [13], in short, under realistic parameters (b>1) a pathological condition is inevitable when even a single MM cell is introduced in the population. In our model, instead, even with b>1 a reduction in the amount of MM cells can induce a change leading to a balance of OC and OB and the extinction of the MM cells (Fig 11). Additionally, the same result can be achieved, in principle, by increasing the fraction of osteoblasts.
Our model suggests, therefore, that therapies that aim at reducing the number of MM cells (for instance, protease inhibitors like Bortezomib [26][27][28], which is currently used in the U.S. for treating multiple myeloma) could be effective against multiple myeloma, a result that is not captured by a pairwise model [13]. In addition, it shows that extinction of malignant plasma cells could be achieved by increasing the number of osteoblasts, a cell therapy approach that might be feasible [29]. Finally, as we have seen, our model reveals that in other scenarios the dynamics of multiple myeloma can be more complex, with additional equilibria in which MM cells persist, or the collapse of MM cells and the return to a healthy balance between OB and OC cells.
Tumor-stroma interactions are, like other frequency-dependent selection processes, not intuitive, and game theory models can help understand the logic of their dynamics and, as a consequence, the feasibility of possible therapies. Our result that the malignant plasma cells may go extinct if they fall below a critical threshold required to model multiple myeloma using a model of collective interactions, which is more realistic than the assumption of pairwise interactions used previously.
Clearly one could add further realism to the model. In particular, we have assumed that the effect of the growth factors secreted by the cells is a linear function of their concentration, while the effect of growth factor is usually a nonlinear function [4,29,30], and it is known that non-linear benefits can lead to different dynamics. We have also assumed that cells interact at random in a well-mixed population, as if multiple myeloma was a purely liquid tumor; in reality, however, multiple myeloma could be considered in part a spatially structured population, and it is known that spatial structure can also affect the dynamics of public goods, especially with linear benefits [31]. It would be interesting to understand how nonlinearities and spatial structure in a collective action model affect the dynamics of tumor-stroma interactions.
The analysis of multiple myeloma made by Dingli et al. [13] has been pivotal in introducing game theory in the study of cooperation between tumor and the microenvironment. By extending their model to a model with collective interactions, we showed some fundamental results that a simple model with pairwise interactions could not reveal. We hope this will stimulate further analysis and further use of evolutionary game theory in the study of cancer dynamics.