# An in-silico study of cancer cell survival and spatial distribution within a 3D microenvironment

Jul 31, 2020

### Validation of SALSA

The validation of SALSA consisted in comparing, over time, a) the cell density measured in-vitro to the results obtained with the in-silico simulation, and b) the ability of the virtual cells of inducing a stiffening of the ECM comparable to what observed experimentally. This latter comparison is particularly important for the considered experimental system, where the two compared cell lines were shown to behave differently39,40. MDA-MB-231, in particular, demonstrated to be able to increase the compressive Young’s modulus of the scaffold by reorganizing its fibers and increasing their density.

As detailed in the methods section, the cell density was defined as the number of living cells with respect to the initial population and it was evaluated for the two different breast cancer cell lines over a period of 10 days. Every condition was simulated 50 times and the average results are shown in Fig. 1a,b. as square markers, while the data obtained in-vitro are shown as circles. A good agreement between the two datasets is observed throughout the whole considered time frame. Indeed, the distribution of the mean average percentage error (MAPE, Eq. 1, Fig. 1c,d.) was shown to have a median within the experimental variability (22% and 25% for the MDA-MB-231 and MCF7 cells respectively) despite an initial slightly larger distance among in-silico and in-vitro data, likely caused by variability in the efficiency of the in-vitro cell seeding procedure, that was not represented in the model.

$$MAPE = frac{{left| {silico – vitro} right|}}{vitro}$$

(1)

The differential behavior of the two cell lines is evident. Epithelial-like MCF7 cells, Fig. 1a, show an initial short proliferative period, followed by a rapid decrease in the number of cells that leads, at the end of the simulation, to a cellular density that is approximately half the initial one. On the contrary, mesenchymal-like MDA-MB-231 cells, in Fig. 1b, proliferate at a nearly constant rate for the first 3 days. Afterwards, the cell density decreases and reaches at day 10 a value comparable to the initial one. This limited proliferation has been observed in other 3D cell culture systems41, where it was associated with the dynamics of diffusion of nutrients and oxygen through the matrix.

The differential behavior of the two cell lines is also reflected in their ability of modifying their environment. The effect breast cancer cells on the stiffness of collagen scaffolds was previously studied39. Briefly we observed that MDA-MB-231 cells can significantly increase its value from 46.9 ± 5.3 to 57.9 ± 7.0 kPa, 10 days after cells were put in culture within the scaffolds (Kruskal–Wallis test, p = 10−6). Conversely, MCF-7 cells did not induce significant modifications of the scaffold stiffness. This was shown to correlate with the expression of Lysil-Oxydase (LOX), an enzyme responsible for collagen cross-linking that is 1,000-fold more expressed in MDA-MB-231 compared with MCF7 cells. Including in SALSA this differential LOX expression leads to the results reported in Fig. 2.

Simulated data are again able to consistently replicate these experimental results (Kruskal–Wallis test p = 10−20), the only difference being the variability obtained at day 10 for the scaffolds where MDA-MB-231 cells were cultured (Fig. 2b). This higher 95% confidence interval is caused by an inherent difference between the two datasets. Indeed in-vitro variability was obtained as the 95% confidence interval of the average stiffness of different scaffolds, while in-silico it was computed as the dispersion of the Young’s moduli distribution, obtained combining all the available data. Using the former approach to elaborate the in-silico data leads to a 15-fold reduction in variability that however does not account for inter-scaffold differences. As such the latter approach, can be considered to be more accurate, especially for MDA-MB-231 cells that exert a significant activity on their environment.

This analysis was completed by the study of how the main parameters of the model influence the results shown in Figs. 1 and 2. This was obtained through a global sensitivity analysis, fully detailed in the methods section, that allowed to compute the first order and total Sobol indices for each parameter (a, b, c, d, e, UGlu, UO2, UYM, s) listed in
Table 1 (see also Fig. 3 and Supplementary Figs. 6, 7 and 8) and the associated temporal evolution of the standard deviation (Supplementary Fig. 9). First order and total Sobol indices were shown to be approximately equal showing that the model’s parameters have minimal interactions. The most relevant change pertains the behavioural parameter c, that drives the transition between proliferation and quiescency and its role in the determination of the Young’s Modulus. This is likely determined by the central role of proliferant cells in the modification of the surrogate matrix properties.

Limiting the analysis to the first order indices, cell density was shown to be mainly determined by a, c, d, and s, albeit with different intensities, while the matrix stiffness was shown to be more influenced by a, b, d, e, UYM and s. As described in the supplementary material, some of these parameters were determined from experimental data (a, e, s, UY M) while others were computationally optimized (b, c, d). This compromise was necessary due to the lack of experimental evidences regarding certain aspects of the model, like the transitions between different cell states. We assume that overfitting risk is minimized having extensively used literature data for parameter estimation. Moreover only part of the experimental data was used for this purpose. The resulting model was indeed able to predict the value of the rest of the experimental results.

Stochasticity, on the other hand, was shown (Supplementary Fig. 9) to progressively increase as the simulation proceeds. This is partly due to the uniform starting condition imposed in our simulations, and is manly determined by the probabilistic nature of the rules and the procedures used to update cell status. Environmental parameters were shown to dominate behavioural ones when considering the Young’s modulus as output, while both classes of parameters equivalently determined cell density.

In the following sections, we will present the use of this computational tool to study important aspects of 3D cultures that are difficult to assess experimentally, such as the relationship between cell localization and viability, the local matrix stiffness and the distribution of oxygen and glucose within the scaffold. Finally an example of how SALSA could drive the in-vitro analysis will be shown. The simulation of three alternative initial cell densities will be presented and the experimental condition capable of granting sustained growth of the virtual population will be identified.

### Study of local variables using SALSA

A fundamental characteristic of 3D cultures, that makes them more physiologically representative than their 2D counterpart, is that distinct locations within the scaffold display differential microenvironments due to the presence of a nutrients gradient from the external layer to the inner scaffold core. Measuring these differences in-vitro, however, is particularly challenging due to the lack of high resolution quantitative techniques.

SALSA can be used to address these limitations as it tracks the location of each cell and the distributions of oxygen, glucose and Young’s modulus with spatial and temporal resolutions of 1 mm and 1 h respectively. This information can be used to complement the experimental analysis and retrieve valuable information difficult to obtain otherwise. This concept is exemplified in Fig. 4, where the results of these simulations are represented highlighting the effect of the distance from the center of the scaffold on each variable. For the simulated results the Manhattan distance was substituted to the euclidean one, since the cubic lattice used for the simulation is not radially symmetric.

When considering cell density, the color scale represents the fraction of living cells normalized with respect to the cardinality of the initial population. Although relatively uniform at the beginning of the simulation, this value rapidly decreases in the scaffold core (∆ center < 6 mm) while it stabilizes (MCF7) or increases (MDA-MB-231) in the most peripherical regions. Albeit coherent with the global results presented in Fig. 1 this prediction points at the existence of two radically different micro-environments within the 3D structure, one compatible with cell growth and survival and another associated with population collapse. Imaging of a scaffold section with a confocal microscope confirmed this result in-vitro40. Indeed the more aggressive MDA-MB-231 cells were capable of migrating toward more favourable environments for survival, while MCF7 cells maintained their original approximately uniform distribution and were unable to proliferate effectively.

The high mortality rate in the scaffold core seems to be connected with glucose availability, as the initial (24–48 h) decrease in the average concentration of this nutrient in the core of the scaffold is synchronized with cell density reduction (Fig. 4). Notably, despite a comparable reduction in glucose level, MDA-MB-231 cells in the external layers of the scaffold are predicted to be in the proliferative status. Their proximity to nutrient and oxygen rich medium is a plausible explanation for this predicted behaviour.

Simulated oxygen levels don’t decrease prominently at the beginning of the simulation. Their lowest concentrations were registered for the MDA-MB-231 cell line toward the end of the experiment. This might be connected with both the higher diffusivity (more than fourfold) and the lower cell uptake ( 104) of this molecule, when compared to glucose. The consequential reduced abruptness in the drop of oxygen concentration might reduce the impact of this variable on cell behaviour.

Finally the average simulated Young’s modulus displays two different behaviours for MCF7 and MDA-MB-231 cells. In the former there is no clear difference with respect to the initial condition, while the latter exhibit a dependence on the distance from the center of the scaffold similar to cell density distribution. Indeed a 10% difference in stiffness between the more rigid external shell and the softer core is predicted for MDA-MB-231 at the end of the simulation. This result suggests that these cells might generate an anisotropic material with more complex mechanical properties and behaviour than the initial structure.

The increased resolution granted by SALSA simulations was here shown to be potentially instrumental for the study of the complex feedback loops that govern the interaction between the cells and their environment, as the simplified dynamics and the immediate access to all the variables of interest could aid the optimization of the experimental setting and the interpretation of population-level results.

### Study of impact of differential initial population cardinality

The analysis of the SALSA simulations described in the previous section highlights how an initial density of 5 M cells might not be ideal, as it was associated with non-viable conditions in the innermost part of the virtual scaffold. To test if a different initial population density could increase the uniformity of cell distribution and stiffness we simulated the same experiment considering starting populations of either 625 K, 1.25 M or 2.5 M cells. The results of this in-silico analysis are reported in Figs. 5 and 6, for MCF7 and MDA-MB-231 cell lines respectively.

In these simulations, a lower initial cell density importantly reduced the glucose depletion in the core region at the beginning of the experiment, especially when considering starting populations of 625 K and 1.25 M cells. This was associated with a more pronounced and uniform increase in cell density. A similar distribution was also observed in the simulated scaffold stiffness, where starting populations with fewer than 1.25 M cells were shown to be associated with a low intrascaffold variability (5% at most).

These predicted features are valuable drivers for the design of in-vitro 3D culture protocols, where more homogeneous cell populations are expected to reduce the inter-experimental variability and to increase the significance of outcome. To confirm these results, we used confocal microscopy for the imaging of the cell spatial distribution within the whole scaffold.

The initial seeding of 5 M MDA-MB-231 cells determined an edge region significantly more populated than the core, after 7 days of culture40, while an initial concentration of 1 M cells was associated with a more homogeneous cell distribution within the scaffold without statistically-significant differences observed over time between core and edge regions (Supplementary Fig. 10).

Additionally, a further analysis of our simulated data revealed that smaller initial cell populations, although being unable to generate the cell densities obtained seeding 5 M cells, achieved final cardinalities that were almost independent from the initial condition (Fig. 7a,b), that is coherent with the experimental results shown in Bitar et al. (2008)42. Room availability within the scaffold appears not to be the reason for this result, that is likely to be rather determined by the establishment of an equilibrium between cellular growth and nutrients availability.

Normalizing with respect to the initial population (Fig. 7c,d) highlights how lower population densities are associated with a more pronounced growth and generate populations that can be sustained by the system for the entire experiment.

The overall results of this section, while mainly qualitative and limited to cell density, show that an initial population of 625 K cells could be the best tradeoff between the amount of living cells and a homogeneous distribution within the scaffold at the end of the experiment. This is also supported by other experimental works43,44,45 that set the initial population cardinality between 300 K and 1 M cells.