**2. Methods**

The present computational model was developed to analyze the interactions between MMCs and their ECM. It couples particle dynamics with fluid dynamics, considering cell–cell interactions by analyzing the cells using the discrete element method (DPM-DEM), implemented in Fluent (Ansys [31]) (ANS). The results were compared with those of a model developed previously in-house, using the commercial software Abaqus (ABQ) [32], which was adapted and calibrated with appropriate parameters to simulate MMCs [33,34]. Cell processes in both models were modulated by the specific cell conditions and implemented via UDFs. Both models were used to analyze the progression and development of MM tumor aggregation.

### *2.1. Multiphase Coupled Fluid-Particle Model (ANS)*

In this article, we present a hybrid computational model, which combines the formulation of two classical methods: continuous and discrete methods. Indeed, due to the

mathematical incompatibility of both models, each one calculates the corresponding phase in parallel, and the bidirectional interactive coupling of the methods is achieved through the intercommunication of the properties and conditions of each phase. Fluid dynamics are included in the cell through the inclusion of velocity, density, and viscosity values of the fluid in the motion equation of the cell. Besides, the perturbation of the continuous phase, due to the motion of the cells, is included as an external contribution to the momentum balance equation. Therefore, although both models are calculated in parallel, consistency and stability are guaranteed through common variables such as cell and fluid velocities.

In the discrete phase, our objective was to determine the cell velocity, **v***c*, at every time step, *t*. In such a case, cell motion was transferred as a moment applied to the particle corresponding to the drag forces applied during the cell–ECM interaction. Other contributions, such as the gravity or collision forces, were also taken into account, with the trajectory of the discrete phase being obtained using the force balance acting on the cell as [31]:

$$m\frac{d\mathbf{v}\_c}{dt} = \mathbf{F}\_{drag} + \mathbf{F}\_r + \mathbf{F}\_{\text{grav}} + \mathbf{F}\_{ij\prime} \tag{1}$$

where *m* is the cell mass, and **<sup>F</sup>***drag*, **F***<sup>r</sup>*, and **<sup>F</sup>***grav* are the corresponding contributions of the drag force, domain motion, and gravity, respectively. As the considered experiments have been developed in static conditions, in this paper, there is no contribution of **F***r* component. Besides, the contribution of the gravity forces is defined by considering differences in the cell and media mass density as:

$$\mathbf{F}\_{\mathcal{S}^{\text{grav}}} = m \frac{\mathbf{g} (\rho\_{\mathcal{E}} - \rho\_f)}{\rho\_{\mathcal{E}}},\tag{2}$$

where **g** is the gravity acceleration, *ρc* is the cell mass density, and *ρf* is the fluid mass density.

Other forces acting on the cell, such as cell–cell, and cell–wall collision forces, are referred to as **<sup>F</sup>***ij*.

The drag force on cells, as calculated using the cell constant spherical shape, is proportional to the difference between the cell and fluid velocity, as follows:

$$\mathbf{F}\_{\text{drag}} = \frac{1}{2} \mathbb{C}\_{D} \rho\_f A\_{\text{c}} \parallel \mathbf{v}\_f - \mathbf{v}\_{\text{c}} \parallel (\mathbf{v}\_f - \mathbf{v}\_{\text{c}})\_{\text{\textquotedblleft}} \tag{3}$$

where *CD* is the drag coefficient, *Ac* is the cell cross-section, and **<sup>v</sup>***f* is the fluid velocity.

In the event of overlap, cell–cell and cell–wall collision forces were taken into account using the discrete element method (DEM) (Figure 1). These collision forces were described by defining the stiffness of the contributing parts (cell–cell or cell–wall) and the overlap distance, *<sup>δ</sup>ij*, which was calculated as follows [31]:

$$\delta\_{\vec{i}\vec{j}} = \left| \left| \mathbf{x}\_{\vec{i}} - \mathbf{x}\_{\vec{j}} \right\|\right| - (r\_{\vec{i}} + r\_{\vec{j}})\_{\prime} \tag{4}$$

where **x***i* and **x***j* are the positions of the participating cells, while *ri* and *rj* are their radii. The cell–cell contact direction, **<sup>e</sup>***ij*, can be defined as:

$$\mathbf{e}\_{ij} = \frac{(\mathbf{x}\_i - \mathbf{x}\_j)}{\|\mathbf{x}\_i - \mathbf{x}\_j\|}. \tag{5}$$

**Figure 1.** Forces acting on the cell. The cell–cell overlap distance and cell stiffness were considered when computing contact forces. Drag forces were determined using the momentum exchange between the cell and the fluidic ECM.

The Hertzian–Dashpot Collision Law was then used to establish the collision normal force, **F***<sup>n</sup>*,:

$$\mathbf{F}\_{\rm n} = (K\_{\rm H} \delta\_{ij}^{3/2} + \gamma(\mathbf{v}\_{ij}\mathbf{e}\_{ij}))\mathbf{e}\_{ij\prime} \tag{6}$$

where *γ* is the damping coefficient, and **<sup>v</sup>***ij* represents the relative velocity of the collision, defined as **<sup>v</sup>***ij* = **v***j* − **<sup>v</sup>***i*. Finally, *KH* is the stiffness coefficient of the collision, which was calculated using the stiffness of the participating cells, as:

$$K\_{H} = \frac{4}{3} \frac{E\_i E\_j}{E\_j (1 - \nu\_i^2) + E\_i (1 - \nu\_j^2)} \sqrt{\frac{r\_i r\_j}{r\_i + r\_j}},\tag{7}$$

where *Ei*, *Ej* and *νi*, *νj* are the stiffness and the Poisson coefficient of the cells involved, respectively.

The tangential force, **Ft**, in the contact area described by the friction coefficient, *μ*, is defined as:

$$\mathbf{F}\_t = \mathbf{F}\_n \boldsymbol{\mu}.\tag{8}$$

The contribution of the collision force, **<sup>F</sup>***ij*, to cell motion (Equation (1)) was then calculated for each cell using tangential and normal forces as follows:

$$\mathbf{F}\_{ij} = \sum\_{i=j}^{n} \left( \mathbf{F}\_n^{ij} + \mathbf{F}\_t^{ij} \right). \tag{9}$$

For its part, in the continuous phase, the objective is to determine the fluid velocity, **<sup>v</sup>***f* , taking into account its active interaction with the cells. In such a case, the momentum exchange between the continuous and discrete phases, **F***ex*, is the change in the particle's momentum as it goes through each control volume at each flow time, Δ*t*. During the momentum balance, this momentum exchange is regarded as a momentum source in the continuous phase and is defined as [31]:

$$\mathbf{F}\_{\rm ex} = \sum \left[ \frac{18\mu C\_D Re}{\rho\_c (2r\_i)^2 24} (\mathbf{v}\_c - \mathbf{v}\_f) \right] \dot{m}\_c \Delta t\_\prime \tag{10}$$

where *μ* is the fluid viscosity, *Re* is the relative Reynolds number, *ρc* is the cell density, and *m*˙ *c* is the mass flow rate of the discrete phase.

Thus, the general equation for the momentum conservation of the continuous phase results in [31]:

$$\frac{\partial}{\partial t}(\rho\_f \mathbf{v}\_f) + \nabla(\rho\_f \mathbf{v}\_f \mathbf{v}\_f) = -\nabla p + \nabla \tau + \rho\_f \mathbf{g} + \mathbf{F}\_{\text{ext}} \tag{11}$$

where ∇ is the divergence operator, *p* is the static pressure, **g** is the gravitational acceleration, and *τ* is the contribution of the stress tensor given by:

$$\tau = \mu \left[ (\nabla \mathbf{v}\_f + \nabla \mathbf{v}\_f^T) - \frac{2}{3} \nabla \mathbf{v}\_f I \right],\tag{12}$$

where *I* is the unit tensor.

Cells evaluate and respond to the specific mechanical conditions to which they are subjected during cellular processes, such as proliferation and differentiation [35,36]. Thus, we defined the cell-cycle progression as the maturation index (MI), which is a time-dependent process, as:

$$\text{MI} = \begin{cases} \frac{t\_c}{t\_{\text{mat}}} & t\_c < t\_{\text{mat}} \\\\ 1 & t\_c \ge t\_{\text{mat}} \end{cases} \tag{13}$$

where *tc* is the cell-cycle progression time, and *tmat* is the time needed for the cell to undergo cellular processes, such as proliferation and differentiation [30]. In addition, the maturation time is dependent on the mechanical conditions of the cell. In this model, due to the consideration of the fluidic environment, the mechanical conditions perceived by the cell were defined in terms of the number of cell–cell (*Nc*) and cell–wall (*Nw*) contacts at each time step increment as:

$$t\_{mat} = \frac{t\_n}{1 + \gamma\_w N\_w + \gamma\_c N\_c} \,\tag{14}$$

where *γw* and *γc* are the proportional time factors for the cell–wall and cell–cell interactions, respectively, and *tn* is the natural maturation time, which is the time needed by the cell to complete the cell cycle without stimulus. As a cell progresses through the cell cycle, it becomes ready to divide into two new cells, thereby implying that the cell expands in size before proliferating. Thus, we assumed that the cell radius increases proportionally to the state of maturation as follows:

$$r\_i = r\_0(1 + \frac{\mathbf{M}\mathbf{l}}{2}),\tag{15}$$

where *ri* is the cell radius for each *i*th cell at each time step, and *r*0 is the nominal cell radius obtained from the literature [37].

#### *2.2. Design of In Vitro Experiments and Data Collection*

Experimental data from MMC lines were used to validate the model. The MMC line RPMI8226 was cultured in vitro in multi-well plates. RPMI8226 was purchased from the American Type Culture Collection (ATCC, Rockville, MD, USA) and cultured in RPMI1640 medium supplemented with 10% fetal bovine serum (FBS), 2 mM L-glutamine and 100 μg/mL penicillin, and 100 μg/mL streptomycin. Cells were cultured at 37 ◦C under a 5% CO2 atmosphere in a Galaxy S incubator (Eppendorf New Brunswick, Hamburg, Germany). A total of 100,000 cells were seeded in 500 μL of medium in p24 multi-well plates. Every 24 h, 50% of the volume of culture medium was renewed. After 24, 72, and 96 h of culture, cell proliferation was determined using the colorimetric MTS assay, according to the manufacturer's instructions, pipetting the supernatant onto a 96-well plate and reading at 490 nm (Victor 3 microplate reader, Perkin Elmer, Waltham, MA, USA). The results obtained are presented in [38].

### *2.3. Model Assumptions*

In order to reduce the computational costs of the model, we made some assumptions. Thus, the cell geometry is considered to be spherical with no changes in its morphology [14,22]. We also considered a reduced, but significant, portion of the whole in vitro experiment to reduce the number of cells simulated. Cell properties were obtained, when possible, for the RPMI-8226 cell line to adequately compare results with the in vitro experiments (Table 1). In general, these properties were considered homogeneous for all participating

cells, although slight discrepancies are possible due to the non-homogeneous nature of MM tumors [9]. Finally, neither oxygen nor nutrient consumption were considered.


**Table 1.** Mechanical parameters considered in the model.
