Targeting the Biophysical Properties of the Myeloma Initiating Cell Niches: A Pharmaceutical Synergism Analysis Using Multi-Scale Agent-Based Modeling

Multiple myeloma, the second most common hematological cancer, is currently incurable due to refractory disease relapse and development of multiple drug resistance. We and others recently established the biophysical model that myeloma initiating (stem) cells (MICs) trigger the stiffening of their niches via SDF-1/CXCR4 paracrine; The stiffened niches then promote the colonogenesis of MICs and protect them from drug treatment. In this work we examined in silico the pharmaceutical potential of targeting MIC niche stiffness to facilitate cytotoxic chemotherapies. We first established a multi-scale agent-based model using the Markov Chain Monte Carlo approach to recapitulate the niche stiffness centric, pro-oncogenetic positive feedback loop between MICs and myeloma-associated bone marrow stromal cells (MBMSCs), and investigated the effects of such intercellular chemo-physical communications on myeloma development. Then we used AMD3100 (to interrupt the interactions between MICs and their stroma) and Bortezomib (a recently developed novel therapeutic agent) as representative drugs to examine if the biophysical properties of myeloma niches are drugable. Results showed that our model recaptured the key experimental observation that the MBMSCs were more sensitive to SDF-1 secreted by MICs, and provided stiffer niches for these initiating cells and promoted their proliferation and drug resistance. Drug synergism analysis suggested that AMD3100 treatment undermined the capability of MICs to modulate the bone marrow microenvironment, and thus re-sensitized myeloma to Bortezomib treatments. This work is also the first attempt to virtually visualize in 3D the dynamics of the bone marrow stiffness during myeloma development. In summary, we established a multi-scale model to facilitate the translation of the niche-stiffness centric myeloma model as well as experimental observations to possible clinical applications. We concluded that targeting the biophysical properties of stem cell niches is of high clinical potential since it may re-sensitize tumor initiating cells to chemotherapies and reduce risks of cancer relapse.


Introduction
Multiple myeloma (MM) and other tumors have a small population of tumor initiating (stem) cells that retain key stem cell properties including self-renewal and tumorigenesis [1][2][3][4][5][6][7][8][9][10][11][12][13]. Recent reports [3,4] showed that a small population of CD138negative B cells with ''side population'' characteristics present in myeloma. These cells have clonogenic potential in vitro and, when engrafted into immunodeficienct/nonobese diabetes (SCID/ NOD) mice, can initiate de novo myeloma lesions of bulk of CD138+ cells in both primary and secondary transplant assays. Additionally, these myeloma initiating cells (MICs) have shown higher resistance to chemotherapeutic agents and thus are more likely to survive despite therapies [1][2][3][4][5][6][7][8][9][10]. These findings have led to the hypothesis that MICs survive chemo-and radio-therapies, regenerate the bulk of tumors, and thus cause the disease relapse. This idea is consistent with the clinical observation that disease relapse in multiple myeloma patients is common even if patients are treated with new therapeutic agents that can initially result in complete clinical responses [14][15][16]. Understanding and controlling MIC drug resistance is critical to the development of new therapies for the cure of myeloma.
Our group pioneered the research of the roles of biophysical properties in blood cancers and established the mechanism of the MIC-stroma positive feedback loop [17,18]. Previous studies on the interactions between BMSCs and myeloma cells, especially MICs, have predominantly focused on biochemical communications such as the stimuli of growth factors, cytokines and chemotactic via paracrine signaling [19]. However, recent studies in solid tumors have indicated that a critical stage of the malignant transformation journey of cancer cells involves marked alterations in the biomechanical phenotype of the cell and its surrounding microenvironment [20,21]. Indeed, it has been suggested that targeting the microenvironments (the ''niches'') of the tumor stem cell could result in a reduction of the tumor burden [22][23][24]. Bone marrow stromal cells (BMSCs), one of the major cellular components in the MIC niches, are in close contact with MICs, and the biomechanical properties of BMSCs, besides chemical communications, also influence the local microenvironment of MICs and hence MIC fates. We have recently demonstrated that Myeloma-associated BMSCs (MBMSCs) from patients are much ''stiffer'' (higher Young's modulus level) and more contractile than Normal BMSCs (NBMSCs). Hydrogels are widely used to mimic the in vivo cellular microenvironments [25,26], so we have utilized hydrogels of various stiffness levels to investigate the impact of such biophysical property on MIC-driven myeloma development. We have shown that stiffer hydrogels support colony formation and adherence of MICs better than softer hydrogels, suggesting that myeloma BMSCs provide myeloma cell-friendly microenvironments for MICs via exerting biomechanical forces [17,18]. We also have demonstrated that MICs over-secrete SDF1 than mature myeloma cells and that treatment of CXCR4 inhibitor, AMD3100, leads to decreased adherence of MICs to MBMSCs, undermined colony formation potential of MICs, and better in vitro and in vivo drug efficacy of Bortezomib(BZM). These discoveries were also consistent with other reports [27].
With the perspective that the biophysical properties of MIC niches in bone marrow may be a promising drug candidate, there is an urgent need to develop mathematical models and tools to estimate drug effects on MICs within the context of their niches during when screening and evaluating drug candidates. Currently pharmaco-industries use the killing efficiency of the bulk tumor as the major in vitro drug screening criterion. Such practices overlook the importance of MICs cell microenvironments contributing to the common failure of translating promising drug candidates discovered in screening into clinical usage. However, mathematical models that involve the cancer initiating (stem) cells and their niches are still rare, except for few examples [28][29][30] including our recent work [31].
Modeling the dynamics of the MIC-derived myeloma lineage is crucial for associating the abnormality at cellular level with the features of the myeloma pathogenesis and pathophysiology at tissue level. We recently developed a myeloma lineage model [31] to systematically analyze the myeloma lineage development and experimentally monitor the major sub-populations in the lineage. Briefly, cells involved in the myeloma lineage were classified as myeloma initiating cells (MICs), myeloma progenitor cells (PCs), and matured myeloma cells (MMs). These three sub-populations were experimentally distinguishable by dual-staining for the expression of plasma cell surface markers CD138 and for side population (SP) using Hoechst staining, each showed unique features related with lineage sub-populations. MICs, recognized by SP staining and negative CD138 staining (SP/CD1382), showed unlimited renewal capability, potential to initiate and fully reestablish the whole myeloma lineage, and enhanced drug resistance. PCs, defined as the non-SP and negative CD138 expression (non-SP/CD1382) population, were only able to passage for limited rounds before beginning to express CD138 and losing the key PC character. The CD138+ sub-population, which was of 90% to 95% of the total myeloma population, was classified as the MMs and exhibit very limited proliferation capability. In this work we used similar definitions except further distinguishing terminal MMs (TMMs) that could no longer proliferate from those that still could divide (MMs).
To comprehensively illustrate the interaction between the myeloma initiating cells and their niches and to develop effective drug treatment, the following key questions need special attentions: (1) How do BMSCs, MICs, myeloma progenitor cells (PC), matured myeloma cells (MMs) and terminal MMs (TMMs) communicate in MIC niches? (2) How do the biomechanical phenotypes of BMSC, in terms of cell stiffness and contractibility, modulate MIC's growth and fates? (3) How do MIC growth and differentiation drive the development of myeloma? (4) How do the molecular level intracellular features of MICs and their stromal counterpart, BMSCs, contribute to the tissue level 3D cancer growth via the intercellular cell-to-cell interactions? And (5) How does the interaction between tumor cells and their niches change the drug treatment effect?
To address these questions, in this study we developed a 3D multi-scale agent-based model (ABM) using Markov Chain Monte Carlo approaches [32] to study the role of tumor-stroma interactions in multiple myeloma tumor progress. The system was classified into three levels: the intracellular, the intercellular, and the tissue level. Their relations were conceptually defined as ''interfaces'' among these levels. The agent-based model integrates events of multiple spatial and temporal scales. (1) Spatial scales: Intracellular signaling pathways were encapsulated into each cell type to determine the BMSC intercellular biomechanical phenotype (cell stiffness) or tumor cell behaviors (migration or proliferation). Cancer cells competed for the best location in 3D extracellular matrix to migrate or proliferate regarding to the change of BMSC cell's stiffness and sensitivity of chemoattractant (SDF1). In turn, chemoattractant cues at the tissue level triggered intracellular signaling pathways inside cells via the agent interfaces (receptors), and the resultant feedbacks were the changes of either cells' properties (change of cell stiffness or sensitivity to outside stimuli) or behaviors (secretion of chemoattractant, proliferation, differentiation, apoptosis, or migration). (2) Temporal scales: The model covered minute-to-hour-level signaling dynamics; day-level cell division, apoptosis, and local migration; week-level drug responses; and month-level tumor growth.
The spatial characteristics of the 3D myeloma cell distributions was described by the local cell metrics [33] we previously developed and the corresponding parameters estimated using the 3D cell co-culture system levitated by magnetic field and nanoshuttles [34].

Model development
We established an agent-based multi-scale model to simulate the development of myeloma in various bone marrow microenvironments in three-dimensional space, and validated the model with experimental data. The central hypotheses are: 1) MICs drive the development of multiple myeloma; and 2) the SDF1-centric biophysical and chemical positive feedback loop boosts MIC growth and colonogenesis, and protects MICs from drug treatments. As shown in Figure 1, myeloma development was simulated at intracellular, intercellular, and tissue levels. In this model we included five types of cells which were represented by five types of agents of encapsulated intracellular signaling events and interfaces through which these agents communicated with their microenvironments: (1) Bone marrow stromal cell (BMSC) agent. This agent modeled the 3D reticular network formed by BMSCs. BMSCs changed its biophysical stiffness by cell contraction in response to SDF-1 hints in microenvironment and thus altered the threedimensional distribution of bone marrow stiffness; (2) Myeloma initiating cell (MIC) agent: This agent represented myeloma stem cells enriched by side population staining [4,35]. MIC sensed the stiffness of its microenvironment and accordingly modulated its proliferation, differentiation, apoptosis, and drug resistance; Please note that although mathematically PC, MM, and TMM can be lumped into one agent, each of them have unique biological behaviors (such as apoptosis rates) and can be experimentally measured by flow cytometry. Therefore we kept them distinguished for better understanding the simulation results for future experimental validations, and for clinical applications.
At the intercellular level the communication among agents were simulated by the biomechanical model representing the MIC-BMSC positive feedback loops and the myeloma lineage model illustrating the dynamics of different myeloma populations. Interactions between BMSC and MIC were: (a) MICs secreted SDF1 into the neighborhood extracellular matrix. The diffusion of SDF1 in the tissue defined the biochemical microenvironment in bone marrow. The BMSCs contracted according the paracrine SDF1 in their neighborhood and altered the biophysical properties of the bone marrow. MICs sensed and preferably attached to stiffer BMSCs. When attached, stiffer BMSCs boosted proliferation and self-renewal of MICs, and thus promoted the expansion of both MIC and myeloma, and drove sustainable multiple myeloma cancer growth. Stiffer BMSCs also protected MICs from the treatment of cytotoxic drugs.
At the tissue level, the SDF1 profile defined by MIC paracrine and diffusion and the drug concentration defined by the administration of chemotherapeutics determined the biochemical microenvironment in bone marrow, the tissue stiffness defined by BMSC contraction determined the biophysical microenvironment, and myeloma proliferation and migration determined the cellular distribution in bone marrow.
Multiple myeloma cells at different stage of differentiation as defined above were initially seeded at the center of the bone marrow, 100 cells for each cell type, to mimic the initiation stage during myeloma spreading to new locations or after implantation into the bone marrow of animal models. Bortezomib (BZM) and AMD 3100 were used as the representatives of cytotoxic drugs and inhibitors against MIC-BMSC interactions, respectively. After 200 hours of tumor seeding these two drugs were delivered intravenously into the simulated extracellular matrix at combinations of various doses to eliminate multiple myeloma cancer cells. The chemotherapy was scheduled that each cycle included a 400hour treatment period followed by 500 hours rest (about 2 week treatment followed by 3 week rest), which was typical in clinics. In our experiments usually the animals were sacrificed and the outcome was examined 200 to 400 hours after the first treatment period or one month after tumor xenografting, so the simulation results of the first 600 to 1,000 hours (about one month) were shown. We only simulated the tumor responses to the first treatment cycle because multiple myeloma quickly develops drug resistance and therefore elimination of all MICs in the first treatment cycle is crucial. Parameters were directly determined or indirectly inferred based on our previous experimental results  [17,18,36], or set according to the best of our knowledge of multiple myeloma, and were summarized in Table S1 and Table  S2 in Supporting Information.
To reduce the computational cost the diffusion of the two drugs in bone marrow was assumed to be instantly. This assumption is valid since the myeloma-enriched regions in bone marrow are well vascularized and the diffusion of small molecule drugs is significantly faster comparing with either typical cell behaviors such as migration, proliferation, and apoptosis, or comparing with the diffusion of SDF-1.
Twelve dose levels (including level 0) for each drug and their full combinations were examined under both myeloma-associated BMSC (MBMSC) and normal BMSC (NBMSC) microenvironments, each condition simulated for 20 times. Time resolution (time step) was 2 hours. Totally 288 conditions was explored by 5,760 simulations using parallel computation on a Dell TM PowerEdge TM 850 server of 48 CPU cores.
Examples of model simulation results of MIC-driven myeloma growth in bone marrow were visualized and compared with in vitro 3D levitated culture image in  Figure 2 G. Stiffness was labeled in color, from blue to red with respect to the increase of stiffness. Since in vivo 3D imaging of BMSC contraction is yet unavailable, we co-cultured myeloma with BMSCs in vitro using the 3D levitation system to mimic the bone marrow microenvironment, and 4-day co-cultured tumor tissue was stained and 3D imaged by confocal microscopy (details see Methods section) and shown in Figure 2 H. BMSCs (white arrows) and myeloma cells (green arrows) were recognized by cell shapes according to F-actin staining (red) as well as nuclei staining (blue). The primed MBMSCs (yellow arrows) were recognized by the formation of stress fibers. The size of the experimentally cultured tumor at day 4 was about 2666492 mm, which was very close to the diameter of the infiltration frontier (about 460 mm) of the simulated results at time point 100 hr (Figure 2 A).
Details of the model can be found in the Methods section. Typical movies as Movie S1, S2, and S3 can be found in Supporting Information, and the corresponding simulation conditions in the caption of Figure 2.

Model validation
The simulation results of the agent-based model of myeloma growth under different biophysical microenvironments were compared with experimental results [18,36] to validate the model. The simulated MIC populations as well as the total tumor sizes under the NMBSC and MBMSC (Blue vs. red lines in Figure 3 A and B) after 4 weeks growth were consistent with experimental results [18] (Blue vs. red bars in Figure 3 A and B).

Effects of MBMSCs in disease development
We first tuned the model to recapture experimental observations and explored the roles of myeloma-associated MBMSCs in disease development using the NBMSCs as controls. The simulation results of tumor growth and the comparison with experimental data were summarized in Figure 4 and

Effects of MBMSCs on chemotherapy outcomes using BZM
To explore the drug resistance to chemotherapies and the risk of tumor relapse in MBMSC and NBMSC microenvironments, we simulated the first BZM treatment cycle at 11 dose levels covering two magnitudes. BZM was delivered 200 hours after tumor initiation and the treatments lasted for 400 hours. The dynamic of total tumor cells as well as MICs during the first 1000 hours were shown in Figure 5. After drug delivery, the tumor populations quickly dropped for both MBMSC and NBMSC (red and blue lines in Figure 5 a, respectively) cases. The drug efficacy, in terms of the elimination of myeloma population, was similar in both cases, and during the first treatment cycle the tumor cell number decreased to a very low level that were clinically not detectable. Two typical simulations were shown in Figure 5, and the rest were discussed in the following sections.

Inhibition of MIC-BMSC communications by AMD 3100 treatment
We then examined if the inhibition of MIC-BMSC communications could deprive MICs from the protection of their stroma and re-sensitize the tumor in the MBMSC context to chemotherapy. AMD 3100, a competitive inhibitor of SDF-1, was delivered at 11 levels crossing two magnitudes, either alone or with BZM, to myeloma growing with either myeloma-associated or normal BMSC stroma. The effects of medium AMD dose were visualized by two simulations in Figure 6 and the rest summarized in Figure 7. AMD-treated and AMD-free cases were denoted in solid and dashed lines, and the myeloma-associated and normal BMSC niches in red and blue colors, respectively. Consistent with

Combined drug effects on both cases
To evaluate if interrupting the MIC-BMSC communications has clinical potentials for multiple myeloma, we tested the combinatorial effects of the two representative drugs, AMD3100 and BZM. For each drug, 12 doses (0 as control and 11 levels from 0.16 to 106 in geometric sequence relative to the original dose) were selected and the full combinations of the two drugs were explored for drug efficacies. The default doses for each, denoted as 16, was defined as the minimum doses that, for BZM, eliminated all MICs in NBMSC case and, for AMD3100, deprived MBMSC support. Each combinatorial condition was simulated for 20 times at 400 time points for 5 cell types and SDF-1 concentration distributions, and totally 13,824,000 data sets were generated. The means of replicates were used for further analysis, and the synergistic effect maps were further smoothed with a two dimensional Gaussian kernel (h = 2) for robust visualization.
Loewe combination index [37,38] (for an introduction and a guide of interpretation of Loewe combination index please refer to Figure S1 in the Supporting Information), defined according to the 50% MIC population reduction by the end of the two week chemotherapy (E50 at time point 600 hr), was utilized to evaluate the synergism of AMD3100 and BZM while treating myeloma in either NBMSC or MBMSC niche. We selected the total myeloma death rate as well as MIC death rate as the indicator of drug efficacy because as demonstrated in Figure 5 (first and second rows, respectively). MIC played essential role in post-treatment disease relapse [4]. The SDF-1/CXCR4 inhibitor showed on moderate effects on the efficacy of BZM when treating multiple myeloma developed in both normal and myeloma-associated bone marrow stromal cell microenvironments in terms of both the shrinkage of total tumor size (red contour lines in Figure 7 a and b, respectively) and the decrease of MIC populations (red contour lines in Figure 7 c and d, respectively). In contrast, dramatic synergistic effects of the two drugs were observed if the total elimination of MIC (E100) was used as the criterion for Loewe combination index in MBMSC cases (green contour lines in

Discussion
This work previewed the druggability of the biophysical features of cancer stem cell niches, and cast new light on the strategies to overcome the drug resistance and relapse of multiple myeloma.
The central hypothesis of this work, that the positive feedback loop between MIC and MBMSC via SDF-1 paracrine and the increase of MBMSC niche stiffness promotes myeloma development and is responsible for drug resistance and cancer relapse, has been successfully realized in this multi-scale agent-based model. Figure 2 showed that MBMSCs close to MICs were primed and activated by SDF-1 secreted by MICs provided a stiffer microenvironment; such pro-ontogenetic niches in turn boosted MICs in terms of proliferation and drug resistance.
It was the first time the dynamics of the bone marrow stiffness during the development of myeloma was simulated and visualized in three-dimensional space at multiple scales. While the in vivo assay technologies of the bone marrow biophysical properties are yet not available, this model provides insights into the mutual communications between cellular and mechanical information. The simulations (Figure 2 and Supporting Information Movie S1 and S2) suggested that both the infiltration process of myeloma as well as the distribution of bone marrow stiffness were highly dynamic and the two types of distributions showed temporal and spatial lags. In vivo monitoring the changes of bone marrow biophysical properties can be a valuable method to directly estimate drug efficacy of myeloma-associated BMSC inhibitors.
Our simulation results also emphasized the crucial role of MICs in disease relapse, which was illustrated in Figure 5. The elimination of cancer stem cell is the essential goal of chemotherapies for multiple myeloma; otherwise the survived MICs will drive the re-initiation of the cancer and the relapse of the disease [4]. However, due to the drug resistance of MICs and the protection of MIC niches, it is challenge for traditional chemotherapies to target and purge all MICs [4,39]. Targeting the MIC niches thus becomes an intriguing strategy.
Besides that the model is conceptually consistent with our knowledge of multiple myeloma, we also confirmed that the simulations reproduced the key experimental findings of MIC/ MBMSC interactions [17,18,36]. Although parameters of this piloting model were roughly and bona fide determined from literatures, our experiments, and the best of our biological knowledge (Table S1 and Table S2 in Supporting Information), the model is capable to re-capture our experimental observations ( Figure 2 H and Figure 3).
(1) Tumor initiation. One major feature of multiple myeloma is its frequent metastasis. Understanding the initiation and early stage development of myeloma at new bone marrow sites is thus critical for clinical intervention. However, it is yet challenging to monitor such events in vivo. We established the 3D levitated co-culture system to mimic the early events of the development of secondary myeloma and to reveal its unique features. The similar growth trends from the simulated results and the 3D co-culture in terms of tumor size (Figure 2 A vs. H) suggested that the model is capable of predicting growth trends after successful metastasis as well as the post-treatment re-initiation of tumors. The geometric shapes of myeloma growth are dynamic and subject to the stochastic feature of cell migration as well as the context of cell growth, though. (2) Pro-oncogenetic MBMSC microenvironments. As mentioned before, accumulating evident suggest that MICs were capable to prime MBMSCs but not NBMSCs in terms of BMSC stiffness, which in turn promotes myeloma growth and enhance drug resistance. The model re-captured these features. MICs population was boosted in myeloma-associated  [4,39], and such drug resistance was enhanced at the presence of MBMSCs ( Figure 5 red vs. blue lines) [18].
The simulation results were consistent with our and others experimental observations, which provided a solid basis for the exploration of the synergistic effects of niche-specific drugs with cytotoxic chemotherapeutics. Please notice that although the in vitro 3D levitated culture has been known to be able to mimic in vivo microenvironment [34], the bone marrow contexts are still significantly different from the in vitro environments. Therefore, the highly consistent results between experimental data and simulation (Figure 2 A vs. H) only suggested that the simulation results were reasonable and the model is sufficient for piloting in silico study. In vivo assays of myeloma growth and drug responses are still necessary for parameter estimation, model training, and accurate predictions, which is future work and is beyond the focus of this study.
Loewe drug combination analysis [40,41] strongly suggested that inhibitors that block the MIC/BMSC communications are of high clinical potentials. Traditional Loewe combination analyses for chemotherapeutics often focus on the bulk tumor and use isobole at effect level 50 (E50) as standard index, which is convenient to monitor in animal models and clinical chemotherapy evaluations using micro-CT, bioluminescence, biopsy, and others. The elimination of the bulk tumor is also used as clinical index. However, this standard may not accurately reflect the drug efficacies against tumor-initiating-cell-driven tumor development. Capability to eliminate of all tumor initiating cells should be the key evaluation index. We thus defined the E100 isobole against MICs in our Loewe analysis and explored drug combinations that could successfully kill all MICs in the first treatment cycle. Simulations of the combination effects of the cytotoxic drug (BZM) and the niche-specific inhibitor (AMD3100) in myeloma-associated stroma, evaluated by the E100 isoboles against MICs, suggested strong synergistic effects (green lines in Figure 7 D). In contrast, E50 isoboles only indicated slight synergism between the two drugs (red lines in Figure 7 A through D). The control groups of the normal BMSC environments showed no responses to niche inhibitors (Figure 7 A and C), which further confirmed that AMD3100 did not directly affect MICs and the strong synergism observed in the MBMSC cases were due to the interruption of the MIC-BMSC communication. Such dramatic difference suggested that for TIC-driven tumors drug effects should be estimated in terms of the total eradication of TICs, and E100 isoboles against  [18] and the total tumor sizes (B) [18] under NBMSC or MBMSC environments as well as the dose responses in terms of total tumor size decrease (C) [36] were compared with our published experiment results. doi:10.1371/journal.pone.0085059.g003 TICs provide a more sensitive evaluation of the effectiveness of treatments. Specifically, for the multiple myeloma treatment evaluation and relapse risk estimation, our results highly suggested using the percentage of MIC population in patient bone marrow aspirate by either side population staining or surface marker immuno-staining (e.g. CD138) and flow cytometry.
Simulations also suggested that niche-specific drugs are strongly synergistic with tumor-specific cytotoxic drugs and re-sensitize myeloma to the otherwise resistant drugs. Cytotoxic drugs alone, though may efficiently kill the bulk of tumor, sometimes are not sufficient to demolish all TICs. As shown in Figure 7 D, without the aid of AMD3100, BZM alone cannot extinguish all MICs and thus the disease will relapse soon after the chemotherapy ( Figure 5). Taken together, the physical properties of the myeloma initiating cell niches is highly targetable, and the inhibitors of the interactions between cancer initiating cells and their stroma are promising co-drugs for traditional cytotoxic agents.
As a piloting research, this computational model is conceptual and some parameters were not directly estimated from experiment data, and some in vivo/clinical data and assays such as drug delivery and metabolism in patients may improve the accuracy of the predictions. However, we have omitted such details due to the following reasons: 1) To keep the focus of this work of the preliminarily modeling of a novel system. It is the first time to simulate the druggability of the biophysical properties, mainly the stiffness, of the bone marrow microenvironments and the major focus of this work is to establish a reliable model according to the  relatively solid in vitro and animal experiment data. At this preliminary stage, it may be not ready to reach a precise predictive model for clinical treatment outcomes by including clinical data about the bone marrow properties as well as the underlying mechanisms (such as the diffusion processes of cytokines and the pharmacokinetics and the pharmacodynamics of bortezomib and AMD3100). 2) To address our major questions and reach conclusions with a reasonably simple model as possible. We have also simplified the detailed temporal profiles of drugs during chemotherapies and adjuvant therapies to constant levels as such detailed dynamic features, though important for providing optimal dosing in clinical treatments in the future, are not our focus at this stage. Indeed such simplification approaches are common in preliminary models of new systems and do not change the main conclusions of this work. 3) Furthermore, detailed molecular mechanisms of underlying signaling pathways are still missing before we can reach an accurate and more predictive model. The major goal of this work is to extend our current knowledge of how MICs modulate the physical properties of their niches to evaluating the therapeutic potential of targeting such remodeling process. For such purpose, the results of this work suggest using the complete eradication of TICs as drug efficacy index for TICdriven cancers, emphasize the importance of targeting TIC niches, and encourage following-up investigations of niche-specific therapeutics.

Conclusions
It was the first time the three dimensional modeling of the biophysical properties of cancer stem cell niches and the therapeutic potential of targeting such mechanical features were attempted at the system level. Our in silico model successfully realized the bio-model of the MIC-BMSC chemical-physical positive feedback loop, and was capable to capture the key features of our experimental observations. Simulations suggest that the isoboles of drug capability to completely eliminate TICs in Loewe combination analysis can better predict treatment outcomes. Drug synergism analysis suggested that intervening the communications between myeloma initiating cells and their niches dramatically enhanced the efficacy of cytotoxic drugs against myeloma initiating cells, re-sensitized multiple myeloma to chemotherapies, and reduced risks of cancer relapse.

Methods
We defined five types of agents in the model to represent BMSC, MIC, PC, MM and TMM. We initialized the bone marrow as a cylinder 3D rectangular framework, with the tumor extracellular matrix (ECM) filled with BMSC agents distributed evenly across the 3D ECM and the mixed MIC, PC, MM and TMM agents at the center of the bone marrow as a sphere. Such initiating cell ratios were according to the myeloma initiating cell hypothesis as well as our previous publications that the MICs can restore the whole myeloma population, including the composition ratios of myeloma cell types. This multi-scale modeling consisted of three scales: intracellular, intercellular and tissue scales, which were illustrated in Figure 1, Figure 8, and described in details in the following sections. The model was generally presented as pseudo-code and listed in Figure S6. Detailed flowcharts of each myeloma agents were illustrated in Figure S2 (for MIC agent), Figure S3 (for PC agent), Figure S4 (for MM agent), and Figure S5 (for TMM agent) in the Supporting Information. Parameters were determined according to our previous work as well as literature [17,18,36]. Mathematical and computational details are elaborated in Support Material.

Stochastic Simulation of Cell Behaviors
The Markov Chain Monte Carlo approach was used to simulate cell behaviors of each individual cell. As shown in Figure 8, cell behaviors were simulated by probability-based rule implementation. A cell sensed the hints of its neighborhood such as stiffness, SDF-1 level, and drug doses, processed them with imbedded pathways, and outputted the changes of probabilities of corresponding cell behaviors including cell proliferation rate, apoptosis rate, differentiation ratio, migration rate, contraction rate, and cytokine secretion rate. Cell decision was then determined by rolling a dice and compared with the probability of a given cell behavior. Cell behaviors in turn remodeled the properties of its niches. Details of each cell behavior for each type of cell agent as well as the corresponding rule was discussed in the following sections as well as the Supporting Information.

Intracellular Scale
Each BMSC agent encapsulated signaling pathway functions to determine its biomechanical phenotype switch, using an agentspecific Hill function [42] to describe the SDF-1/CXCR4 signaling pathway which regulated the local stiffness in response to the in-situ relative SDF1 concentration. The effects of the SDF1/CXCR4 inhibitor AMD3100 was simulated by competing with SDF1 for CXCR4 and thus reducing the effective SDF1 concentration.
Meanwhile, the responses of tumor agents (MIC, PC, MM, and TMM) to local ECM stiffness in terms of the possibilities of cells to enter the proliferation, apoptosis, and migration status were also simulated using Hill functions in similar ways. Cell decisionmaking process was defined by agent rules with such probabilities as the major inputs. The stochastic feature of the decision of an individual cell was realized by die casting simulation.
Stem cell fate determination. Once a MIC decided to enter cell cycle, its fates were further determined according to its micronenvironment. A MIC either generate two MIC daughter cells, known as self-renewal, or to two PC cells, know as differentiation, or to one MIC and one PC, known as asymmetric division. The probability of each fate was determined using Hill functions, and the decision of each MIC agent was also realized by die casting simulation as mentioned above.
Proliferation fates of PC and MM agents. The fates of intermediate cell agents were determined by their passage ages as well as the probabilities of proliferation which represent the effects of cell neighborhood characters such as stiffness and cytokine concentration, and the current cell cycle status of the cells. PC agents were different from MM agents for their limited selfrenewal capability. When maximum renewal limit reached, a PC agent differentiated to an MM agent. TMM agents did not proliferate. According to the myeloma initiating cell hypothesis, only MICs can self-renew and proliferate without limits, so defined is the LGN, which is the maximum passage number a PC cell can self-renew, or a MM cell can proliferate.

Intercellular Scale
Once a myeloma cell (MIC, PC, MM, or TMM) had determined its responses to the biomechanical properties of its microenviroment, it would proliferate, migrate, being quiescent, or undergo death, which were described in the intercellular scale.
Migration. A non-M-phase cell would migrate if it could find free space nearby. The migration was governed by rules representing space availability, stroma preference, migration speed, and stochastic effects using Hill functions and die-casting simulation as mentioned above. Division of MIC, PC and MM agents. If M-phase cell found at least one free location within searching distance, it would divide, following the same migration rules but with a much smaller migration distance (i.e., slower migration speed) so the de novo daughter cells were always next to the parental cells. If no space was available, cells would remain in M-phase and try in next round.
Apoptosis. The decision-making of cell apoptosis was simulated using a pre-defined threshold for the apoptosis rate reflecting cell microenvironment, especially the local drug (BZM and AMD) concentrations and stiffness. The whole apoptosis process took about 20 hours. RPMI 8226 myeloma cell apoptosis rates under drug treatments were determined according to our previous study [36].

Tissue scales
In the tissue scale of this ABM, the secretion of SDF1 from MIC agents and the diffusion of SDF1 in the 3D ECM defined the dynamic 3D distribution of SDF1 concentration [43][44][45]. SDF1 was uniformly initialized at the start with Dirichlet boundary.
In summary, intracellular signaling pathways were encapsulated into each cell to determine either the BMSC intercellular biomechanical phenotype (cell stiffness) or tumor cells' (MIC, PC and MM) behaviors (migration, differentiation, proliferation, or apoptosis). Cancer cells competed for the best location in 3D extracellular matrix to migrate or proliferate regarding to the change of BMSC cell's stiffness and cell density. In turn, chemoattractant cues (SDF1) at the tissue level triggered BMSC cell's intracellular signaling pathways by receptors. And the resultant feedbacks were the changes of either cancer cells' properties (change of BMSC cells' stiffness) or behaviors of cancer cells (secretion of cytokines, proliferation, differentiation, apoptosis, or migration). Thus, the 3D dynamics of bone marrow stiffness and tumor growth were simulated at multiple temporal and spatial scales.

Myeloma cell model
We chose RPMI 8226, one of the most widely used human myeloma cell lines, as the myeloma cell model because it had been shown representative of key myeloma pathogenesis features of interest in this work. In vivo and in vitro evidence from our [17,18,36] and other research groups [3,4,46] suggested the clinical translational values [46] of the MICs from RPMI 8226, which re-capture major clinical myeloma stem cell characters such as drug resistance and initiating myeloma in bone marrow by generating the whole myeloma population, both of which are boosted by the positive feedback loops with their associated stroma cells via the SDF-1/stiffness chemo-physical interactions. Specifically, RPMI 8226 has been shown to be representative among other frequently used human myeloma cell lines in terms of Bortezomib treatment responses [47,48].

Cell isolation and culture
Myeloma-associated and normal stroma cells were isolated and expanded from myeloma patients similarly as previously described [18], depending if they were diagnosed of bone marrow involved myeloma or not. Briefly, after lysis of red blood cells (RBC Lysis Solution, QIAGEN), patient bone marrow aspirate was cultured in MesenPRO RS TM medium (Invitrogen) in 5% CO 2 humidified atmosphere at 37uC for one week and unattached cells discarded. Attached MBMSC was expended under same culture conditions. Usage of these samples has been approved by the Institutional Review Board of The Methodist Hospital Research Institute (TMHRI). RPMI8226 multiple myeloma cell line was purchased from ATCC and cultured in RPMI 1640 medium (Mediatech) supplemented with 8% fetal bovine serum (Invitrogen) 100 units/ ml penicillinand 100 mg/ml streptomycin (Life Technologies) at 37uC in 5% CO2. This study was granted for Consent waiver by IRB.

Staining and imaging
At day 4 the assembled MIC/MBMSC co-cultured tissue was collected, rinsed in PBS, fixed with 3.7% paraformaldehyde (Fisher) for 10 min, blocked with 1% bovine serum albumin (Fisher), and stained with Alexa FluorH 594 phalloidin (Invitrogen) using 1:40 dilution for 30 min, followed by counter staining using 300 nM DAPI (Invitrogen) for 5 min, rinsed with DI water, and mounted with ProLongH Gold Antifade Reagent (Invitrogen) for imaging. The process was under the protection of the Magnetic Drive to avoid lost of cells. The Nikon A1 Confocal Imaging System was used to image samples for DAPI and Alexa FluorH 594 signals with optical slicing distance of 8 mM and the 3D images were reconstructed using NIS-Elements Microscope Imaging Software (Nikon). (ZIP) Figure S4 The flowchart of the MM agent.