Quantifying the role of mechanics in the free and encapsulated growth of cancer spheroids

Mechanics is of primordial importance to understand cancer; hence, experimental 14 and mathematical models providing quantitative information and play a pivotal role on the 15 development of new therapies. Within this context, encapsulated spheroids are emerging as 16 exceptional in vitro tools to investigate the impact of mechanical forces on tumor growth, since 17 from the deformation of the alginate capsule the internal pressure can be retrieved. We show 18 that bio-chemo-poro-mechanics is a suitable theoretical framework to understand, explain and 19 design these in vitro experiments. Such mechanistic models are based on a set of coupled partial 20 differential equations solved within an open-source framework (FEniCS). Through sensitivity 21 analyses, our mathematical model suggests that the main parameters determining the 22 encapsulated and free growth configurations are independent, this observation indicates that 23 radically different phenomena are at play. This demonstrates potentialities of reactive porous 24 media mechanics in oncophysics and can serve to supplement in vivo clinical data. 25


27
As a tumor grows, it deforms surrounding living tissues, which causes strains and associated stresses.   (2012). Alginate 44 capsules ensure favorable conditions for cellular growth, because its permeability permits the free 45 flow of nutrients and oxygen without requiring additional molecules that could be potentially toxic 46 for the MCTS. By surrounding a growing tumor by an alginate capsule and monitoring the defor-47 mation of the capsule, whose mechanical properties can be robustly estimated, it is possible to 48 estimate the forces exerted by the growing tumor on the surrounding capsule. More generally, 49 alginate capsules offer an approach to interrogate, in vitro, the interplay between tumor growth 50 and mechanical forces. As measuring strains in vivo on tissues is either complex or outright im-51 possible Nia et al. (2016), this reproducible and robust in vitro approach is helpful to investigate 52 the pressure exerted by growing tumors onto the surrounding tissues which constrain its growth. 53 These mechanical forces were suggested to contribute to tumor growth regulation Helmlinger et al. 54 (1997), Paszek et al. (2005), hence evaluating their magnitude, distribution and evolution over time 55 could be a key step in devising successful methods to prevent the spread of tumors.

56
The experiments allow us for monitoring the capsule strains and the MCTS radius during growth. 57 We analysed tumor cell density and evaluated whether the cells are in a proliferating or necrosis 58 state. These experimental results offer invaluable information for the evaluation of any mathemat-59 ical model for tumor growth.

60
Modeling approach 61 The understanding of the physics and mathematical modeling in oncology have made significant 62 progress due to our improved ability to measure physical quantities associated with the devel-63 opment and growth of cancer. These experimental investigations have shown that mechanical 64 stresses influence biology. In particular, it has been shown that such mechanical forces inhibit 65 tumor growth Helmlinger et al. (1997), influence cell differentiation (mechanotransduction) and 66 hence the acquisition of a cancerous phenotype Paszek et al. (2005). 67 Health research centers have been collaborating with engineers, mathematicians and physicists to 68 introduce mechano-biology within clinical practice. This is particularly true for biochemical and ge-69 netic approaches which have been validated on patient cohorts, e.g. for the prediction of surgical  Three physics-based modeling approaches are possible to model cancer: discrete, continuous, hy-73 brid. The reader is referred to more detailed descriptions in the work of Lowengrub et al. (2010). 74 The most advanced continuum models consider tumors as a composite system, fluids and solids,   This tumor growth modeling approach is a subset of a more general set of methods within the 85 theory of porous media and was also used to describe ulceration and necrosis of human plantar 86 tissues, for patient with diabetes Sciumè et al. (2014b).

87
In this article, we specialize this reference model to numerically reproduce the experiment of  ical capsule, a few hundreds of microns in diameter. These alginate capsules are built using a 101 microfluidic co-extrusion device: the outer sheath is made of a sodium alginate solution; an in-102 termediate compartment contains a calcium-free medium; and the inner core is composed of the 103 cell suspension (CT26 mouse colon carcinoma). Performing extrusion in the air, the liquid jet is 104 fragmented into droplets (due to Plateau-Rayleigh instability) which upon contact with a calcium 105 bath, readily crosslink as shells encapsulating cells (alginate undergoes gelation in the presence 106 of divalent ions). The capsule allows convection of interstitial fluid (IF) and diffusion of nutrients 107 species, growth factors and drugs through its surface; however, thanks to alginate pores' size (of 108 the order of 20 nm), cells cannot escape. The capsule therefore serves as micro-compartment for 109 the 3D cell culture. During growth, the tumor deforms its surroundings which reacts with a con-110 finement pressure. The mechanism is similar to what happens in the CCT experiment: during cell 111 proliferation, the fraction of capsule volume occupied by cells increases until the capsule is filled 112 (confluence); then the tumor spheroid starts to strongly interact with the capsule and deforms it.

113
After confluence, the alginate capsule, deformed by the MCTS, responds with a confinement pres-114 sure due to action-reaction principle. This confinement pressure and non-optimal oxygenation of 115 the MTCS core areas generate important measurable heterogeneities (necrosis, local increase in 116 cell density, etc.) along the spheroid radius.

117
CCT allows generating capsules with desired size and shell thickness. This can be achieved by  Before confluence, in both thick and thin capsules, the growth rate of the CT26 spheroid is almost the same as that of the free MCTS case, indicating that access to nutrients is not compromised by the presence of the alginate shell. After confluence, the behavior strongly deviates from that of the free MCTS case. In the previous article Alessandri et al. (2013), we observed qualitatively the same phenomenology in thick and thin capsules. However results are quantitatively very different due to the different overall stiffness of the alginate capsules. We monitored the evolution of MCTS radius, capsule strain and cells states (proliferating or not) over several days. Images were taken at regular intervals by time-lapse phase-contrast microscopy ( Figure 1D-G). MCTS volumes were computed from the measured radii, assuming a spherical geometry. To label the necrotic core of the tumor, we used sulforhodamine B (SRB) with two-photon microscopy ( Figure 1B,C). The proliferating cells were examined using Ki-67 staining (Figure 1D,E). We observed that the proliferating cells were uniformly distributed in the free growing MCTS Figure 1D) while after confluence, in the encapsulated MCTS, cell division principally occurred in the peripheral rim ( Figure 1E). Regarding the dilatation of the alginate capsule, we characterized it as an elastic deformation with negligible plasticity and no hysteresis. Young's modulus was measured by atomic force microscopy indentation and osmotic swelling, giving a range of 68 ± 21kPa Alessandri et al. (2013). Thanks to the identified Young's modulus, capsules can be used as a biophysical dynamometer as a relation can be constructed which relates the variation of the inner pressure with radial deformation, monitored using video-microscopy. The pressure exerted by the MCTS is calculated using the formalism of thick-walled internally pressurized spherical vessels. Assuming that the alginate gel is isotropic and incompressible the radial displacement of the inner wall, ( in ), reads where is the internal pressure, is the Young's modulus, and in and out are the inner and outer radii of the capsule, respectively. Alginate incompressibility also implies volume conservation of the shell. This gives the following constraint equation Using this equation, the two time variables in ( ) and out ( ) can be separated and pressure, ( ), written as a function of in ( ) only Experimentally, we thus need to measure only the initial outer and inner radii and track the evo- We have four primary variables: three are scalar, namely the pressure of the medium/interstitial 156 fluid, , the pressure difference between the cell phase and the medium/interstitial fluid , , the 157 mass fraction of oxygen, ̄ ; one is vectorial, the displacement of the solid scaffold, . We have 158 two internal variables: the porosity and the TC necrotic fraction . 159 We introduce two kinds of closure relationships for the system: mechanical and mechano-biological.

160
Details about derivation of the governing equations and this constitutive relationships are provided 161 in appendix Appendix A. The Multiphase System.

162
The parameters of these closure relationships are of critical importance and will be studied through  Computational framework 177 We implemented the above model in Python and C++ within the FEniCS framework Alnaes et al. (2015), with an incremental monolithic resolution of the mixed finite element (FE) formulation and assuming cylindrical symmetry. The simulations have been run with composite Taylor-Hood element 3 (ℝ 2 ), [ 2 (ℝ)] 3 (one vectorial and three scalar unknowns), a mesh cell size of ℎ = 5 m and an implicit Euler scheme with = 1200 s. An updated lagrangian approach is adopted to account for geometrical nonlinearity. All the details and analytical verification of the FE formulation can be found in Appendix B. Computational framework. For all the sensitivity or optimization topics, we measured the error by the root mean square error (RMSE) relatively to a reference, specified each time. The error on the numerical quantity num relative to a reference ex , evaluated at points is: Presentation of the two configurations 178 Figure 2 shows the two modeled configurations of MCTS. We assumed cylindrical symmetry in the 179 simulations: each mesh is a half of a sphere, because we also exploit symmetry with respect to a 180 diametrical plane.

181
For the three scalar variables, we prescribed Dirichlet boundary conditions along the outer radius  The following model outputs have been compared to the experimental results: volume, capsule 214 strain, TC density and necrotic fraction.

216
Based on a detailed sensitivity study, the optimized set of parameters was tested and crossvali-   Pa Initial guess 5944 Numerical simulation also provides a wide output of qualitative results (e.g. IF flux, necrotic fraction 219 and TC saturation, etc.) which are presented and interpreted at the end of the section.

231
For the encapsulated configuration, the most important parameter is the critical inhibitory pres-232 sure crit (43.5%). Only two parameters gather 54.9% of the cost function variation: crit and 1 233 (11.35%).

234
Optimization 235 We identified the three governing parameters , 0 , for the free MCTS configuration by the Nelder-

236
Mead simplex algorithm and fitted to the experimental data with a = 0.031. To be physi-237 cally relevant, the same parameters set should be shared by the two configurations, thus we have 238 injected these three parameters within the encapsulated configuration and we identified its two 239 governing parameters 1 , crit by the same algorithm. We fit the experimental data of the encap-240 sulation with a = 0.124. Figure 3A shows the two configurations fitted with the following 241 set of parameters: = 3.33 .10 −2 , 0 = 6.65 .10 −4 , = 890, 1 = 1432, crit = 5944 (see Table 1). This  Figure 3A shows the range of modeling possibilities of the optimized parameters, 252 respectively to the alginate stiffness range. Figure 3B shows the validation on 2 thick capsules 253 with an alginate stiffness at = 52.5kPa and = 70kPa respectively. The modeling range is in accordance with these experimental results. Figure 3C shows validation on one thin capsule and 255 the other is fitted with the lower stiffness value = 47kPa, as the linear elastic hypothesis applied 256 to alginate is limited assumption for this level of deformation (> 20%). However, we showed that the 257 model could adapt to different geometries without lost its relevance: Figure Figure 5A and B show the interplay between oxygen consumption and necrosis. Indeed as men-272 tioned in the experiments, 85 hours after confluence, the viable space remaining for TC is a 20 m 273 thick rim. This is explicit in Figure 4, upper right circle, NTC quarter. Figure 5C shows the IF pressure 274 evolution, after confluence, we see a sucking phenomenon due to IF absorption by growing TC. As 275 the cells activity decrease at the MCTS inner core, IF accumulates and its pressure becomes positive 276 during a certain time (see Figure 5C, 28h). After confluence the initial gradient of IF pressure (green 277 line in Figure 5C) inverses sense: cells in the proliferative peripheral areas move toward the core 278 so IF has to go in the opposite direction (due to mass conservation). After 2 days of quick growth, 279 experimentally and numerically, the MTCS reaches a state of linear and slow evolution and from 280 that point onwards, the IF flux will not qualitatively change. This overall dynamic is clearly visible 281 on Figure 5D and E, as the capsule displacement is driven by TC partial pressure. The comparison 282 of Figure 5F and B show a relation between the saturation of TC and their necrotic fraction. This 283 is a basic experimental fact that, when the cells bodies collapse in a necrotic core, the aggregate 284 density increase accordingly. 285 We reproduce another case presented in Alessandri et al. (2013) to study the quantitative level 286 of TC saturation and necrotic fraction in a 50 m radius capsule. Figure 6A and Figure 6B Figure 4. Qualitative results, experimental and 6 physical quantities from the modeling: oxygen, necrotic tumor cells, interstitial fluid pressure, radial displacement, partial tumor cells pressure and tumor cells saturation. Left, at confluence. Right, 85 hours after confluence. The mass of growing tumor cells is directly taken to interstitial fluid, this implies a sucking phenomenon at the inner capsule radius (down right, upper left quarter). According to the experiment, a 20 m thick viable rim is obtain 85 hours after confluence (upper right, right half), the tumor cell saturation depends both on partial pressure and necrotic fraction (down right, down right quarter).

307
The sensitivity analysis we performed shows that the parameters which govern the growth in free 308 and encapsulated conditions are almost independent, which suggests the physical phenomena 309 should be independent too. is not suitable when the strain regime is relatively large (Figure 3C). The rich output of the mathe-316 matical model provides information about phenomena not yet observable or quantifiable in vitro 317 (see Figure 4) and allows for quantification of living and necrotic cells along the capsule radius (Fig-318   ure 6). 319 This poromechanical approach offers a unified framework to model both the growing MTCS and non-proper constitutive relationship (or starting hypothesis) that is identified and adjusted or fully 328 modified. Following a deductive scientific methodology, we have systematically tested the math-329 ematical model not trying to validate it but conversely trying to highlight its limitations. This has 330 allowed to progressively improve it and finally to achieve a predictive (non-phenomenological) for-331 mulation.

332
In 2020s mathematical modeling in oncology begins to enter a stage of maturity; today mathemati-333 cal models of tumor growth tend to clinical applications and therefore must be really predictive and 334 funded on measurable or at least quantifiable parameters having, as much as possible, a sound 335 physical meaning. This motivated this paper which presents not only a mechanistic bio-chemo-336 poromechanical model but also a modus procedendi to achieve a certain predicative potential and, 337 with intercession of sensitivity analysis, to quantify relative relevance of mechanisms underlying 338 tumor growth phenomenology.     We give in this section all the details about governing equations and the constitutive relationships.

418
According to the different phases, solid scaffold , medium/interstitial fluid and tumor cells phase 419 , described in the main article, section Methods, which constitute the multiphase system, at each 420 point in the domain, the following constraint must be respected where is the volume fraction of phase . Defining the porosity as Equation 1 can also be expressed in terms of the saturation degree of the fluid phase, = ∕ 423 (with = , ) 424 + = 1. (3)

425
Mass conservation equations 426 We express the mass conservation equation for each phase. We use a material description for the 427 motion of the solid phase and a spatial description for the fluid phases, whose reference space 428 is that occupied by the solid scaffold. As the solid is deformable, this reference space is not fixed where is the material time derivative with respect to the solid phase, is the density of phase 436 , ̄ is the velocity vector of the solid phase, ∑ ∈ → is the total mass exchange (water, oxygen and 437 other nutrients) from the interstitial fluid to the tumor due to cell growth and metabolism, ̄ is 438 the relative velocity of cells, and ̄ is relative velocity of the interstitial fluid.

439
The tumor cell phase is a mixture of living (LTC) and necrotic tumor cells (NTC), with mass fraction 440 ̄ and ̄ , respectively. The following constraint applies Mass conservation equations for each fraction, assuming that there is no diffusion of both necrotic 442 and living cells, read where is the death rate of tumor cells. Note that only one of Eqs 8-9 is independent: actually, 445 one can be obtained subtracting the other from Eqn 5 and accounting for the constraint Eqn 7.

446
Oxygen is the only nutrient which we consider explicitly. Another mass balance equation is intro-447 duced which governs the advection-diffusion of oxygen, , within the interstitial fluid where ̄ is the diffusive velocity of oxygen in the interstitial fluid and → the oxygen consumed 449 by tumor cells due to their metabolism and proliferation rate.

450
Momentum conservation equations 451 We neglect here the effect of gravitational body forces as their contribution is negligible compared 452 to that of other forces. Furthermore, as we assume quasi-static processes and small difference in where ̄ is the stress tensor of phase , is the set phases connected to and → is the inter-457 action force between phase and the adjacent phases. Summing eqn 11 over all phases gives the 458 momentum equation of the whole multiphase system as where ̄ is the total Cauchy stress tensor acting on the multiphase system.

460
Assuming that for relatively slow flow, the stress tensor for a fluid phase, , can be properly approx- where is the averaged fluid pressure and the unit tensor, Eqn. 11 which apply for a generic 463 phase (solid or fluid) can be expressed in an alternative form for fluid phases as Sciumè et al. 464 (2013) where is a symmetric second order resistance tensor accounting for interaction between the 466 fluid phase and the solid phase, . Eqn. 14 can be rewritten as

474
Effective stress principle and closure relationships 475 We assume here that fluid phases are incompressible and the solid phase is almost incompressible.

476
However, the overall multiphase system is not incompressible, because of the presence of porosity 477 that evolves according with the scaffold deformation. As phases are incompressible, their densities 478 with = , , are constant. As the solid phase is quasi incompressible, the Biot's coefficient is 479 set to unity. With these premises, the total Cauchy stress tensor appearing in eqn 12 is related to 480 the Biot's effective stress as follows  The chosen closure relationship for the effective stress ̄ is linear elastic: with ( ) = 1 2 (∇ ̄ +(∇ ̄ ) ) and̄ ( , ) the fourth order elasticity tensor depending on the Young 485 modulus of the solid scaffold and its Poisson ratio.

486
The experimental measurement of cells density inside the capsule revealed a strong dependency 487 to necrotic fraction ̄ . Hence, the pressure-saturation closure relationship has been improved with pressure difference between tumor and interstitial fluid (i.e. = − ). The saturation is 491 directly linked to the partial pressure of the phase and a constant parameter , which represents 492 refinement of the porous/fibrous network. Its influence is offset by the necrotic fraction of tumor 493 cells, ̄ (see Fig.Figure 7), which allows us to model necrotic areas of very high cell density accord-494 ing with experimental evidence.

495
Fick's law is adopted to model diffusive flow of oxygen.

496
Tumor cells growth, metabolism and necrosis 497 ∑ ∈ → , the total mass exchange from interstitial fluid to the tumor cell phase is defined as with the tumor growth rate parameter, cell-line dependent.  and  are regularized step func-499 tions varying between 0 and 1, with two threshold parameters 1 , 2 :  =  ( , 1 , 2 ). When the 500 variable is greater than 2 ,  is equal to 1, it decreases progressively when the variable is between 501 1 and 2 and is zero when the variable is lower that 1 .  represent the growth dependency to 502 oxygen: env , the optimal oxygen mass fraction, is set to 4, 2.10 −6 which corresponds, according to Henry's 504 law, to 90mmHg, the usual oxygen mass fraction in arteries (see Ortiz-Prado et al. (2019)). crit , the 505 hypoxia threshold, is cell-line dependent, for tumor cells, it has been set to a very low value: 10 −6 506 (≈ 20mmHg, for common human tissue cells, hypoxic level is defined between 10 and 20mmHg 507 Khan et al. (2007)) The function ( ̄ , crit , env ) is plotted Figure 8A. 508 Function (1 −  ) represents the dependency on pressure: We set crit to 6kPa as initial guess (in Helmlinger et al. (1997), they found a inhibitory pressure 510 at 10kPa) and 1 , the pressure threshold when the inhibitory process starts, at 2kPa. The function 511  ( , 1 , crit ) is plotted Figure 8B. 512 The tumor oxygen consumption → follows the law with the consumption rate due to proliferation, 0 the consumption rate due to quiescent 514 metabolism and: The model does not discriminate between proliferating and quiescent cells, but the growth is sub-516 ject to ( ̄ , crit , env ). To make possible the comparison with the experimental proliferative cell 517 quantities, the following relationship has been set:  is also used in the definition of hypoxic necrosis rate which reads where = 0.01 is the necrotic growth rate. Its value has not been changed in this article. The mixed FE problem has been computed on 5 different meshes, with uniform cell sizes ℎ = 543 50, 20, 10, 5 and 2.5 m. To measure the FE solution degradation the primary variable ̄ , the oxy-544 gen mass fraction, has been monitored at the spheroid center during 4 days (see Fig Figure 10) .

545
The thinner mesh of cell size ℎ = 2.5 m has been used as reference for the RMSE. Despite an 546 important increase of the computation time, the mesh of cell size ℎ = 5 m has been chosen to 547 restrict the relative degradation of the FE solution to = 0.01 (see Table 2). With the characteristic time of the consolidation̄ , equal to 2 , sets to 100 m and , the consolidation coefficient: where and are Lamé constants of the solid scaffold, is its intrinsic permeability and the fluid 559 dynamic viscosity. The addition of the RMSE of the 4 samples at̄ = 0.01, 0.1, 0.5, 1 (see Fig.Figure 11) 560 with the analytical solution as reference gives ∑ = 0.0028. The surface error for different 561 cell sizes ℎ and time steps is in Fig.Figure 11(right). 562 Appendix C. Sensitivity analysis 563 The seven parameters of the model have been set to generic values, gathered Table 1  responding cost has been computed (see Figure 12 and Figure 13).
To analyse the influence of their gradient on the finite element solution, unit-norm auxiliary parameters have been considered:  To study the cell-line specific parameters, the experimental data of the free MCTS have been used. However, the only physical quantity given for this configuration is the MCTS volume. This sparse experimental data in front of a multiphase system have made difficult to set a well defined cost function. To build a defined and differentiable one, a simulation noted goal has been run freely until it reaches the experimental volume of days 1,2,3 and 4. At the corresponding iterations, the numerical quantities as porosity and tumor phase saturation have been stored to be used has goal. For the day , the tumor volume is equal to: Where Ω is the whole computation domain, as = 0 outside of the tumor zone.

573
Then a simulation with the same parameters, has been run for 4 days and its volume has been 574 compared to goal at the time steps corresponding to 1, 2, 3 and 4 days, noted ( = 1, 2, 3, 4) . One