A Hybrid Computation Model to Describe the Progression of Multiple Myeloma and Its Intra-Clonal Heterogeneity

: Multiple myeloma (MM) is a genetically complex hematological cancer that is characterized by proliferation of malignant plasma cells in the bone marrow. MM evolves from the clonal premalignant disorder monoclonal gammopathy of unknown signiﬁcance (MGUS) by sequential genetic changes involving many different genes, resulting in dysregulated growth of multiple clones of plasma cells. The migration, survival, and proliferation of these clones require the direct and indirect interactions with the non-hematopoietic cells of the bone marrow. We develop a hybrid discrete-continuous model of MM development from the MGUS stage. The discrete aspect of the model is observed at the cellular level: cells are represented as individual objects which move, interact, divide, and die by apoptosis. Each of these actions is regulated by intracellular and extracellular processes as described by continuous models. The hybrid model consists of the following submodels that have been simpliﬁed from the much more complex state of evolving MM: cell motion due to chemotaxis, intracellular regulation of plasma cells, extracellular regulation in the bone marrow, and acquisition of mutations upon cell division. By extending a previous, simpler model in which the extracellular matrix was considered to be uniformly distributed, the new hybrid model provides a more accurate description in which cytokines are produced by the marrow microenvironment and consumed by the myeloma cells. The complex multiple genetic changes in MM cells and the numerous cell-cell and cytokine-mediated interactions between myeloma cells and their marrow microenviroment are simpliﬁed in the model such that four related but evolving MM clones can be studied as they compete for dominance in the setting of intraclonal heterogeneity


Introduction
Multiple myeloma (MM) is initiated through the acquisition of genetic changes that transform the plasma cells from normal to malignant.These changes often result in the development of selective advantage, leading to the excessive proliferation of myeloma cells.The development of MM leads to several harmful clinical conditions including anemia, renal failure, recurrent infections, hypercalcemia, and osteoporosis with bone fractures.These pathological conditions can frequently result in the death of the patient.
The evolution of myeloma from the clonal premalignant plasma cell condition termed monoclonal gammopathy of unknown significance (MGUS) through the intermediate stage of smoldering myeloma to the malignant MM stage is mediated by multiple sequential genetic changes, including chromosomal translocations, hyperdiploidy, and mutations, which permit independent growth and spread of MM in the bone marrow [1,2].A wide variety of genetic changes involving many different genes have been documented in MM cases.In addition to these sequential genetic changes, MM cells require specific interactions with the non-hematopoietic cells of the bone marrow including the stromal cells (BMSCs), osteoblasts, osteoclasts, and cells associated with vascular supply of the marrow [1,2].These marrow microenvironment interactions with the MM cells include direct binding of cell surface adhesion proteins and their binding partners on the other cell types as well as diffusible, soluble molecules that are secreted by one cell type and are bound and internalized by another cell type.These soluble molecules are often chemokines and cytokines produced by other marrow cells, and their receptors are expressed on the MM cells.However, specific molecules produced by MM cells indirectly affect growth of the malignancy by inducing localized resorption of bone, death of hematopoietic cells, and expansion of vessels supplying blood to the MM.The combined effects of sequential genetic changes within the evolving myeloma cells and interactions of the MM cells with the marrow microenvironment affect the migration, proliferation, and survival of the MM cells.
There are three types of mathematical models of tumor growth and intra-clonal heterogeneity.The first type is the continuous model in which partial and ordinary differential equations are used to describe the density of premalignant and malignant cells [3][4][5].Some continuous models study the effects of nutrients and inhibitors on tumor growth that may or may not be accompanied by necrosis [6,7], while others focus on the nonlinear three-dimensional spatial aspects of tumors [8].The question of clonal competition within tumors has been studied by other continuous models [9,10], and the role of clonal heterogeneity in drug resistance has been previously explored [11].While these continuous models sustain mathematical analysis and provide key insights in the development of tumors, they do not describe the cell-cell interactions and intracellular regulation of single cells.The second type of mathematical models describing tumor growth is the discrete model which usually represents cells as individual objects that interact with each other.These interactions determine the fate of each cell [12][13][14][15][16][17].Cellular automatons provide a simple yet powerful tool to understand the mechanisms of tumor growth [18,19].The discrete models usually describe cells in on-lattice or off-lattice grids.They focus on the interactions between cells and how they lead to the development of tumors, but discrete models do not link the different processes taking place at different scales.The third type of tumor growth model, which is known as a hybrid model, combines the features of both continuous and discrete models, thereby benefiting from the advantages of both of them [20][21][22].In hybrid models, the concentrations of intracellular and extracellular proteins are described by partial and ordinary differential equations while the cells are introduced as individual objects.The heterogeneous nature of tumors has been studied by some of these models [23][24][25][26].The models describing MM development usually study the impact of MM on different biological systems such as bone remodeling [27] and erythropoiesis [28].
In this work, we develop a hybrid model describing the development and intra-clonal heterogeneity in MM as described in Walker et al. [29].In this model, malignant cells are represented by elastic spheres that can move, grow, interact, divide and die by apoptosis.Their fate is determined by intracellular and extracellular regulation networks.The growth of the cellular radii between cell divisions and the formation of two daughter cells with each division lead to the expansion of the tumor.To investigate the role of the bone marrow (BM) microenvironment in the growth of myeloma tumors, we specifically consider another population of cells, the bone marrow stromal cells (BMSCs).These cells secrete cytokines such as stromal cell-derived factor 1 (SDF-1), insulin-like growth factor 1 (IGF-1),and interleukin 6 (IL-6) which are necessary to the survival, homing and growth of malignant plasma cells.The concentration of each of these cytokines is modelled by a reaction-diffusion equation.The intracellular regulation of malignant cells is represented by a system of ordinary differential equations which depends on the concentration of extracellular cytokines.We characterize each cell by a specific genotype, and we consider that after its division, daughter cells will either inherit the same genotype or acquire a slightly different one due to random mutations.As the result, aggressive clones emerge during tumor progression in a parallel pattern.Although experiments show the possibility of clones emerging in both linear and branching evolutionary patterns, we restrict this study to the parallel evolutionary case.The MM cells consume cytokines, leading to competition between clones for the available resources.We study the dynamics of clonal competition and its role in the progression of MM.

Hybrid Model of MM Development
We developed a hybrid discrete-continuous model of MGUS progression to MM.The discrete aspect of the model is observed at the cellular level: cells are represented as individual objects which move, interact, divide, and die by apoptosis.Each of these actions is regulated by intracellular and extracellular processes described by continuous models.The hybrid model consists of the following submodels: cell motion due to chemotaxis, intracellular regulation of plasma cells, extracellular regulation in the bone marrow, and acquisition of mutations upon cell division.It is an extension of a previously developed simpler model [30].While the previous model considers the cytokines in the extracellular matrix to be uniformly distributed, the present study provides a more accurate description by considering that these cytokines are produced by the BMSCs and consumed by the myeloma cells.Furthermore, a more detailed description of the intracellular pathways regulating the fate of plasma cells is provided in the present model.
The model is based on the direct effects of sequential genetic changes and marrow microenvironmental chemokine and cytokine activity that influence the chemotaxis, proliferation and survival of MM, but does not include the MM effects on bone resorption, hematopoietic cell loss, or development of specialized vasculature.The complex multiple genetic changes in MM cells and the numerous cell-cell and cytokine-mediated interactions between myeloma cells and their marrow microenviroment are simplified in the model so that four related but evolving clones develop in a process termed intra-clonal heterogeneity [31,32] (See Figure 1).Competition among these four MM clones is based on differences in cellular growth and survival rates and interactions with the marrow microenvironment.This competition results in predominance of the more fit clones and decline and ultimate extinction of the less fit ones.An early event in the MM model (Figure 1a) is a standard-risk genetic change, a t (11;14) translocation that involves the immunoglobulin heavy chain switch region on chromosome 14, inducing overexpression of the gene encoding cyclin D1, a regulator of cell cycle progression located on chromosome 11 [33].This t(11;14) clone has deregulated cell cycle progression due to the translocation, and it undergoes two separate secondary genetic events: mutations involving the oncogenes N-RAS and K-RAS, which are common secondary events in the development of multiple myeloma [2].The resultant clone with t(11;14)/mutant N-RAS has significantly enhanced proliferation compared to the parent t(11;14) clone, while the t(11;14)/ mutant K-RAS clone has much less of an increase in proliferation relative to the t(11;14) parental clone.Thus, the two descendant clones have differing degrees of increased RAS activity, resulting in a proliferative advantage for the t(11;14)/mutant N-RAS clone compared to both the t(11;14) parental clone and t(11;14)/ mutant K-RAS clone.A subsequent genetic event in the t(11;14)/K-RAS clone is, however, a mutation in the gene encoding interferon regulatory factor 4 (IRF4).Mutant IRF4 is a protein that can enhance survival and proliferation of MM cells.IRF4 mutation has been associated with mutations of N-RAS or K-RAS [32], and it allows the t(11;14)/K-RAS mutant/IRF4 mutant clone in the model to compete more successfully with the t(11;14)/N-RAS clone than either the parent t (11;14)   Among the many possible microenvironmental factors that may influence the growth of MM in the marrow, the simplified model includes the chemokine stromal cell-derived factor 1 (SDF-1).It is produced by multiple cell types but mainly by BMSCs.Other extracellular cytokines included in the model are interleukin 6 (IL-6) and insulin-like growth factor 1 (IGF-1), which are produced by other types of non-hematopoietic cells in the marrow as well as BMSCs.The marrow stromal cells are considered as sources of SDF-1, IL-6, and IGF-1 in the model (Figure 1a).SDF-1, through its receptor CXCR4, mediates the homing of circulating myeloma cells to the marrow and their migration within the marrow space [34,35].IL-6 and IGF-1 induce multiple effects after binding to their respective specific surface receptors on the MM cells, but in the model we restrict our study to their effects on the RAS/Extracellular signal-regulated kinases (RAS/ERK) pathway that promotes proliferation and the phosphatidylinositol-3 kinase/protein kinase B/Forkhead in rabdomyosarcoma (Akt/FKHR) pathway that regulates apoptosis [2] (Figure 1a).

Cell Motion
We represent each cell as an elastic sphere.It contains two parts: an incompressible part and a compressible one.As the cell grows, its radius increases and hence it pushes the neighboring cells.The motion of each cell is determined by Newton's second law.Let us denote the cell number by i.Then we have the following equation for the coordinate x i of the center of the i − th cell: where m is the mass of cell, µ is the friction coefficient, f ij is the interaction force between the cells i and j, f ch is the chemotactic force which depends on the concentration of SDF-1.We consider f ij in the following form: Here h ij is the distance between the center of the cells i and j, h 0 is the sum of their radii, K is a positive parameter and h 1 represents the incompressible part of each cell.The numerical implementaton of Newton's equation is presented in Appendix A.1.The chemotactic force represents the motion of the cell in response to SDF-1 stimulus.Let us denote the concentration of SDF-1 by S. Then the expression of f ch is given by: where κ is a positive constant and ∇S denotes the SDF-1 gradient.

Extracellular Regulation
We consider a square two-dimensional computational domain with the side equal to 1000 microns.The hybrid model contains two types of cells: myeloma cells and BMSCs.The former have a diameter of 10 microns while the latter are considered to have six-fold larger radii as observed in experimental data.Furthermore, they secrete cytokines that are necessary for the survival and proliferation of myeloma cells such as SDF-1, IL-6, and IGF-1.Let us denote the total concentration of the last two cytokines by I and the concentration of SDF-1 by S. Assuming that these cytokines have the same diffusion coefficient, production rate, consumption rate, and degradation rate, we describe their concentrations in the extracellular matrix by: where D is the diffusion coefficient, W the production factor, λ the consumption rate by each myeloma cell, and σ the degradation rate.We set the Dirichlet boundary condition I = 0 and S = 0 at all boundaries.We have prescribed zero Dirichlet boundary condition to represent the local dynamics of tumor growth in a bone marrow site surrounded by adipose tissue.The finite difference method was used to implement the reaction-diffusion equations.The used numerical methods and implementation algorithm are presented in Appendix A.2.

Intracellular Regulation
Although myeloma cells are genetically complex [2], we will restrict the intracellular regulation to two pathways (Figure 1a).These pathways are the RAS/ERK pathway which is responsible for cell proliferation and the Akt/FKHR pathway which regulates its survival.Other pathways such as the janus kinase/signal transducer and activator of transcription (JAK/STAT) pathway are not considered in the model.We have chosen these two pathways in order to study the role of aquired RAS mutations in the progression of MM and how it is affected by the extracellular matrix.Let us denote the concentrations of ERK, Akt, and FKHR as e, a, and f , respectively.We describe these concentrations inside each cell as follows: where α i and β i , i = 1, 2, 3 and γ 3 are positive constants.The coefficient α 1 (z) depends on the genotype z, and this coefficient varies depending on the effect of RAS mutations.The cell will die by apoptosis if f > f * during its lifetime cycle.If the myeloma cell survives, then it will self-renew, if e > e * by the end of its life cycle.Otherwise, the cell will die by apoptosis.We consider that the proliferation threshold e * depends on the IRF4 gene expression.The values of both intracellular and extracellular regulation parameters are provided in Appendix B.

Myeloma Cells Division and Mutations
In the present study, we investigate the progression of MGUS into MM through RAS mutations [29].We characterize each cell by an aggressiveness phenotype which represents the cumulative effects of acquired RAS mutations.We introduce an aggressiveness phenotype function f (z) resulting from the RAS genotype for each cell.After its division, the daughter myeloma cell will keep the genome of the mother cell with a probability of one-third.Otherwise, it will acquire a mutation and either increases or diminishes its aggressiveness by a positive value .We represent the genotype and resultant RAS activation that determine aggressiveness as due to describing the frequency and effects of RAS mutations in the resulting clones in Figure 2a.Furthermore, we assume that the ERK threshold for cell division is reduced upon the mutation of the IRF4 gene as shown in Figure 2b.We consider four clones of MM denoted by c 1 , c 2 , c 3 , and, c 4 (Figure 1b).The clone c 1 is the initial clone which does not harbour any RAS mutations.The clone c 2 is more aggressive than clone c 1 .It consists of MM cells that acquired the N-RAS gene mutation.The clones c 3 and c 4 belong to a lineage that is independent from the clone c 2 .They emerge when the cells of c 1 acquire the K-RAS gene mutation for clone c 3 and an additional IRF4 gene mutation for clone c 4 .The latter clone is therefore more aggressive than the former.We show the relationships between these different clones in Figure 1b.

Results
Initially, there are 47 myeloma cells in a bone marrow site containing a single BMSC, although the initial number of these cells is not essential for the results of the simulations.All of these malignant cells are considered to be tumor-initiating cells because they belong to the same initial clone.The motion of the malignant cells depends on their distance from the BMSC.Closer cells move faster than farther ones due to the increase in SDF-1 gradient near the BMSC.During their motion, myeloma cells consume the cytokines necessary for their survival such as IL-6 and IGF-1.These cytokines promote survival by activating the RAS/ERK pathway, and dowregulating the Akt/FKHR pathway.Depending upon its initial distance from the BMSC, the myeloma cell will either survive and get closer to the BMSC or it will die by apoptosis.The surviving myeloma cells will surround the BMSC and form the tumor niche.Initially, this niche consists of cells belonging to the same clone.After some time, other subclones will emerge in the course of tumor growth.We have shown the different steps of myeloma cells homing to a BMSC in Figure 3.The initial condition for the simulation is shown in Figure 3a.Due to the low number of myeloma cells in the beginning of the simulation, the BMSC secretes a high concentration of the chemokine SDF-1 which attracts the myeloma cells located close to it (Figure 3b).The surviving myeloma cells adhere to the BMSC and surround it while other cells divide and remain farther away (Figure 3c).During this process, cells whose aggressiveness phenotype belongs to the intermediate state between the clones c 1 − c 2 and c 1 − c 3 emerge (Figure 3d) and the global concentration of SDF-1 starts to reach stability due to its consumption by myeloma cells (Figure 4a).MGUS progresses to aggressive MM through the acquisition of random mutations by myeloma cells.In this respect, more aggressive subclones emerge in the course of MGUS progression due to the acquisition of RAS-associated mutations.The aggressive subclones need lower concentrations of IL-6 and IGF-1 to survive, resulting in the expansion and the persistence of the tumor (Figure 4b).Because of their rapid proliferation, they crowd out other clones and consume most of the available resources.Ultimately, these less aggressive subclones remain relatively stable in growth pattern.Still, a few of them manage to survive in proximity to the BMSC where there is a higher concentration of the IL-6 and IGF-1 cytokines, while the more aggressive subclones occupy the outer region of the tumor.The different steps of tumor development and intraclonal competition are shown in Figure 5.As the tumor progresses, the concentrations of IL-6 and IGF-1 become stable due to their consumption by the growing number of myeloma cells.As a result, the aggressive subclones expand to the detriment of the less aggressive ones.By the end of the simulation, we see that the tumor size is stable but the competition between clones is still in progress (Figure 5c vs. Figure 5d).The emergence of more aggressive clones not only reduces the populations of the less aggressive ones but also increases the total number of malignant cells.These clones survive in areas farther from the BMSC, making the tumor expand.We show the total number of malignant cells over time in Figure 4b.The proportions of each clone in the total population show that in the end of the simulation, the clone c 1 only represents ≈5% of the total population.The subclone c 2 is predominant, with ≈60% of cells (Figure 6).Overall, the global population of malignant cells increases with the emergence of more resistant subclones.The emergence of the subclone c 2 leads to an increase in the number of malignant cells because the cells belonging to this subclone can survive far from the BMSCs.The subclone c 3 emerges few hours after c 2 but barely manage to survive because of limited resources.It gives rise to the subclone c 4 , which is as aggressive as the subclone c 2 due to the additional IRF4 mutation.As the tumor progresses, these two subclones c 2 and c 4 become predominant because they are better adapted to survive in sites with limited resources.By the end of the simulation, the number of malignant cells oscillates around a stable value as well as the total concentration of cytokines in the domain (Figure 4).We show the percent of the malignant cell population occupied by each clone over time in Figure 6.

Discussion
Bone marrow stromal cells play an important role in the pathogenesis of MM.After extravasation into the bone marrow, MGUS and myeloma cells migrate and home to the areas surrounding BMSCs.This homing process is mediated by the chemokine SDF-1.It promotes the chemotaxis and homing of myeloma cells upon the interaction with its receptor CXCR4 [36].However, it is crucial for myeloma cells to be sufficiently close to BMSCs to be attracted by the secreted SDF-1.Otherwise, they will die if they are left without enough resources to survive and divide.The progression of MGUS to MM with t(11;14) translocation is marked by the emergence of new aggressive subclones which is crucial for the expansion of the tumor because it allows its adaptation to limited resources.As the tumor progresses, the aggressive subclones become predominant because lower concentrations of IL-6 and IGF-1 allow them to survive.
To study the development of MM tumors and its intraclonal heterogeneity, we developed a specific hybrid discrete-continuous model.The model was able to reproduce the experimental results presented in [29].We have used it to investigate the central role of the BMSCs in myeloma cell homing and the progression of MGUS to MM. BMSCs participate in the homing of myeloma cells and provide them with the necessary cytokines for their survival and proliferation.Numerical simulation results suggest that the initial distance between the BMSC and infiltrating myeloma cells is of paramount importance for the survival of the latter.After their homing, new, more aggressive myeloma subclones emerge due to the acquisition of RAS mutations.Other gene mutations such as the one acquired by the IRF4 gene further increase the aggressiveness of some of these subclones.They compete with each other for the available cytokines, resulting in the predominance of the more fit subclones.To better quantify the heterogeneity of clones during MM progression, we show the kernel density plots at different stages of MGUS to MM progression in Figure 7.The x-axis shows the aggressiveness phenotype function scaled from 0 to 1.These results show the predominating subclones at the different stages of tumor progression.Our findings suggest that the total number of malignant cells oscillates around a stable value after a few weeks of MM development in agreement with the in vivo experiments conducted in [37].However, we speculate that MM tumors can further expand due to other mechanisms such as stimulation of IL-6 and IGF-1 production by BMSC, the migration to other BMSCs [38], or the emegence of more resistant clones due to other mutations.The multi-scale model follows a systems biology approach [39] by integrating different interfering biological processes in one model.To understand complex phenomena such as MGUS progression to MM, it is important to properly use available data in describing events at the single cell level where processes such as gene expression and mutations take place as well as the larger tissue level where each cell interacts with its environment.Systems biology approaches focus on the whole system rather than the sum of its components.Similarly, this study is more focused on the behavior of the global model of MM homing and tumor growth than the various individual processes which regulate it.
Although the model presented here was able to reproduce the experiments describing MM intra-clonal heterogeneity, it has several limitations.Firstly, the intracellular regulation network was limited to two pathways while in reality MM is more biochemically complex [2].Furthermore, the number of considered mutations is also reduced compared to those found in MM cases.Interactions of the myeloma cells and BMSCs with other components of the bone marrow such as the osteoblasts and osteoclasts were not included in the present model.Other mechanisms affecting myeloma cell proliferation were introduced implicitly as a random perturbation of the cell cycle in the myeloma cells.The simplification of the myeloma cell transduction pathways and mutations included was due to limitations of computational power and the difficulty of studying complex models.We have limited the model to two pathways of intracellular regulation because we wanted to study the effect of cytokines IL-6 and IGF-1 on the fate of the myeloma cells and how that fate is affected by RAS mutations and the IRF4 mutation that is associated with mutant RAS genes.Also, we assumed that IL-6 and IGF-1 play the same role in the activation of the ERK/RAS and Akt/FKHR pathways.This hypothesis was made because it is difficult to distinguish the interfering actions of the different cytokines in the cell culture experiments.Finally, the study was restricted to a site of the bone marrow containing a single BMSC.Other configurations with more BMSCs and other marrow cell types can be considered as well.In future works, we will introduce the action of treatment and study the role of the intraclonal heterogeneity in MM drug resistance.
We split the first equation of the problem (P 1 ) into two sub-steps as follows: We solve the first equation for each fixed j to obtain S n+1/2 .Next, we solve the second to obtain S n .Let us consider the first equation: .
Therefore, we can write the first equation of the system (A2) in the form: a i,j S n+1/2 i−1,j + b i,j S n+1/2 i,j + c i,j S n+1/2 i+1,j = f i,j , For each fixed j, j = 1, 2, ..., N y − 1, we solve numerically: with the boundary conditions S n i=0 = 0, S n i=Nx = 0. We solve the Equation (A3) using Thomas algorithm.For that, we write the left boundary condition S i=0 = 0 as follows: where L 1/2 = 0 and K 1/2 = 0.Then, from the Equation (A3) for i = 1: .We continue for i = 2, 3, ..., N x − 1: where . We first obtain the coefficients L i+1/2 , K i+1/2 using the Equation (A3).Next, we find the solution S n+1/2 for the sub-step n + 1/2 by backward sweep using the Equation (A4).We apply the same procedure on the second equation of the system (A2) to compute the next step solution S n .

Appendix B. Parameter Values
Because of the simplifications made to the regulation network of myeloma cells, the parameters were fitted to reproduce the experimental results.The average cell cycle of the cancerous cells was considered to be stochastic.Both the extracellular and intracellular proteins concentrations are considered to be nondimensional.

Figure 1 .
Figure 1.(a)The intracellular regulation of myeloma cells as described in the model.Bone marrow stromal cells (BMSCs) secrete the cytokines insulin-like growth factor 1 (IGF-1), interleukin-6 (IL-6), and chemokine stromal cell-derived factor 1 (SDF-1) which are necessary to the survival, homing and proliferation of myeloma cells.Via their respective receptors, IGF-1 and IL-6 activate the RAS/ERK pathway which promotes the cell proliferation.They inhibit apoptosis through the phosphatidylinositol-3 kinase/protein kinase B/Forkhead in rabdomyosarcoma (Akt/FKHR) pathway.The cell migrates and homes to BMSCs through the SDF-1/CXCR4 axis.The IRF4 mutation, which has been associated with concomitant RAS mutations, promotes survival and proliferation.BMSCs, which are much larger cells than the myeloma cells, are shown in reduced in size in this figure as well as the receptors of both IGF-1 and IL-6; (b) The parallel evolution pattern of multiple myeloma clones resulting in intra-clonal heterogeneity.More aggressive clones result from a more aggressive N-RAS mutation in clone 2 or the acquisition of IRF4 mutation in addition to the less aggressive K-RAS in clone 4. Each clone is shown by its corresponding color in the model.

Figure 2 .
Figure 2. (a) The activation rate of RAS protein as a function of the genotype function z; (b) The ERK threshold for division as a function of the genotype function z, it decreases due to the IRF4 mutation found in the clone c 4 .

Figure 3 .
Figure 3.The sequential steps of myeloma cells homing to a BMSC.Myeloma cells are represented by the cells with the smaller radii and the BMSC is the large cell in the middle.Each clone of multiple myeloma (MM) cells is denoted by a specific color.The red color gradient represents the summed concentrations of the cytokines SDF-1, IL-6, and IGF-1 shown with white (0) to red (1).(a) After their infiltration, the myeloma cells move towards the BMSC; (b) The surviving cells surround the BMSC; (c) The myeloma cells divide and form the tumor niche; (d) The tumor expands and more aggressive subclones start emerging.

Figure 4 .
Figure 4. (a) The total concentration of the IL-6 and IGF-1 cytokines over time in log scale; (b) The total number of malignant cells over time.

Figure 5 .
Figure 5. Snapshots of a simulation showing the competition between clones.(a) The tumor consists mainly of the clone c 1 (shown in yellow) with few cells in the intermediate state (in purple) between clones and the emergence of the clone c 2 (in cyan) in the sides.The size of the tumor remains limited because clone c 1 cells need relatively high concentration of I to survive; (b) Compared to the clone c 1 , the clone c 2 cells expand and survive in areas with lower concentration of cytokines; (c) The clone c 2 cells surround the tumor and crowd out the cells of the clone c 1 leading to the reduction of their population.The clones c 3 (in blue) and c 4 (in magneta); (d) The subclone c 4 is as aggressive as the subclone c 2 due to the additional IRF4 mutation and it manages to coexist with it in the remote areas with fewer cytokines.

Figure 6 .
Figure 6.The proportion of each clone in the total population over time.

Figure 7 .
Figure 7. Gaussian kernel density plots indicating the frequency of cells acquiring each specific mutation at the different stages of the simuation.(a) The initial population composed of cells belonging to the subclone c 1 ; (b) c 2 is the most predominant subclone after 16 hours of the beginning of the simulation; (c) The K-RAS subclones c 3 and c 4 emerge and start expanding; (d) The clones c 2 and c 4 represent the majority of the population because they are equally aggressive.

Table B1 .
Values of intracellular regulation parameters for myeloma cells.δ is an arbitrary length unit.

Table B2 .
Values of extracellular regulation parameters.δ is an arbitrary length unit.