Predicting the impact of combined therapies on myeloma cell growth using a hybrid multi-scale agent-based model

Multiple myeloma is a malignant still incurable plasma cell disorder. This is due to refractory disease relapse, immune impairment, and development of multi-drug resistance. The growth of malignant plasma cells is dependent on the bone marrow (BM) microenvironment and evasion of the host's anti-tumor immune response. Hence, we hypothesized that targeting tumor-stromal cell interaction and endogenous immune system in BM will potentially improve the response of multiple myeloma (MM). Therefore, we proposed a computational simulation of the myeloma development in the complicated microenvironment which includes immune cell components and bone marrow stromal cells and predicted the effects of combined treatment with multi-drugs on myeloma cell growth. We constructed a hybrid multi-scale agent-based model (HABM) that combines an ODE system and Agent-based model (ABM). The ODEs was used for modeling the dynamic changes of intracellular signal transductions and ABM for modeling the cell-cell interactions between stromal cells, tumor, and immune components in the BM. This model simulated myeloma growth in the bone marrow microenvironment and revealed the important role of immune system in this process. The predicted outcomes were consistent with the experimental observations from previous studies. Moreover, we applied this model to predict the treatment effects of three key therapeutic drugs used for MM, and found that the combination of these three drugs potentially suppress the growth of myeloma cells and reactivate the immune response. In summary, the proposed model may serve as a novel computational platform for simulating the formation of MM and evaluating the treatment response of MM to multiple drugs.


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 10, cell behaviors were simulated by probabilitybased rule implementation. A cell sensed the hints of its neighborhood such as stiffness, SDF-1 and TGFβ level, and drug doses, processed them with signaling pathways (via OED systems or Hill functions), and outputted the changes of probabilities of corresponding cell behaviors including cell proliferation rate, survival (apoptosis) rate, differentiation ratio, migration 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 updated the properties of its niches. Details of each cell behavior for each type of cell agent as well as the corresponding rule was introduced in the following sections.

Intracellular scale
Each BMSC agent determines its biomechanical phenotype via the ODE system of SDF-1/CXCR4 signaling pathway (see Materials and Methods). Each grid point of the BMSC represented a section of the BMSC reticular. Each BMSC agent senses the local concentration of SDF-1, then activates the SDF-1/CXCR4 pathways, and finally output the stiffness.
Similarly, the oncogenic effect of the niche stiffness on the proliferation rates of myeloma initiating cells (MIC) and mature multiple myeloma (MM) were described with a Hill function (Eq.1).
Where , l prol P l = 1, 2were the proliferation rates of MIC, and MM, respectively. 0 l P denotes the initial proliferation rate, l pathway P is the maximal proliferation rate of the embedded signaling pathway in BMSC cells. K E , E,and E 0 were the maximal stiffness, local stiffness, and base stiffness value. We assumed that the contraction of a BMSC section altered the stiffness of its local extracellular matrix, and thus defined the stiffness a myeloma cell sensed E, as the maximal stiffness of its neighbor BMSC sections within 3 grid distance.
The proliferation rate of CD8 + T cell was not only regulated by local concentration of TGFβ and the dose of Lenalidomide, but also was affected by the population of Tregs (Eq. 2).
Where 0 8 CD Prol + is the base proliferation rate of CD8 + T cells. PL denotes the relative LEN dose and K L the sensitivity of the CD8 + T cells to LEN. T ijk denotes the relative TGFβ concentration at the gird point (i, j, k); and K T denotes the threshold TGF level to suppress CD8 + T cells. The effect of Tregs on CD8 + was described with R, which determined by the population ratio of Treg against CD8 + . Based on the Haart's experimental data [1], we inferred the formulation of function R according to different ratios of Treg to CD8 + (Eq. 3).
According to formula (3)(4), we can clearly know that increase of Treg population will suppress the proliferation of CD8 + .
Finally, the proliferation rate of Treg cell was defined as Eq. 5.
Prol is the base proliferation rate of regulatory T cells (Treg). Formula (5) indicates that the production of Tregs can be promoted by TGFβ and suppressed by IMiD drug (such as, LEN).
Cell decision-making process was defined by agent rules, and the stochastic feature of the decision of an individual cell was realized by die casting simulation. Here, we used the cell decision-making of entering cell cycle as an example to elaborate the algorithm. Cell proliferation rules were defined according to proliferation rates as probabilities of two cell decisions with a die If the die fell into the interval [0, prol P ], the cell entered cell cycle and started to proliferate; otherwise, stayed quiescent. Four cell cycle phases [2,3] (G0/G1, S, G2, and M) were defined. Cells kept migrating during the first three phases, while tried to find a location to divide after entering the M phase, which was defined in the following sections. The duration of stage G0/G1 was determined by the Monte Carlo simulation (die casting) process, the S phase lasted for 20 hours (10 time steps), and the G2 phase would not exceed 10 hours (5 time steps).

Stem cell fate determination
Once a MIC decided to enter cell cycle, its fates were further determined according to its microenvironment. A MIC either generates two MIC daughter cells, known as self-renewal, or to two MM cells, known as differentiation. The probability of each fate was denoted as r r . Rules for MIC fate determination were represented in Eq.7: Where ,0 r r r r r = + ∆ . Δr describes the effects of stiffness-trigged MIC signaling pathway on cell fates, and was determined by Eq.8: Where denotes the signaling pathway activation threshold.

Intercellular scale
Each type of cells would proliferate, migrate, being quiescent, or undergo death, which were described in the intercellular scale.

Migration
A non-M-phase cell at position p would migrate if it would find free space nearby. Except BMSC agents, other four types of cell agents can migrate from one place to another free space. MICs and MMs prefer to migrate towards the stiffer BMSCs. Once they attached on BMSC sections, they prefer to stay there rather than move to non-BMSC position (Supplementary Figure S2). For the CD8 + T cells, they tried to migrate towards Myeloma cells (MICs and MMs) and to induce the lysis of tumor cells. As to Tregs, they prefer to migrate to CD8 + T cell population, and to induce the cell cycle arrest or apoptosis of the effector cells. Now, we firstly introduce the migration of MIC. In this study, if a MIC cell will migrate depends on what's its current position. If a MIC cell is not attaching on BMSC, it will definitely implement the procedure of migration; otherwise, we firstly calculate the adhesion rate (P adh ) of MIC with local stiffness via ODE system, then make dice rolling to generate a random value md and to determine if it will migrate to another place (P adh < md). According to Supplementary Figure S2, we can find that MIC prefers to continuously stay on BMSC rather than move to non-BMSC position. The migration was defined in the following rules. All the unoccupied positions within a radius r max (Eq. 10) ( ) Were listed as candidates, where D was the basic migration speed index, Δt is the time step (2hours). All candidate locations were ranked as Eq. 11.
R l was the ranking score of each candidate position P ijk .Variable r l is the distance from the candidate location to the original position P ijk , and the subscript ijk is the coordinates of the candidate locations. p (rl) is the visiting chance of P ijk (Eq. 12).
, 1, , 6 n l C n = … were the preference of neighbor at P ijk . The position P ijk had six immediate neighbors. The value of n l C was defined by Eq. 13.

( ) if the n-th neighbour is BMSC
The coefficients F E,0 and F E,max were the weights of Eq. 14 which described the degree of preference of an MIC attaching to an unprimed BMSC and the increase of such preference to a primed MBMSC. n xyz E was the stiffness at P ijk 's n-th neighbor with the exact coordinate (x,y,z). E 0 was the minimal sensed stiffness and E,max the maximal.
a neighbor as CD8 + T cell. If yes, the value of C CTL is small to decrease the chances of MICs migration to CD8 + ; otherwise, C CTL is large (Eq. 15).

at least a neighbor is CD8
Variable V l was defined as Eq. 16. ijk ijk ijk ijk 1 P has 5 6 neighbor cells 8 1 P has 3 4 neighbor cells 4 1 P has1 2 neighbor cells 1 P has no neighobr cells 16 Eq. 16 represents that myeloma cells try to avoid loneliness as well as crowdedness, though they are not as sensitive as solid tumor cells are. Ranks of all candidates were normalized so that the sum of all normalized ranks is 1 (Eq. 17).
To facilitate die-casting, all normalized ranks were incorporated to form a scale S, in which each candidate is corresponding to a range . For the CD8 + T cell, it will migrate towards MIC or MM rather than empty position and BMSC. The CD8-mediated lysis of myeloma cell might probably occur the positions where are not BMSC. In other word, myeloma cell adhere to BMSC will reduce the probability of apoptosis and escape from immune system [1]. The candidate rank was calculated with Eq.20.
In Eq.20, the values of p(r l ) and V l are also calculated by Eq.12 and Eq.16. The value of n l C for CD8 + T cell was defined by Eq.21.
As to the regulatory T cell (Treg), it will migrate towards CD8 + T cell and try to suppress the proliferation of CTL. The candidate rank was also calculated with Eq.20. In Eq.20, the values of p(r l ) and V l are also calculated by Eq.12 and Eq.16. The value of n l C for Treg was defined by Eq.22.

CD8 + mediated Lysis of myeloma cells
During the migration of CD8 + T cells, they moved towards to target cells (MICs or MMs) and tried to kill them. When a CD8 + T cell located at a position where there was a myeloma cell, the target cell was labeled with a status variable and will put into apoptosis list in the next time step.

Suppression effect of Tregs on CD8 +
When regulatory T cells were activated, they would sharply colon themselves. Tregs would move towards CD8 + T cells and tried to suppress the generation of CD8 + T cells by inducing cell cycle arrest or direct apoptosis.

Apoptosis
At each time step, if the apoptosis rate of the cell (MIC, MM, CD8 + , and Treg) was lesser than the threshold; all four types of cells in this system would start apoptosis. Cell took 10 time steps to finish apoptosis and was then absorbed. As we discussed in this study, cytotoxic drug BTZ increased the apoptosis rate of cancer cell as Eq. 23 (for MIC) and Eq.24 (for MM).
Where S a is the ODE-inferred survival rate based on the local stiffness and relative concentration of BTZ. S b is the CD8-induced survival ratio of myeloma cell [1]. As to MM cell, its apoptosis rate was defined as Eq. 26.
Where S c is the BTZ-induced survival ratio inferred from literature [4]. At last, the apoptosis rates of CD8+ and Treg were defined as Eq.28-29.

Tissue scale
In the tissue scale of this HABM model, the secretion of SDF-1 (from MIC) and TGF (From both MM and BMSC) and the diffusion of them in the 3D ECM defined the dynamic 3D distribution of SDF-1 and TGFβ, which were both described as Eq.30.