Digital twinning of Cellular Capsule Technology: Emerging outcomes from the perspective of porous media mechanics

Spheroids encapsulated within alginate capsules are emerging as suitable in vitro tools to investigate the impact of mechanical forces on tumor growth since the internal tumor pressure can be retrieved from the deformation of the capsule. Here we focus on the particular case of Cellular Capsule Technology (CCT). We show in this contribution that a modeling approach accounting for the triphasic nature of the spheroid (extracellular matrix, tumor cells and interstitial fluid) offers a new perspective of analysis revealing that the pressure retrieved experimentally cannot be interpreted as a direct picture of the pressure sustained by the tumor cells and, as such, cannot therefore be used to quantify the critical pressure which induces stress-induced phenotype switch in tumor cells. The proposed multiphase reactive poro-mechanical model was cross-validated. Parameter sensitivity analyses on the digital twin revealed that the main parameters determining the encapsulated growth configuration are different from those driving growth in free condition, confirming that radically different phenomena are at play. Results reported in this contribution support the idea that multiphase reactive poro-mechanics is an exceptional theoretical framework to attain an in-depth understanding of CCT experiments, to confirm their hypotheses and to further improve their design.


32
As a tumor grows, it deforms surrounding living tissues, which causes strains and associated stresses. 33 Mechano-biology focuses on these mechanical forces and their interplay with biological processes 34 which has been extensively studied experimentally (Helmlinger et al. (1997)). Within this context, 35 current mathematical models of tumor growth are becoming more and more reliable, comple-36 ment experiments and are useful tools for understanding, explaining and building upon these

Methods and Model
= 34 m) monitored for 8 days ( Figure 1D); 145 The reliability of the mathematical model has been tested with the denoted validation data set this fictitious solid phase is two orders of magnitude lower than that of the alginate solid scaffold 185 ( Figure 1A). A unique physical model is defined for the three domains, with some penalty parame-186 ters (e.g. a low intrinsic permeability) avoiding cell infiltration in the alginate shell domain. Oxygen 187 advection-diffusion within the medium/interstitial fluid phase is considered, oxygen acting as the 188 limiting nutrient of TC, as prolonged hypoxia leading to the cell necrosis.
Computational framework 243 We implemented the above model in Python and C++ within the FEniCS framework Alnaes et al.
All the details of the process, auxiliary parameters and cost functions can be found in Appendix C.

265
For both configurations, the optimization procedure was based on sensitivity profiles, that is to say, . Sobol indices of the solution sensitivity on 7 parameters: the TC dynamic viscosity, the parameter accounting for the joint impact of the ECM thinness and cell surface tension, the TC growth rate, and 0 the oxygen consumption rate due to growth and quiescent metabolism respectively, 1 and crit , two thresholds which govern the pressure-induced inhibition of the TC proliferation. Free MCTS configuration,A, First order analysis: Only 5 parameters remain, the governing parameter is , the tumor cells growth rate, the sensitivity of the solution on the pressure parameters, 1 and crit , is 0. B Interaction: among 10 parameters tuples, one is significant ( , ) 14.5% of the solution variance. Thus, these two parameters are not independent and should identified together. The total variance of the solution shows that, considering all the interactions, the influence of each parameter alone is not qualitatively changed: from 81.1% to 66.4%, 0 from 7.4% to 6.1%, from 6.2% to 5.1%. Encapsulated MCTS configuration, C, First order analysis: the governing parameter is crit the inhibitory pressure of tumor cells growth (73.5% of the solution variance). D Interaction: the sum of 21 parameters tuples represents 3.6% of the solution variance (the detail of 21 tuples can be found in Appendix C. Sensitivity analysis, table 6).
we optimized the parameters that gathered 90% of the variance of the solution. The selected pa-

277
Based on a detailed sensitivity study, the optimized set of parameters was tested and cross-validated

298
For all parameters perturbations in the first order and second-order interaction analyses, the pres-299 sure of the TC phase = + was less than 1 KPa, thus the first threshold of growth inhibitory 300 due to pressure 1 was never reached and, a fortiori, the critical threshold of total inhibition crit .

301
Thus, the sensitivity of the FE solution to 1 and crit was 0. The 3 parameters , and 0 has there-302 fore been optimized for the free configuration. For the encapsulated configuration, the governing 303 parameter is the critical inhibitory pressure crit (first order crit = 73.5%, with interaction 70.9%). , 304 and 0 has already been identified for the free configuration, the only non negligible parameter 305 remaining is 1 (first order index 1 = 3.4%, with interaction 3.3%). The difference between Sobol 306 indices of first order and interactions is weak, indeed, the 21 parameters tuples gather only 3.6% of 307 the solution variance. Thus, in the encapsulated configuration the parameters can be considered 308 non-correlated and be identified separately.

309
Such results allow us to highlight that in the encapsulated configuration the mechanical constraint 310 is the phenomenon that determines the overall growth phenomenology provided by the mathe-311 matical model.

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

314
Mead simplex algorithm and fitted to the experimental data with a = 0.031. To be physically 315 relevant, the same parameters set should be shared by the two configurations. Hence, these three  parameters were injected within the encapsulated configuration and its two governing parameters 317 1 , crit were identified using the same algorithm. We fitted the experimental data of the encapsu-318 lation with a = 0.124. Figure 3A shows the two configurations fitted with the following set 319 of parameters: = 3.33 .10 −2 , 0 = 6.65 .10 −4 , = 890, 1 = 1432, crit = 5944 (see Table 1). This set is 320 cell-line specific, only relevant for CT26 mouse colon carcinoma. 329 Figure 3B shows that the modeling range on two thick capsules is in accordance with the experi-330 mental results. Two fits with an alginate stiffness at = 52.5kPa and = 70kPa respectively are 331 proposed.

332
The Figure 3C shows results relative to two thin capsules. The dynamics is properly reproduced 354 Figure 5A and B show the interplay between oxygen consumption and necrosis. Indeed as men-355 tioned in the experiments, 85 hours after confluence, the viable space remaining for TC is a 20 m 356 thick rim. This is explicit in Figure 4, upper right circle, NTC quarter. The comparison of Figure 5F 357 and B shows a relation between the saturation of TC and their necrotic fraction. This is a basic 358 experimental fact that, when the cells bodies collapse in a necrotic core, the aggregate density 359 increases accordingly. Figure 5D and E allow to 'visualize' the overall dynamics of the process:  analysis (see Appendix C. Sensitivity analysis eq.39). At the confluence time, is importantly lower 369 than the pressure in the tumor cell phase since at that time the MCTS consists also of 40% of IF.

370
After confluence the saturation of tumor cells increases progressively, so , becomes closer to the 371 pressure sustained by the cells.

372
In the presented physics-based approach, mass conservation is prescribed, so the growing MCTS, 373 which increases in density and size, have to lead to a decreasing mass of interstitial fluid. This re-374 sult, which cannot be measured experimentally, is shown in Figure 5 where a sucking phenomenon 375 due to IF absorption by growing TC can be observed. The Figure 5C shows that after confluence 376 the interstitial fluid pressure becomes positive during a while (see plot relative to 28h). Indeed, 377 after confluence the initial gradient of IF pressure (green line in Figure 5C) reverses since cells in 378 the proliferative peripheral areas move toward the core so IF has to go in the opposite direction, identified on a training data set (reported in Figure 1C and Figure 1D), has been used for all the 398 numerical simulations performed.

399
In the frame of the parameter identification process a second order sensitivity analysis has re-400 vealed that the parameters of the model become almost independent under confinement (see 401 Figure 2D). Results of sensitivty analyses also demonstrate, that if the tumor is free to grow the      According to the different phases, solid scaffold , medium/interstitial fluid and tumor cells phase 505 , described in the main article, section Methods, which constitute the multiphase system, at each 506 point in the domain, the following constraint must be respected where is the volume fraction of phase . Defining the porosity as Mass conservation equations 512 We express the mass conservation equation for each phase. We use a material description for the 513 motion of the solid phase and a spatial description for the fluid phases, whose reference space 514 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 where is the death rate of tumor cells. Note that only one of Eqs 11-12 is independent: actu-531 ally, one can be obtained subtracting the other from Eqn 8 and accounting for the constraint Eqn 532 10.

533
Oxygen is the only nutrient which we consider explicitly. Another mass balance equation is intro-534 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 536 by tumor cells due to their metabolism and proliferation rate.

537
Momentum conservation equations 538 We neglect here the effect of gravitational body forces as their contribution is negligible compared 539 to that of other forces. Furthermore, as we assume quasi-static processes and small difference in 540 density between cells and aqueous solutions, inertial forces and the force due to mass exchange 541 can also be neglected. These assumptions simplify the general form of the linear momentum 542 balance equation given in Gray and Miller (2014) which becomes where ̄ is the stress tensor of phase , is the set phases connected to and → is the inter-544 action force between phase and the adjacent phases. Summing eqn 14 over all phases gives the 545 momentum equation of the whole multiphase system as where ̄ is the total Cauchy stress tensor acting on the multiphase system.

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

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

564
However, the overall multiphase system is not incompressible, because of the presence of porosity 565 that evolves according with the scaffold deformation. As phases are incompressible, their densities 566 with = , , are constant. As the solid phase is quasi incompressible, the Biot's coefficient is 567 set to unity. With these premises, the total Cauchy stress tensor appearing in eqn 15 is related to 568 the Biot's effective stress as follows where = + is the so-called solid pressure, describing the interaction between the two 570 fluids and the solid scaffold.

575
the Young modulus of the solid scaffold and its Poisson ratio.

576
The experimental measurement of cells density inside the capsule revealed a strong dependency 577 to necrotic fraction ̄ . Hence, the pressure-saturation closure relationship has been improved 578 with respect to that proposed in Sciumè et al. (2014a), to be more physically relevant and adapted 579 to confinement situation with pressure difference between tumor and interstitial fluid (i.e. = − ). The saturation is 581 directly linked to the partial pressure of the phase and a constant parameter , which accounts for 582 the effect of cell surface tension and of the refinement of the porous network (see Sciumè et al. 583 (2014d) for the biophysical justification of the proposed equation). Its influence is offset by the 584 necrotic fraction of tumor cells, ̄ (see Fig.Figure 7), which allows us to model necrotic areas of 585 very high cell density according with experimental evidence.

597
Tumor cell growth is related to the exchange of nutrients between the IF and the living fraction of 598 the tumor. The total mass exchange from IF to the tumor cell phase is defined as Note that 1 − ̄ is the living fraction of the tumor. is the tumor growth rate parameter, 600 cell-line dependent.  and  are regularized step functions varying between 0 and 1, with two 601 threshold parameters 1 , 2 , that is to say  =  ( , 1 , 2 ). When the variable is greater than 2 , 602  is equal to 1, it decreases progressively when the variable is between 1 and 2 and is equal to 603 zero when the variable is lower that 1 .  represent the growth dependency to oxygen:  Figure 8A. 609 Function (1 −  ) represents the dependency on pressure: An example of the function  ( , 1 , crit ) is plotted Figure 8B, we have set crit to 6kPa as initial 611 guess (in Helmlinger et al. (1997), they found a inhibitory pressure at 10kPa) and 1 , the pressure 612 threshold when the inhibitory process starts, at 2kPa.

613
As tumor grows, nutrients are taken up from the IF so that the sink term in eq.13 takes the following 614 form: Nutrient consumption from IF is due to two contributions: the growth of the tumor cells, as given by 616 the first term within the square brackets in eq.29, the metabolism of the healthy cells, as presented 617 in the second term. Thus, is related to the cell proliferation, as discussed above; whereas the 618 coefficient 0 relates to the cell metabolism. is an adaptation of the previous step functions for 619 the cell metabolism: Appendix 0 Figure 8. Two mechano-biological laws. A ( ̄ , env , crit ). The TC growth and nutrient consumption are dependent to the oxygen mass fraction ̄ . If it is lower than crit , the TC growth is stopped and the nutrient consumption is reduced to the metabolism needs only. If it is greater or equal to env , the growth and the nutrient consumption are maximum. B  ( , crit , 1 ). The TC growth and nutrient consumption are dependent to the TC pressure. If it is greater than 1 , the 2 processes begin to be strongly affected and if the TC pressure reaches crit , they are totally stopped.
The model does not discriminate between proliferating and quiescent cells, but the growth is sub-621 ject to ( ̄ , crit , env ). To make possible the comparison with the experimental proliferative cell 622 quantities (see Figure 6), 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. As the experimental data on necrosis were to sparse 626 for this parameter identification (only a few stained cells imaging), we have kept its generic value.   Fig.Figure 9). To avoid this phenomenon, the composite Taylor-Hood element has been set to  The mixed FE problem has been computed on 5 different meshes, with uniform cell sizes ℎ = 658 50, 20, 10, 5 and 2.5 m. To measure the FE solution degradation the primary variable ̄ , the oxy-659 gen mass fraction, has been monitored at the spheroid center during 4 days (see Fig Figure 10) .

660
The thinner mesh of cell size ℎ = 2.5 m has been used as reference for the RMSE. Despite an 661 important increase of the computation time, the mesh cell size of ℎ = 5 m has been chosen to 662 restrict the relative degradation of the FE solution to = 0.01 (see Table 2).  The fluid is free to escape only at the loaded boundary, this boundary condition is known as drying 672 condition. The analytical solution of this problem is: 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 674 dynamic viscosity. The addition of the RMSE of the 4 samples at̄ = 0.01, 0.1, 0.5, 1 (see Fig.Figure 11) 675 with the analytical solution as reference gives ∑ = 0.0028. The surface error for different 676 cell sizes ℎ and time steps is in Fig.Figure 11( Table 1). Experimental, green dot ; Model, red ; Sensitivity evaluation, blue x. Left: free growth ; Right: encapsulated growth Appendix C. Sensitivity analysis 678 For the sensitivity analysis, the experimental input data were:

679
• for the free MCTS, the volume monitored over a time span from day 1 to day 4. These data

685
We performed a variance-based sensitivity study of the FE solution on the parameters, both on the 686 free and encapsulated MCTS, as follows:

687
• A first order analysis, the 7 parameters are disturbed one at a time respectively to a 8 points 688 grid.

689
• Interaction analysis, the 21 parameters tuples are evaluated at the 2 extreme points of the 690 grid.

691
All the results were interpreted with a polynomial model in order to quantify their weights in the In order to build a the cost function for the multiphase system with the sparse data available, a simulation noted goal has been run freely until it reaches the experimental volume corresponding to exp free , the experimental data of days 1,2,3 and 4. At the corresponding iterations, the numerical quantities have been stored. At day , the tumor volume is equal to: where Ω is the whole computation domain, as = 0 outside of the tumor zone.

695
A second simulation with the same parameters has been run for 4 days and its volume has been 696 compared to goal at the time steps corresponding to 1, 2, 3 and 4 days, noted ( = 1, 2, 3, 4) . One where conf is the internal pressure, is the Young's modulus, and in and out are the inner and The numerical approximation ( , Φ, ) has the two corresponding components ( , ) in Capsule . 710 We compared the FE solution with the experimental data 1 day after confluence. One can write 711 this cost function explicitly 2 : Two simulations were run with the 7 parameters at their generic value (see Figure 12). We denoted 713 the respective cost functions 0 free and 0 conf . Appendix 0 Table 3 In order to quantify the impact of each parameter, the following linear model was set: where is an auxiliary parameter ∈ [−1, +1] representing the perturbations of the ℎ parameter 720 along the grid and the slope of the variation.

721
In a first order analysis, the influence of the ℎ parameter is given by the Sobol indices: The results for the free and encapsulated configurations are reported in Tables 3 and 4.

724
As the independence of physical phenomenons involved in both configuration is one our major 725 modeling assessment, the interaction between parameters has also been studied. The 21 tuples with the respective Sobol indices: