Multi-Scale Agent-Based Multiple Myeloma Cancer Modeling and the Related Study of the Balance between Osteoclasts and Osteoblasts

Research Background Currently, multiple myeloma is the second most common hematological malignancy in the U.S., constituting 1% of all cancers. With conventional treatment, the median survival time is typically 3–4 years, although it can be extended to 5–7 years or longer with advanced treatments. Recent research indicated that an increase in osteoclast (OC) activity is often associated withmultiple myeloma (MM) and that a decrease inosteoblast (OB) activity contributesto the osteolytic lesions in MM. Normally, the populations of OCs and OBs are inequilibrium, and an imbalance in this statecontributes to the development of lesions. Research procedures A multi-scale agent-based multiple myeloma model was developed to simulate the proliferation, migration and death of OBs and OCs. Subsequently, this model was employed to investigate the efficacy of thethree most commonly used drugs for MM treatment under the following two premises: the reduction in the progression of MM and the re-establishment of the equilibrium between OCs and OBs. Research purposes The simulated results not only demonstrated the capacity of the model to choose optimal combinations of the drugs but also showed that the optimal use of the three drugs can restore the balance between OCs and OBs as well as kill MMs. Furthermore, the drug synergism analysis function of the model revealed that restoring the balance between OBs and OCs can significantly increase the efficacy of drugs against tumor cells.


Introduction
Previous studies [1] stated that multiple myeloma is the second most common hematological malignancy in the U.S. (after non-Hodgkin lymphoma), constituting 1% of all cancers. Multiple myeloma treatmentcan be classified intothe following three methods. Thefirst is high-dose chemotherapy with autologous hematopoietic stemcell transplantation, which can prolong overall survival and evoke complete remission, but it is not curative. The second is allogeneic stem cell transplantation, which can cure MM in a small percentage of patients with significant side effects [1]. The third is chemotherapy withthe following drug combinations: 1,bortezomib, melphalan, andprednisone, with an estimated overall survival of 83% at 30 months [2]; 2,lenalidomide plus low-dose dexamethasone,with 82% survival at two years [3];and 3,melphalan, prednisone and lenalidomide, with90% survival at 2 years [4]. Patients over 65 years old and those with significant concurrent illness can only receive the thirdtreatment, but these drugs have significant side effects, and the treatment effect is not obvious. To identify novel therapeutic options for the treatment of multiple myelomascientists are investigating the multi-scale pathogenesis of multiple myelomaat the intracellular, intercellular and tissue scalesand employing molecular drugs to treat MMs.
Overall, 80-90% of myeloma patients develop bone lesions during their disease course [1]. Multiple myeloma bone disease is characterized by dysfunction of both OB-mediated bone formation and OC-mediated bone resorption [5]. Bone homeostasis is maintained by the balance between the synthesis of new bone by OBs and the removal of old bone by OCs. In MM, there is an imbalance in the proportion of OCs and OBs. OB activity is markedly decreased or absent, and OC bone resorption is activated [5,6]. In this study, the balance is defined by two standards: one, the ratio of OCs to OBs; and two, the absolute difference in the number of OCs and OBs within a reasonable interval.
Multiple interactions in the myeloma bone marrow microenvironment are responsible for myeloma bone disease. A recent study [7] demonstrated that the DKK1-Wnt-OPG/RANKL intracellular signaling pathway can mediate the balance between OBs and OCs, which has becomeone of the most important factors in the pathogenesis of multiple myeloma. There are four major scenarios for the multi-scale pathogenesis of multiple myeloma (Fig 1). I: The Wnt signaling pathway stimulates the growth, differentiation and activity of osteoblasts [8]. II: Dickkopf (DKK1) is secreted by MMs. Because DKK1 is a Wnt inhibitor, it inhibits the phosphorylation of beta-catenin to prevent its degradation [1]. Higher DKK1 expression has been found in myeloma patients and has shown a positive correlation with the advanced stages of myeloma [9]. III: DKK1 directly increasesRANKL and decreasesosteoprotegerin(OPG) expression in OBs [10]. The ratio of OPG/RANKL is negatively related to the number of OCs. IV: OCsproduce TNFα, which directly stimulates the formation of MMs and induces stromal cells to secrete factors, such as RANKL, that drive OC formation. TNFα is a potent inducer of OCs that blocks OB differentiation and promotes MM growth. MMs inhibit the growth of OBs and stimulate OCs to evoke a vicious cycle that promotes the imbalance between these two cell types. If OB formation is simultaneously inhibited by scenarios I and II and the growth of MMs is stimulated by scenario IV, the ratio of OPG/RANKL will decrease markedly, thereby escalating the generation of OCs [11].
Based on the pathogenesis of MM, several drugs have been developed to treat this disease [12]. Glucocorticoidshave been used for hematological cancer therapy, but they are associated with multiple adverse outcomes, such as the suppression of OCs and OBs [13]. BHQ880inhibits DKK1 and thus promotes bone formation, which in turn inhibits myeloma-induced osteolytic disease and MMgrowth [14]. Lidamycin accelerates the apoptosis of MMs.
However, most of these drugs do not workas well in vivo as they do in vitro due to absorption, distribution, metabolism, and toxicity (ADME) problems. For this reason, we employed multiple drugs to treat MM with an optimal drug combination plan to solve ADME problems and increase in vivo drug efficacy. However, it is impractical to quantitatively evaluate the optimaldrug combination in vivo. Hence, we employed a 3D multi-scale agent-based model that encompasses the intracellular, intercellular and tissue scales to address the following three specific aims: to develop a platform to describe the relationships among MMs, OBs and OCs; to investigate the pathogenesis of MM using this platform; and to employ this platform to identify the optimal drug combination for MM treatment.
The simulated results demonstrated that our model can be used not only to simulate the proliferation, migration and death of OBs, OCs and MMs but also to investigate the optimal use of these three drugs to inhibit MM growth and restore the balance between OCs and OBs.

Materials and Methods
To describe tumor growth with the imbalance between osteoblasts and osteoclasts and to study the response of multiple myeloma to particular drug combinations, three types of agents are denoted in the model: MMs, OCs and OBs.Our multi-scale model consists of three biological scales: intracellular, intercellular and tissue. The intracellular scale describes the fundamental mechanisms for cellular phenotypic switches, and the intercellular scale bridges the tissue and intracellular scales as follows: (a) cytokines and drugs are delivered to the tumor microenvironment at the tissue scale; (b) MMs, OBs and OCs undergo phenotypic switching in response to stimulation by specific cytokines at the intercellular scale; and (c) DKK1 secreted by MMs affects the balance between OBs and OCs as well as their migration in response to secreted cytokines and drugs at the tissue scale. A 100 Ã 100 Ã 100 three-dimensional cube with four sub-compartments was employed to represent a slice of the virtual tumor extracellular matrix (ECM). The lattice size is 5 μm, which is approximately the same as the radius of a MM. Approximately one hundred MMs, OBs and OCs were mixed and initialized in the center of the lattice, forming a sphere.The age of the MMs, OBs and OCs was randomly initialized between 0 and 24 hours.

Intracellular Scale: Phenotype switching of tumor cells as "agents"
At every simulation step (Δt = 2 hours), each MM, osteoclast orosteoblast determines its phenotype according to the following rules (Fig 2).
Apoptosis. At each time step, if the apoptosis probability of the cell (OBs, OCs or MMs) is less than the threshold, Apop rate the cell will initiate apoptosis. Any given cell takes 10 time steps to complete apoptosis and is then absorbed. As reported in previous studies [15], the cell's apoptosis rate is defined according to an exponential distribution (Eq 1): where 4t denotes the time step(which is two hours), andλ is a positive number representing the average apoptosis frequency of each cell type. Eq 2 describes the impact of Lidamycin on the apoptosis rate of MMs based on the Hill function [16]: where L ijk denotes the relative Lidamycin concentration at the grid point (i, j, k), β mm L is the maximum increase in the apoptosis rate, λ 0 mm is the basal apoptosis rate of MMs, and (K L mm ) is the threshold of Lidamycin that promotes MM apoptosis.
To restore the balance of OCs and OBs, the apoptosis rates of OCs and OBs should be regulated by drugs. Eq 3 describes how glucocorticoids (GCs) promote the apoptosis of OBs and OCs based on the Hill function [17]: where n = 1 for OBs, n = 2 for OCs, G ijk denotes the relative GC concentration at the grid point (i, j, k), and K GC n is the threshold at which GCs promote the apoptosis of OCs and OBs.
Proliferation. Eq 4 describes the proliferation probability of multiple myeloma cells based on a Hill function [18]: where T ijk denotes the relative TNFα concentration at the grid point (i, j, k), K T MM is the threshold at which TNFα promotes the proliferation of MMs, P prol MM represents the impact of TNFαon the proliferation rate of MMs, P 0 MM denotes the basal proliferation rate, and β MM T is eamum level of proliferation of MMs. Migration: A proliferating cell will migrate in the first three phases of cell cycle (G0/G1, S and G2), whereas it will look for an empty location to divide after entering the mitotic M phase.Quiescence: There are two possibilities for the cell to be in a quiescent state: the cell cannot go through the cell cycle, or the cell cannot find an appropriate free location for division. doi:10.1371/journal.pone.0143206.g002 Eq 5 describes the proliferation probability of OBs based on the Hill function [19]: where P prol ob represents the proliferation rate of OBs after responding to the positive effects of Wnt signaling, P 0 ob denotes the basal proliferation rate, and β ob is the maximum proliferation level. Eq 6 describes the proliferation probability of osteoclasts [20]: where P prol oc is the proliferation probability of osteoclasts, P 0 oc denotes the basal proliferation probability, p pathway oc is the maximum inhibition in relation to O_RL ijk , K O_RL is the threshold of O_RL ijk that inhibits the proliferation of osteoclasts, and O_RL ijk denotes the concentration ratio of OPG to RANKL at the grid point (i, j, k). Eq 7 utilizes a Bernoulli function to determine whether the cell enters the cell cycle: If C rand falls in the interval [0,p prol ), the cell enters the cell cycle and starts to proliferate; otherwise, the cell stays quiescent and awaits the next round. In the M phase of the cell cycle, a cell will divide if it finds at least one free location within the search distance.
Migration. Non-proliferating cells will migrate at each step. Proliferating cells will migratein the first three phases of the cell cycle (G0/G1, S and G2), whereas they will look for an empty location to divide after entering the mitotic M phase [21]. This will be discussed in detail in the next section.
Quiescence. There are two possibilities for the cell to be in a quiescent state: one, the cell cannot go through the cell cycle as a result of Eq 5; and two, the cell cannot find an appropriate free location for cell division.
Intercellular Scale: The dynamics of molecules in signaling pathways for each cell and the rules for choosing the "most attractive" location EWnt ijk is the impact of DKK1 on Wnt as described by Eq 8: where D ijk denotes the relative DKK1 concentration at the grid point (ijk), K D is the threshold of DKK1 that inhibitsWnt signaling, EWnt ijk represents the impact of DKK1 on Wnt signaling, and E 0 and β w are the basal and maximum inhibition levels in relation to Wnt signaling. As previously reported [22], an imbalance between OBs and OCs results when OCs outnumber OBs in the context of MM. BHQ880 inhibits DKK1 to reduce Wnt signaling, and Eq 8 is replaced by Eq 9: where D ijk e ¼ D ijk 2 D ijk þB , represents the current effective DKK1 concentration at the grid point (i, j, k), and B denotes the dose of BHQ880.
Each living MM chooses the "most attractive" location to proliferate or migrate based on the following rules: 1) A non-M phase cell at position p 0 will always search for a position with greater nutrition to migrate or divide. The six nearest (Moore) neighbors of p 0 are the candidate locations. Each one is ranked by Eq 10: where R l is the ranking score of each candidate position p ijk l , r l is the distance from the candidate location p ijk to position p 0 , and p(r l ) is the probability that the cell moves to the candidate location p ijk .
where D is the distance squared over time.
MMs try to avoid isolated or over-crowded locations, but they are not as sensitive as solid tumor cells.The preference for neighborhoods (V l ) is denoted by Eq 12: 1=16 P ijk has no neighbor cells The ranks of the candidates are normalized in Eq 13: To facilitate die-casting, all the normalized ranks are incorporated to form a scale S in Eq 14, in which each candidate corresponds to a range S l S is an ordered set of S l . Each S l is a region in the interval [0,1] and relates to the lthcandidate site. Die rolling generates a random value d 2 [0,1). If d falls in S l , the corresponding candidate locationR l will be chosen as the next migration or proliferation site.
2) If no space is available, the cell will become reversiblyquiescent and will await the next round.
C ijk (t + 1) is the cytokine value at the location P ijk in the t + 1 time step, λ C is the diffusivity, and C l ijk ðtÞ; l ¼ 1; 2; . . . ; 6 are the cytokine values of P ijk 's six nearest neighbors in the current time step. If the cytokine is DKK1 or RANKL and there is aMM at P ijk , the time-dependent characteristic function χ mm (t, P ijk ) is set to 1; otherwise, this value is set to 0.If the cytokine is OPG and there is an OB at P ijk , the time-dependent characteristic function χ mm (t, P ijk ) is set to 1; otherwise, this value is set to 0.If the cytokine is TNF-α and there is an OC at P ijk , the time-dependent characteristic function χ mm (t, P ijk ) is set to 1; otherwise, this value is set to 0. Se G and DEG represent the cytokine secretion rate andthe degradation rate, respectively.
As discussed previously, BHQ880, GCs and Lidamycin were utilized in our model to restore the balance between OCs and OBs. Additionally, we developed this novel using in silico algorithms to identify the optimum multi-drug combination to restore this balance. BHQ880 regulates the proliferation of OBs and OCs by inhibiting DKK1 [23]. GCs promote the apoptosis of both OBs and OCs. Because MMs play an important role in generating an imbalance between OBs and OCs by secreting DKK1, Lidamycin can restore the balance by increasing the apoptosis of MMs [24]. Eq 16 describes how these multiple drugs diffuse into the surrounding tissue before being taken upby MMs, OBs and OCs.
where DR ijk (t + 1) represents the concentration of any one of these three drugs at the location P ijk in the (t+1) time step, and λ d is the drug diffusivity. DR l ijk ðtÞ; l ¼ 1; 2; . . . ; 6 are the drug concentrations of P ijk 'snearest six neighbors in the current time step. Pe d is the vessel permeability of the drug, and U d represents the drug uptake rate. Fig 3 shows a flowchart of the algorithm.At the intracellular scale, Eqs 1-7 describe the phenotypic (migration, proliferation, quiescence or apoptosis) switches of the MMs, OBs and OCs. At the intercellular scale,the dynamics of molecules in the cellular signaling pathways after receiving cytokine stimulation from other cellsare represented by Eqs 8 and 9. In addition, MMs compete for the best location to migrate or proliferate, as defined by the distance to candidate locations and the cell density(Eqs [10][11][12][13][14]. At the tissue scale, a set of simplified reaction-diffusion equations (Eqs 15 and 16) describe the spatial concentration changes in cytokines(DKK1, OPG, RANKL, and TNF-α) and drugs(BHQ880, Lidamycin, and glucocorticoids). These changes not only remodel the tumor microenvironment but also greatly influence the behavior of MMs, OCs and OBs (cytokine secretion, proliferation, migration, or apoptosis) at the intracellular scale. The details of the model parameters are listed in Table 1.

Results
This model was implemented using Visual C++ [25] and was used to predict the response of multiple myeloma cells to drug combinationsas well as the related balance between osteoclasts and osteoblasts in response to treatment with drug combinations. Tissue scale behavior: The simulation time encompassed time steps 0 to 120 and each time step represented 2 hours.   BHQ880 was injected at time step 20 to restore the balance between OBs and OCs, when the populations of the three cell types began to differ, asshown in  Fig 4; however, this drug could not change the trend of the three curves. Fig 7 shows the 3D snapshots of the tumor system with the single-agent BHQ880 treatment at the following time steps: t = 40, t = 50 and t = 120.
The MM proliferation rate served as the target to verify the predictive power of the model after injecting BHQ880. Fig 8 shows a small difference between the simulated and experimental data [16] after 120 time steps under similar initial conditions.
GCswere employed to restore the balance by increasing the apoptosis rates of OCs and OBs. Fig 9 shows that when GCs were injected at time step 20, the populations of all three celltypes quickly decreased before time step 60. Thereafter, the MM and OC populations quickly increased, whereas the number of OBs decreased to nearly zero during the rest of the simulation. This phenomenon is also illustrated in Figs 4 and 6. Fig 10 shows the spatial information of these three cell types at time steps t = 60, t = 65 and t = 120. The MM proliferation rate served as the target to verify the predictive power of the model after the injection of GCs. Fig  11 shows a small difference between the simulated and experimental data [26] after 120 time steps undersimilar initial conditions.
Lidamycin promotes the apoptosis of MMs. Fig 12 shows that the number of MMs and OCs markedly decreased before time step 65, whereas the number of OBs significantly increased before time step 70. After time step 65, the MMand OC populationsconsiderably increased. Around time step 85, the OB and OC populations intersected. After time step 120, the number of OBs decreased to almost zero. In particular, the final number of MMs and OCs after treatment with Lidamycin was greater than that without treatment (Fig 4). Fig 13 shows the spatial information of these three cell types at time steps65, 70 and 85.TheMM proliferation rate is employed as the target to verify the predictive power of the model after GC injection. Fig 14 shows a small difference between the simulated and experimental data [18] after 120 time steps undersimilar initial conditions.  Hence, multi-drug combination treatment plans were developed. Each combination of two or three drugs was injected into the patient between time steps t = 20 andt = 60, and we observed the treatment effect until time step t = 120. Similar to previous research [17], the ratio of these drugs was 1:1. The diagram of the upper left corner in Fig 15 shows that the growth rates of both OCs and OBs were greatly inhibited by two drugs, whereas the MM population increased from time step t = 20 to t = 60.
The diagram of the upper right corner in Fig 15 demonstrates that the combination of BHQ880 and Lidamycin markedly inhibited MMs between time steps 20 and 60 while increasing the number of OBs, but this combination did not affect OCs. After time step 60, the  The diagram of the lower left corner in Fig 15 shows that the combination of GCs and Lidamycin clearly regulated MMs and OCs between time step 20 and 60 and increased the number of OBs. However, without these drugs, the number of MMs and OCs quickly increased, whereas the number of OBs decreased.
The diagram of the lower right corner in Fig 15 shows that the combination of GCs, BHQ880 and Lidamycin not only significantly regulated MMs and OCs between time step 20 and 60 but also increased the number of OBs. After discontinuing the drug (time step 60), the number of MMscontinued to decrease, and the balance between OCs and OBs was restored. Because one of the major aims of this study was to investigate which drug combination efficiently restores the initial balance between OCs and OBs and kills MMs to avoid relapse, it was necessary to develop Eq 17 to evaluate the performance of the multi-drug combinations.  Table 2) indicate that neither single drug nor twodrug combinations could restore the initial balance between OBs and OCs and avoid relapse. Nevertheless, the three-drug combination (Fig 15) showed significant potential for realizing this aim. For this reason, we ascertained which combination of these three drugs had the best performance as measured by the smallest R drug in Eq 17.

Cytokine analysis in response to combination therapy
In total, 20 doses of the three drugs were evaluated, from 0.1X to 10X in geometric sequence relative to the original dose. Then, we explored the efficacy of various combinations of the three drugs. Although conventional and high-dose chemotherapy can evoke frequent responses in patients with multiple myeloma, it remains an incurable disease, largely due to drug resistance. The bone marrow microenvironment not only promotes the survival and growth of myeloma cells but alsoparticipates in the development of resistance to single conventional and novel agents. To prevent recurrent multiple myeloma and induce complete remission, it is important to develop combination therapy strategies that inhibit certain interactions between myeloma cells and the bone marrow microenvironment [13,27]. R drug represents the drug efficacy. The minimum R drug value occurred at the doses of 10, 6, and 4 for the combination of GCs, Lidamycin and BHQ880 (Fig 16).
The Loewe additivity [28,29] evaluates whether the combination effect of BHQ880, GCs and Lidamycin is synergistic.The combination index of Loewe synergy is defined as the ratio of the total effective drug dose (combination versus single drugs) required to achieve a given effect, as stipulated in Eq 18: x ð18Þ Fig 11. A comparison of the experimental [17] and simulated data after GC treatment. where d 1 (BHQ880), d 2 (GC) and d 3 (Lidamycin) are the drug combination doses located in the combination isobologram with respect to R drug . GC ð1Þ x , GC ð2Þ x and GC ð3Þ x represent the single-agent concentrations of BHQ880, GCs or Lidamycin with respect to R drug , respectively. CI Loewe < 1, CI Loewe = 1 and CI Loewe > 1 indicate Loewe synergy, antagonism and addition, respectively. In Eq 18, R drug = 960is the threshold for evaluating the drug effect. The simulated results show that the combination of BHQ880, GCs and Lidamycin is synergistic (Fig 17).

Parameter sensitivity analysis
To evaluate the impact of the parameter values on the behavior of the multiple myeloma cancer modeling system, we analyzed the sensitivity of our model based on the following parameters: β ω , β mm T , K D , K w , λ C , DEG, λ d , Se G , Pe d , and U d We varied each parameter individually over the ranges shown in Table 1, while fixing other parameters at their base values. The ranges of the parameters were obtained from the literature, and these values are listed in Table 1. We were limited by the relatively high computing cost of the ABM, so we performed 10 simulations for each set of parameters. To assess the influence of the parameters, we calculated the Spearman rank-order correlation [30]of each parameter versus the number of total MMs, OBs, and OCs. Although high-dose chemotherapies (GCs and Lidamycin) can significantly increase the apoptosis of multiple myeloma cells, the disease inevitably relapses due to the acquisition of drug resistance. One of the reasons for drug resistance is that the bone marrow microenvironment is associated with the resistance to single-agent chemotherapy. Combination therapies will be required to increase cytotoxicity and overcome drug resistance. It has been reported that OCsinduce myeloma proliferation [31]. OB stimulation may also have beneficial effects on decreasing MMgrowth and bone destruction [32]. The OB-OC balance might play a crucial role in supporting myeloma relapse. Therefore, we consider that multi-drug combinations are necessaryto effectively inhibit MM growth and relapse. The simulated results demonstrate that two-drug combinations (Fig 15) offer a better treatment effect than a single drug (Figs 6-14), but MMs kept growing at the late stage of treatment due to the imbalance between OBs and OCs. Although two-drug combinations are not ideal for MM treatment, the results inspired us to evaluate three-drug combinations to treat MM. As shown in Fig 15, MM growth decreased, and the OB-OC balance was restored by a three-drug combination.To optimize the best ratio for the three-drug combination, Eq 16 was developed to evaluate drug efficacy. Fig 16 shows that a 1:2:4 ratioof BHQ880, GCs and Lidamycin is the best for killing MMs and restoring the balance. Lastly, Fig 17 illustrates that the three-drug combination has the strongest synergistic effect.
In general, this research details how to employ a 3D multi-scale agent-based cancer model to describe the communication among MMs, OBs and OCs to the best of our knowledge, and it also investigates the optimal use of drug combination therapy to cure multiple myeloma and restore the balance between OBs and OCs. These simulated results demonstrate that restoring the balance of OBs and OCs not only markedly inhibits the growth rate of MMs but also ensures the stable metabolism of bone tissue.