A mechanobiological model for tumor spheroids evolution: application to glioblastoma

Spheroids are in vitro spherical structures of cell aggregates, eventually cultured within a hydrogel matrix, that are used, among other applications, as a technological platform to investigate tumor formation and evolution. Several interesting features can be replicated using this methodology, such as cell communication mechanisms, the effect of gradients of nutrients, or the creation of realistic 3D biological structures. In this paper, we propose a continuum mechanobiological model which accounts for the most relevant phenomena that take place in tumor spheroids evolution under in vitro suspension, namely, nutrients diffusion in the spheroid, kinetics of cellular growth and death, and mechanical interactions among the cells. The model is qualitatively validated, after calibration of the model parameters, versus in vitro experiments of spheroids of different glioblastoma cell lines. This preliminary validation allowed us to conclude that glioblastoma tumor spheroids evolution is mainly driven by mechanical interactions of the cell aggregate and the dynamical evolution of the cell population. In particular, it is concluded that our model is able to explain quite different setups, such as spheroids growth (up to six times the initial configuration for U-87 MG cell line) or shrinking (almost half of the initial configuration for U-251 MG cell line); as the result of the mechanical interplay of cells driven by cellular evolution. Indeed, the main contribution of this work is to link the spheroid evolution with the mechanical activity of cells, coupled with nutrient consumption and the subsequent cell dynamics. All this information can be used to further investigate mechanistic effects in the evolution of tumors and their role in cancer disease. Author summary Spheroids structures of cell aggregates are an available experimental platform to analyze the evolution and drug response of solid tumors. In particular, the dynamics of different glioblastoma cell lines have been studied in this work using spheroids. Interestingly, very different behaviors were observed, from a half of the initial configuration shrinking for U-251 MG cell line to six times the initial configuration growth for U-87 MG cell line. These results were replicated by means of a coupled mathematical model which accounts for nutrients diffusion in the spheroid, kinetics of cellular growth and death, and mechanical interactions among the cells. Tumor growth or shrinkage can be explained from a continuum mechanics view driven by cell activity and nutrients availability. This modeling put the focus on mechanistic effects and is aligned with novel experimental techniques to quantify the mechanical microenvironment in tumors. These techniques may be combined with the approach presented in this work to further investigate the role of mechanics in cancer disease.

Spheroids have been used to study the mechanisms of cell differentiation, cell contact 18 and the different cell phenotypes that give rise to different regions within the 19 spheroid [11]. Furthermore, Friedrich et al. [12] developed platforms to study anticancer 20 treatments in spheroids. Nath et al. [13] presented several spheroid generation 21 techniques and assays for their characterization. On the other hand, Mark et al. [14] 22 studied the relationship between the cellular mechanics of the spheroid and the 23 deformation of the extracellular matrix; and Ayuso et al. [15] established that spheroids 24 represent a good tool in the study of migration by chemotaxis, the main mechanism of 25 cell migration. In another study, Guillaume et al. [16] suggested that the stresses that 26 originate in a spheroid during its growth affect its physical properties and therefore the 27 response to treatment. 28 From a different perspective, Wu et al. [17] developed microfluidic devices for the 29 study of chemotaxis, through the generation of chemical gradients in both simple and 30 complex environments, and studied chemotaxis as a function of the cell type. The 31 devices designed by Cheng et al. [18] allow a quick reconfiguration for use in different 32 applications, including the analysis of tumor first development and evolution. Ayuso et 33 al. [19] performed a protocol for the confinement of the gel in microfluidic chips and, mathematical model of cell dynamics, growth and death, in a microfluidic device with 45 application to glioblastoma evolution. Laird [24] performed a mathematical model of 46 tumor proliferation with different growth stages. Furthermore, Pettet et al. [25] implemented a mathematical model of cell migration in spheroids, where chemotaxis 48 plays an important role. Sendoel et al. [25] analyzed how hypoxia-inducing factors are 49 related to increased treatment failure and mortality, while Curtis et al. [26] studied the 50 relationship between cell mechanics and proliferation, concluding that compression 51 increases proliferation. Additionally, Bull et al. [27] recently proposed a discrete 52 (agent-based) approach to model the dynamics of cells in spheroids, accounting for the 53 diffusion of nutrients and oxygen from the medium to the interior of the spheroid, which 54 determine the cellular actions of proliferation, differentiation and cell death. In the 55 same context, Amereh et al. [28] present a mathematical model of tumor formation in 56 spheroids based on reaction-diffusion equations. 57 Our application is devoted to analyze glioblastoma spheroids. Glioblastoma is the 58 most malignant brain tumor, with a life expectancy after the standard treatment 59 (surgery plus radio and/or chemotherapy) of about 14 months from diagnosis [29]. It 60 originates in the glial cells and its initial evolution is controled by thrombosis, hypoxia, 61 cell migration events that occur successively [30] and is characterized by microvascular 62 hyperplasia and necrosis. The area where the latter takes place (necrotic core), is 63 usually surrounded by a high cell density, known as 'pseudopalisade'. The hypoxic tissue 64 (more resistant than a healthy one [31]) increases tumor glycolysis to provide energy to 65 the cell, promotes angiogenesis, that is the creation of blood vessels and nourishes the 66 tumor invasion and metastasis by the action of hypoxia-inducible factors [32][33][34]. 67 In this article, we develop a phenomenological mechanobiological model to 68 investigate the in vitro evolution of suspended glioblastoma spheroids. The model takes 69 into consideration the main physics observed during the process, such as nutrients   Medical Technology) with 5% CO2. Spheroids were generated by hanging drop method. 98 Cells were trypsinized and resuspended in growth medium supplemented with 20% 99 methylcellulose solution to reach the concentration of 40 000 cells/mL. In order to form 100 1000-cell spheroids, 25 uL drops were placed on the lid of a Petri dish. The bottom part 101 of the dish was filled with distilled water to prevent evaporation. Plates were placed 102 within the incubator and left for 48h to assure spheroid formation. Afterwards, 103 spheroids were transfered to growth medium in suspension 96 well plates (Sarstedt 104 83.3925.500) treated with anti-adherence solution (Stemcell 07010). Their behavior was 105 followed by phase contrast microscopy (Nikon TiE) during two weeks. Images were 106 analyzed with a Fiji plugin SpheroidJ, which allows automated quantification of the 107 spheroid area [35]. Growth curves were obtained by normalizing those values to the area 108 that spheroids occupied at the beginning of the experiment (day 0).    The diffusion of nutrients (considering nutrients as oxygen and glucose concentration) is 120 assumed, as a first approach, to follow a Fickean diffusion model, analogously to other 121 works in a similar context [26,[36][37][38].
The right-hand side of eq (1a) also includes a second (source) term which accounts 123 for the consumption rate of nutrients by cells. k c is the diffusion coefficient of the 124 spheroid aggregation, r is the rate of nutrients consumption by cells, and c 0 the initial 125 concentration of nutrients within the spheroid. On the other hand, the function f represents the consumption of nutrients available in the suspension experiment, defined 127 as: with α a factor (shape) coefficient which accounts for the amount and availability of 129 nutrients and configuration of the suspension experiment for diffusion of nutrients. This 130 coefficient is introduced as a correction factor to the full and plenty availability of 131 nutrients hypothesis in the experiment. c Γ represents the nutrients concentration in the 132 suspension at the initial time of the analysis. The solution of eq (2) yields: Cell dynamics 134 The dynamics of cellular growth/death is represented, similar to other works [24,39,40], 135 by means of the following first order equation: In eq (4), and as a simplification of other models [23,32,41,42], the kinetics of 137 cellular growth/death is characterized by the function k n (c), which is assumed to be 138 dependent of the nutrient concentration [23]. k n (c) is the function shown in Fig 1,   139 where two regions can be identified: (i) a cell death (or dormancy) region which appears 140 for low nutrient concentration, and (ii) a cell growth zone for high levels of nutrients. Therefore, eq (4) can be split into two different conditions for death and growth 146 dynamics as follows: where k d and k g represent death and growth constants, and the superscripts mean d 148 (death) and g (growth), respectively. Finally, the total cell concentration can then be 149 computed as, Cell mechanics

151
For the sake of simplicity of the numerical framework, we will use an updated cytoskeletal behaviour [43][44][45]. Hence, The cell contractility ε cell is assumed to be isotropically dependent on the cell 162 concentration n as follows: where k ε is the contractility constant, which depends on the available room for cells to 164 exert contractile forces.
The remaining equilibrium, compatibility and constitutive equations (linearized for 178 time t around the current configuration Ω(x, t)) read as, with C being the fourth order elasticity tensor which characterizes the mechanical 180 behaviour of the spheroid. σ 0 is a certain initial pre-stress in the spheroid, andt n and 181 u the prescribed tractions and displacements on the spheroid boundary. These Space and time are now non-dimensionalized as follows: with D 0 the diameter of the spheroid, and T 0 a characteristic time of analysis. On the 188 other hand, dimensionless nutrients concentration is defined with reference to the initial 189 concentration of nutrients in the suspension experiment: while the cell concentration is non-dimensionalized versus the initial cell concentration 191 of the spheroid: Finally, the stresses are referred to the elastic modulus (E) of the spheroid: Introducing the dimensionless quantities in eqs (12)- (15) in the nutrient diffusion eq (1) 195 yields, where the following dimensionless parameters in eq (17) are defined: Cell dynamics 198 The dimensionless quantities defined in eqs (12)- (15) are used in eqs (5) and (6), getting 199 appearing the following dimensionless parameters in eqs (19) and (20): Now, taking the time derivative of eq (21) yields, Considering eqs (19) and (20), eq (23) may be rewritten as: with, 203    δ 1 = 1 ifc ≤c n and δ 1 = 0 elsewhere δ 2 = 1 ifc n <c ≤c sat and δ 2 = 0 elsewhere δ 3 = 1 ifc >c sat and δ 3 = 0 elsewhere (25) Eq (24) will be used to compute the cell evolution in the spheroid. Dead and alive 204 cells can be computed afterwards from eqs (19) and (20).

216
Therefore, the governing equations of the model are solved at each j-step at time t j for 217 the updated configurationx j as described in Box 1.

218
First, after initialization of the field variables and initial domain, a backward-Euler 219 scheme is followed to discretize in time the nutrient diffusion and the cell dynamics

227
Eqs (32) and (33) are spatially discretized following a finite element (FE) numerical 228 framework. Thus, these eqs are firstly written in their weak form (the reader is referred 229 to [48][49][50] for the basics of FE analysis). Then, the field variables in eqs (32) and (32) 230 are approximated as follows: where Nc, Nñ and Nũ are shape function (interpolating) matrices for the continuum 1. Set j-step j = 0.

Parametric analysis 246
In order to investigate the influence of the different model parameters on the evolution 247 of the field variables and results, the dimensionless model shown in the previous sections 248 was repeatedly run for a set of parameters. Therefore, we vary up to 4 parameters from 249 October 7, 2021 12/28 their reference value as shown in Table 3. The selected parameters have a clear 250 physiological relevance such as nutrient consumption of cells (r), kinetics of cell death 251 and growth (k d andk g , respectively); and cell compression and expansion constants 252 (k ed andk eg , respectively). From the reference value of the parameters, the variation 253 was set for all cases from half to double the reference value (see Table 3). Sincek d −k g 254 andk ed −k eg are varied at the same ratio, they are referred ask andk e , respectively. 255 Table 3. Dimensionless values of the model parameters for the parametric analysis.

271
In this section, we validate the proposed mechanobiological model with the  Table 4), is also included in Fig 10 for shown for a spheroid in Fig 12. These results are compared for differed spheroids versus 289 the model outcome in Fig 13. In this case, the model was calibrated for a different set 290 of parameters, as seen in Table 4, since they are referred to different cells and behavior. 291 Using this calibration, we can see in Fig      Diffusion of nutrients and proliferation/death of cells in spheroids, are well-known and 300 important phenomena referred in different works [24,26,[36][37][38][39][40]. However, our main 301 contribution was to link the spheroid evolution with the mechanical activity of cells, 302 coupled with nutrient consumption and the subsequent cell dynamics. As a result, our 303 model predicted a higher proliferation rate for abundant levels of nutrients, and cell 304 death/dormancy for low levels of nutrients. This behavior is aligned with other models 305 and observed experimental behaviors [23]. Moreover, a material point exerts a high  3). This 310 hypothesis is similar to the behavior observed in agent-based models [22,27]: in these 311 models cells are considered as (discrete) contractile elements which evolves until a 312 (contracted or expanded) equilibrium state that accommodates the new daughter cells. 313 This behavior is also aligned with experimental evidences [46,47].

314
The developed mechanobiological model contains a number of phenomenological 315 parameters with physiological meaning. In fact, these parameters may be measured (or 316 calibrated) using standard experiments of proliferation, diffusion [20], or traction force 317 microscopy [57,58] to account for the contractile behavior of glioblastoma cells. The 318 effect of model parameters on the results was investigated in a parametric analysis, 319 using the same rate of variation (from half to double) for all the parameters. It can be 320 seen that contractility constantsk e andñ * have a minor effect on the living cells consumption rates and a high rate of death cells, as expected. In these situations, dead 329 cells rise predominantly in the interior of the spheroid at the last steps of the analysis 330 (see Fig 5). As seen in Fig 6, the most influential parameters on nutrient diffusion are 331 the nutrient consumption and kinetics of cell death/growth constants. Specifically, there 332 is more availability of nutrients for low values of nutrient consumption and a low  Fig 4) along time.

340
The trend of the evolution of the spheroid area is qualitatively similar for all the 341 analyzed cases (see Fig 8): there is an initial decay in the first steps (sooner or later 342 depending on the parameters), followed by an increasing growth trend until the end of 343 the analysis. The highest spheroid growth at the end of the analysis appears for the 344 higher kinetic constant (up to 150%), followed by the cases of higher 345 contractile/expansion constantk e and lowerñ * (volume availability or spheroid 346 compaction at the initial step). The highest spheroid contraction in the first steps of the 347 analysis is produced for higher contractile/expansion constantk e and highñ * . It is also 348 interesting to remark that the minimum of growth decay (spheroid compression) is  Table 4. We maintained the reference values of k c , c n and c sat , while we 357 calibrated the consumption of nutrients, proliferative capacity or contractile behaviour 358 for the different analyzed cell lines since they are related to specific characteristic and 359 physiology of the cell. For example, it was found during the calibration that U-87 cells 360 are more proliferative with a higher rate of nutrient consumption, according to the 361 fitted parameters in U-87. Interestingly, the fitted parameters for U-251 cell line 362 resulted into flat gradients of nutrients, and consequently a more homogeneous cell 363 distribution and absence of pressure along the spheroid (see Fig 14). Differences 364 between these cell lines were already studied and it was shown that U-87 cells are more 365 proliferative than U-251, and that U-251 cells can survive better in conditions of low 366 nutrients and oxygen concentrations [59,60]. In our experiments, U-251 spheroids were 367 small and did not grow, their diameter was around 200 µm, while the diameter of U-87 368 spheroids reached 800 µm in 14 days. It was observed in tissues that oxygen can reach 369 cells that are 100-200 µm away from blood vessel [38,61]. This was confirmed in 370 spheroids, as gradients of oxygen and nutrients, and consequent formation of necrotic 371 core depend on the spheroid size. Spheroids whose diameter was bigger than 500 µm 372 developed physicochemical gradients sufficient to induce necrotic core formation [13,62]. 373 Different cellular behavior, diffusion distances reported in literature and spheroid sizes 374 obtained in our experiments can explain our fitted parameters. As U-251 spheroids were 375 small, all the cells could get enough oxygen and nutrients, and, as these cells barely 376 proliferated in spheroid culture, gradients did not change during the experiment. On 377 the other hand, U-87 cells are highly proliferative, spheroids increased their size, so 378 gradients of oxygen and nutrients were formed. Cells in the interior of the spheroid did 379 not get components necessary for survival, so cell death was activated and necrotic core 380 was formed. The distribution of nutrients in the spheroid for this case can be seen in  Table 4.