## Abstract

Characterization of tumor growth dynamics is of major importance for cancer understanding. By contrast with phenomenological approaches, mechanistic modeling can facilitate disclosing underlying tumor mechanisms and lead to identification of physical factors affecting proliferation and invasive behavior. Current mathematical models are often formulated at the tissue or organ scale with the scope of a direct clinical usefulness. Consequently, these approaches remain empirical and do not allow gaining insight into the tumor properties at the scale of small cell aggregates. Here, experimental and numerical studies of the dynamics of tumor aggregates are performed to propose a physics-based mathematical model as a general framework to investigate tumor microenvironment. The quantitative data extracted from the cellular capsule technology microfluidic experiments allow a thorough quantitative comparison with in silico experiments. This dual approach demonstrates the relative impact of oxygen and external mechanical forces during the time course of tumor model progression.

## INTRODUCTION

Although tumor growth is indisputably controlled by biochemical factors and genetic alterations, it is now well accepted that physical forces play a pivotal role in tumor development (*1*, *2*). From a physical viewpoint, tumor behavior is expected to obey conservation laws of energy, mass, and momentum, like all other physical systems studied in science and engineering. A physical and quantitative approach for modeling of tumor behavior should provide an original picture of growth and invasion mechanisms that could bring complementary insight as compared with purely biochemical or genetic approaches (*3*). Many numerical physics-based models have emerged to describe the kinetics of tumor growth at the tissue and organ scales (*4*–*8*). Nevertheless such models remain of limited predictive interest for a systematic clinical use (*9*). Many reasons may explain this pitfall: Tumor expansion is governed by a number of coupled biological and physical processes occurring at different spatial scales and often having very different characteristic time; patient-specific parameters are too difficult to assess, etc. Within this context, controlled in vitro models are instrumental to pinpoint the physical processes at play during tumor growth.

The influence of mechanical stresses exerted by the environment on tumor progression was first quantitatively demonstrated in a pioneering experiment by Helmlinger *et al.* (*10*). By embedding in vitro tumor models, namely, multicellular spheroids, in an agarose matrix, they could derive an estimate for the growth inhibitory stress (5 to 15 kPa). This threshold value was then theoretically formalized and termed “homeostatic pressure” that corresponds to the pressure at which cell apoptosis and division rates, both affected by stress in an opposite fashion, balance each other. Following the same strategy, some of us (K.A., B.G., and P.N.) have devised a new approach to form multicellular tumor spheroids (MCTS) and monitor their growth inside hollow porous elastic capsules of alginate gel (*11*). This methodology has been named cellular capsule technology (CCT). By contrast with cells embedded in a gel, spheroids first grow freely until they reach three-dimensional (3D) confluence inside the capsule. At this stage, further spheroid growth dilates the capsule, which conversely exerts a restoring mechanical stress onto the spheroids. Inverse analysis of the measured strain of capsule shells allows for the evaluation of the MCTS internal pressure. In both experimental procedures reported above (*10*, *11*), growth inhibition and formation of a necrotic core were observed and hypothesized to arise from pressure buildup. Nonetheless, because of confinement growth condition, cell density also drastically increases; hence, it remains unclear whether the MCTS core is necrotic because of restricted oxygen diffusion through the denser tissue or because of mechanical stress only. In the case of a hypoxia-induced necrotic core and growth inhibition, the measured pressure would not be the cause but only the consequence of confinement. To test this hypothesis, a new experimental configuration has been investigated here. It consists of a cylindrical hollow capsule, generated by an adaptation of the CCT and which alleviates the effect of compaction by releasing confinement along the main direction of the tubular capsule.

These two experimental approaches not only allow a qualitative investigation of the potential coupling between cell density and pressure, but they also provide two simple tumor cell aggregate models that can be quantitatively investigated theoretically. To do so, we have developed a biophysical mathematical model to reproduce in silico the two experimental configurations. Fluid analogy is used to simulate the behavior of cell aggregates (*12*). The experimental protocol permits a rigorous control of input parameters and an efficient decoupling of biological to mechanical effects on cells. Therefore, it constitutes a relevant benchmark for tailoring the numerical model and to better gain insight into the biomechanical mechanisms at play. Numerical simulations are expected to give valuable information not measurable in the laboratory by relying on robust experimental data. Wealth and accuracy of the experimental data allowed us to validate and critically test and refine the mathematical model, which could be used as a general framework to investigate other systems, both in vitro and in vivo.

### The experimental protocol

Multicellular tumor spheroids (widely used as in vitro 3D tumor models) are embedded within porous alginate capsules allowing to mimic generation of mechanical cues arising from the host environment of the tumor. Alginate is used because it forms a porous hydrogel, with a pore size of the order of 20 nm, thus allowing diffusion of oxygen and nutrients through the alginate shell while retaining cells entrapped.

Cellular capsules are prepared using a microfluidic coextrusion device that consists of three coaxial glass capillaries (see Fig. 1A) (*11*). The cell suspension (CS) flows into the innermost capillary; an alginate solution (AS) is injected into the outer capillary; and the intermediate capillary contains a calcium-free culture medium (CM) that avoids alginate gelation inside the chip. Coextrusion is performed in the air at high flow rates (~50 ml/hour) to generate a liquid jet. However, this jet is fragmented downstream into droplets due to Plateau-Rayleigh instability. Upon contact with a calcium chloride bath, the composite droplets readily cross-link as shells encapsulating cells (Fig. 1A and fig. S1A). By dipping the nozzle of the chip directly into the calcium bath, surface tension mismatch is cancelled out, and the capillary instability does not occur any longer. The composite liquid cylinder cross-links and leads to the generation of centimeter-long alginate tubular capsules that enclose cells (Fig. 1B and fig. S1B). Reduction in the flow rates also ensures more regular wall thickness (*13*).

By adjusting the flow rate inside the capillaries, the shell thickness, and thus the aspect ratio of the spherical and tubular capsules, can be tuned (*11*). Once the capsules are transferred into the culture medium at 37°C, tumor cells (here, CT26 colon carcinoma cells) are in standard conditions to go through their cell cycle. In spherical capsules, they usually form individual spheroids (movie S1), while several aggregates are most often found to be dispersed along the length of the tubular capsules (fig. S1B). While they are observed to eventually merge into long cellular cylindroids, we will only focus on the growth process of individual and isolated aggregates (see fig. S1C and movie S2) and discard the coalescence process. Growth kinetics is then monitored by time-lapse phase-contrast microscopy, and a home-made contour detection software is used to derive the average geometrical configuration of the aggregate at different times. Fluorescence confocal microscopy was also used with stable fluorescently labeled cell lines when required. When the growing cellular aggregate fills up the spherical capsule (corresponding to confluence and subsequent confinement) or when a cell aggregate comes in contact with the inner alginate wall of the tubular capsule, it starts to interact mechanically with the capsule and to deform it. As a result of the capsule-cellular aggregate mutual interactions, the internal pressure of the aggregate increases. Pressure quantification is experimentally achieved from the measurement of the displacement of the inner alginate wall and the knowledge of the Young’s modulus alginate hydrogel (*E* = 68 ± 21 kPa) (Fig. 1C) (see Materials and Methods). Detailed cell culture, encapsulation, microscopy, and image analysis procedures are given in Materials and Methods.

The spherical configuration was used to investigate the impact of confinement on the growth of MCTS in a closed compartment. Tubular capsules allow the release of the confinement condition in the longitudinal direction. The cell aggregate can grow freely along the tube axis. This second configuration is expected to reduce the relative importance of the pressure on growth dynamics and, consequently, to highlight the effect of hypoxia, as compared with the spherical geometry. Analysis of both experimental model systems is complemented with mathematical physics-based modeling to assess the combined effect of pressure- and confinement-induced hypoxia on cell proliferation rate.

### The mathematical model

Considering tumor cell dynamics from a biophysical perspective permits to describe its evolution in terms of conservation equations, which must be complemented by constitutive relationships encompassing biological aspects like cell division and metabolism. When the alginate capsule is formed, the content is a mixture that is mostly composed of culture medium (*m*) and tumor cells (*t*) as shown in Fig. 2A. Oxygen and other nutrients (*n*) are also present in the mixture with a minor mass fraction with respect to dominant species *t* and *m*. This mixture has a liquid nature and obeys conservation laws, namely, mass, momentum, and energy conservation equations for the mixture and each mixture species (*14*). Evolution of mass fraction of the main components, *t* and *m*, is modeled using the Cahn Hilliard theory (*15*). Mass conservation of a species includes the advective and the diffusive species transport mechanisms. In the present case, once capsule is formed, mixture velocity is extremely low and is assumed to have a negligible impact on the overall behavior of the system. Hence, diffusive transport is only considered with its driving force originating from the variation of mixture chemical potential.

Heterogeneity of species mass fraction within the capsule allows to identify the area occupied by cell aggregate as that where local mass fraction of tumor cells, ω* _{t}*, is equal to or higher than a certain equilibrium value,

Following a standard description of multicellular spheroids (*12*, *16*–*19*), the cellular aggregate is modeled as a viscoelastic fluid characterized by a surface tension and rheological parameters, namely, a non-Newtonian dynamic viscosity and a relaxation time. Since mixture velocity is assumed negligible, cell aggregate viscoelasticity does not affect the mixture behavior. Hence, surface tension is the sole physical property that, together with growth rate and metabolism, influences the numerical solution. Surface tension, σ, and local mass fraction of tumor cells, ω* _{t}*, appear explicitly in the mixture chemical potential (

*20*)

*f*(ω

*), bulk free energy density of the mixture, has the following form*

_{t}The equilibrium value is set at the density of a free spheroid, * _{t}* = 0, and the other at

*21*) and not as a pure phase. In Eq. 1, the first right-hand side term accounts for cell-cell attraction or mutual repulsion and tends to push the system toward its equilibrium configuration (see Fig. 2C), whereas the second term originates from the gradient part of the free energy density. Details about the assumed form of the free energy density and derivation of the potential are reported in Materials and Methods.

The chemical potential times the gradient of tumor cells mass fraction leads to a body force that accounts for the capillary forces induced by the surface tension of the cell aggregate. Such body forces introduced in the momentum conservation equation of the mixture allows computing the mixture pressure within the cell aggregate before and after confluence.

In the mass conservation equation of the tumor cell species, a reaction term allows to account for cell division. Oxygen is assumed to be the primary “nutrient” species with the most drastic impact on tumor growth, metabolism, and necrosis: Cells need a minimum amount of oxygen to survive and an extra amount to proliferate (*4*, *22*). Consequently, accounting for mass conservation equation of oxygen allows us to compute its mass fraction within the domain. Cell growth rate is assumed to depend on oxygen concentration at the considered point. If the necessary resources for mitosis are not available, cells cannot replicate and remain in a quiescent state.

The decrease in cell activity does not only depend on the local oxygen concentration. Experimentally, the increase in cell aggregate compactness is observed to lead to an increase in the internal pressure, which subsequently reduces the proliferation rate. We integrate these observations into the mathematical model by introducing a pressure threshold in the order of 2 kPa, above which mechanical stress inhibits cell division. Growth rate follows a functional dependency with oxygen and pressure qualitatively depicted in fig. S2 (see Materials and Methods for more details).

The resulting system of nonlinear partial differential equations (PDEs) (momentum conservation equation of the mixture and the mass balance equations of cells and oxygen species) can be applied to cell aggregate before and after confluence and, thus, fully describes the impact of confinement on cell aggregate growth. The pressure of the mixture, the mass fraction of tumor cells, and the mass fraction of oxygen are the primary variables of the proposed physics-based mathematical model.

In this work, for the sake of simplicity without loss of generality, we restrict our numerical analysis to relatively thick-shelled spherical and tubular capsules to ensure that shell deformation is small with respect to the characteristic length of the capsule (i.e., the radius), which allows us to solve the system of PDEs in a fixed geometry. On the other hand, from an experimental viewpoint, shell deformation needs to be detected and measured as accurately as possible, so the capsule cannot be too thick; we found that an optimal compromise condition was obtained for *t*/*R*_{i} ≈ 0.3, where *t* is the alginate wall thickness and *R*_{i} is the inner radius of the capsule. As a consequence, alginate strain and the confinement release induced by the spherical capsule deformation are not accounted in the numerical simulation and so should be negligible to warrant consistency between the mathematical model and the experimental one. To ensure this hypothesis, the inner radius of the spherical capsule is systematically retrieved a posteriori (see Materials and Methods for details). As long as this numerical estimate matches the experimental data, it ensures that the coupling between spherical capsule inner radius deformation and tumor cell growth dynamic has a relatively weak impact on the overall observed behavior, and so, the experiment remains in the scope of validity of the model. It is also interesting to highlight the dual relationship existing between the experimental and computational procedures. In the experimental model, pressure increase is calculated by postprocessing of measured capsule deformation; conversely, in the mathematical model, aggregate pressure comes out from the thermodynamic description of the system and is subsequently used to compute capsule strain.

In addition, boundary conditions (BCs) set by the experimental capsule system need to be defined so that they are physically consistent with the assumption and equation coupling of the mathematical model. First, we exploit the axial symmetry of the capsule systems. Second, the geometries of capsules and tubes are also symmetric with respect to the plane going through the origin of the cylindrical coordinate reference system and orthogonal to the *z* axis. Considering this plane of symmetry allows axisymmetric analysis of one-half of the experimental system (see Fig. 3, A and B). These geometrical symmetries are very convenient from the computational viewpoint because they reduce importantly computational costs of the mathematical model. Symmetry conditions are assumed at the *z* axis and at the middle symmetry plane (BC1 and BC2 in Fig. 3). At the remaining boundaries of the intracapsular domain (BC3 in Fig. 3), we make the following hypotheses: (H.a) The mass fraction of tumor cells is equal to zero; (H.b) the gradient of the chemical potential has no normal component; (H.c) oxygen transport is associated with convective-type BCs (diffusive flow proportional to the difference between mass fraction of oxygen at the intercapsular bound and the mass fraction of oxygen in the external medium); and (H.d) the mechanical stress is fixed and equal to the hydrostatic pressure in the external medium.

(H.a) and (H.d) imply that even after cell aggregate confluence, a very thin film of medium persists at this boundary (BC3 represents the physical limit between the intercapsular space, where tumor cells and medium coexist in the mixture, and the porous alginate domain, which contains only medium); (H.a) together with the assumption of the absence of normal gradient of the chemical potential, (H.b), ensures that tumor cell species do not “escape” the intracapsular/tubular domain. Note that these conditions do not exclude advective flow of the medium at BC3, which is indeed permitted in the experimental setup. Convective-type condition for oxygen (H.c) allows to take into account a resistance to oxygen diffusion across the spherical/tubular capsule porous walls. This is more precisely modeled by the identification of a convective-type exchange coefficient that integrates the thickness and transport properties of alginate walls. Last, with (H.d), it is assumed that the difference between internal capsule pressure at BC3 and the external medium pressure is not notable despite the hydraulic resistance of the alginate shell.

Material and geometrical parameters used for the numerical simulations are reported in Materials and Methods together with the governing and constitutive equations of the model.

## RESULTS

To gain insight into the mechanisms that govern cell aggregate growth, experimental findings and numerical results are quantitatively analyzed and compared. The same set of parameters is used for both the spherical and tubular configurations (see Materials and Methods for values). As mentioned above, we consider that the growth of the multicellular aggregate is primarily regulated by mechanical pressure and oxygen mass fraction, which are initially set equal to 0 Pa and 9 × 10^{−6}, respectively. To account for hypoxia-induced effects, we assume that growth rate decreases when oxygen mass fraction drops below a threshold value

For the spherical capsule, the overall behavior is drastically different depending on whether cells grow freely or in stressed conditions. Before and after confluence, results are analyzed separately in the following two paragraphs. For tubular capsules, we are interested in the growth of cell aggregates of initial size of the order of the tube diameter and in their subsequent deformation. These results are presented in a distinct third paragraph.

### Growth of a cellular aggregate in a spherical capsule: Preconfluent stages

Figure 4A shows the radius and pressure of a multicellular aggregate growing within an alginate capsule as a function of time. Confluence takes place at *t*_{c} = 66 hours. The inner radius of the alginate capsule is 90 μm, the wall thickness is 30 μm, and the radius of the cell aggregate is around 40 μm at *t* = 0 (experimentally taken arbitrarily at the beginning of image acquisition). Capsule radius was selected to be smaller than the diffusion length of oxygen to ensure optimal penetration of oxygen inside the spheroid (*11*) and thus avoid the effect of oxygen shortage on cell proliferation rate. In Fig. 4A, the evolution of the tumor aggregate radius is depicted together with numerical results (see also movie S3). The agreement is fine and corresponds to a progressive increase in the radial diffusive flow (see fig. S3A) and an overall exponential growth regime of the aggregate volume (see fig. S4, spheroid capsule case). In this preconfluent stage, kinetics of encapsulated cell aggregates is very similar to that of a freely growing spheroid (see Fig. 4A). This indicates that oxygen penetrates effectively the alginate capsule during this phase. A few hours before confluence, a slight deviation from the free growth behavior can be detected, as magnified in the upper left insert of Fig. 4A. MCTS growth is slightly affected before the aggregate physically hits the capsule wall. A crossover phase between the free and confined growth regimes that extends over a few hours around confluence is observed both numerically and experimentally. Numerically, at *t* ~1 day, oxygen mass fraction in the central area of the cell aggregate drops below *t* = 24 and 54 hours).

The pressure difference between cell aggregate and surrounding medium emerges from the fact that cell spheroids are characterized by an effective surface tension. According to the Young-Laplace equation of capillarity, the numerical model provides a pressure difference decreasing as 1/*R*. This is observable in Fig. 4E: The tumor aggregate pressure at *t* = 54 hours is slightly lower than at *t* = 24 hours. However, this preconfluent variation is not measurable experimentally and almost invisible compared with the pressure increase amplitude in the postconfluence stage (Fig. 4A).

### Growth of a cellular aggregate in a spherical capsule: Postconfluence stages

As shown in Fig. 4A, postconfluence growth of a multicellular spheroid is experimentally characterized by a kinetics that exhibits two regimes: In the first one, which lasts about 24 hours, an average velocity, *v* = 2 μm/day, which is detected if recorded over long times with sufficient temporal resolution. This residual growth was suggested to be driven by a remaining cell activity (i.e., proliferation) inside the alginate capsule in the peripheral rim of the spheroid associated with a delay in the lysis of dying cells in the core of the aggregate (*11*). We obtain an excellent agreement between experimental and numerical data. As shown in Fig. 4E (*t* = 84 hours), we may remark that (i) the pressure inside the spheroid has reached 2086 Pa, which is higher than the threshold value

To gain additional insight in the detailed mechanism, we shall now focus on the evolution of the cell mass fraction. Intuitively, by contrast with free growth, there is an increase in the cell mass fraction upon spatial confinement of the cell aggregate (because cell division is not instantaneously impeded throughout the whole aggregate). Hence, residual cell division in the bordering areas induces a relevant augmentation MCTS cell density accompanied by a relative slow increase in radius. As shown in Fig. 4C, the mathematical model predicts an overall increase in the tumor cells mass fraction ω* _{t}*. from the equilibrium value 0.65 to 0.86. This increase in ω

*(Fig. 4A, bottom right insert) induces the rapid rise in MCTS pressure (Fig. 4, A and E).*

_{t}Last, let us examine the evolution of the third main variable of our model, that is, the oxygen mass fraction ω* _{n}*. Qualitatively, since confinement is accompanied with a drastic increase in the tumor cell mass fraction and that peripheral cells remain proliferative, we expect a significant increase in oxygen consumption inside the capsule. In addition, from a purely biophysical viewpoint, although the spheroid radius remains lower than the “nominal” oxygen and nutrient diffusion length [~300 μm in the absence of any glucose limitation (

*23*)], the increase in the cell density may also affect the penetration of nutrients as the cell aggregate forms a more tortuous porous medium. Whereas this aspect was overlooked in the original experimental paper (

*11*), the current numerical work allows us to assess its relevance and the coupled contribution to the other parameters. Conversely, Fig. 4D shows that oxygen concentration increases in the capsule once confluence is reached (between

*t*= 54 hours and

*t*= 84 hours). A rationale for this observation is that cell growth inhibition induced by pressure increase leads to a decrease in oxygen consumption and, thus, to an increase in the internal mass fraction of oxygen. Because of high internal pressure, most of the cells located in the core of the spheroid are in a quiescent state, which required a minimal amount of oxygen to maintain a metabolic activity. Model predictions are that final oxygen concentration may reach 90% of the initial value, which suggests that the impact of hypoxia on proliferation is marginal.

As shown in Fig. 4A, the multicellular spheroid internal pressure derived from capsule deformation measurements reaches ~3 kPa in the late stages. However, as mentioned above, there are two regimes; the pressure first increases from 0 Pa to 2 kPa within 1 day and more slowly in the later stages. This temporal pattern is reproduced numerically (see Fig. 4A): When the cell mass fraction ω* _{t}* exceeds the equilibrium value (

### Growth of a cellular aggregate in a tubular capsule

For long and open-ended cylindrical capsules containing cells, end effects can be neglected, and confinement is released along the axis of the tube. Tumor cell aggregates have the shape of cylindroids and are characterized by the end-to-end length after monitoring by optical microscopy (movie S2 and fig. S1C). To derive the pressure exerted by the cellular cylindroid on the alginate walls, we also monitored the variation of the tube radius at different positions as a function of time during growth (fig. S6). The experimental data are compared with the in silico model (Fig. 5A and movie S4). In early stages, the initial tumor cell aggregate, which has a cigar shape of relatively low cell density due to random aggregation of cell biased by the cylindrical geometry, undergoes a significant (by 50%) and rapid (within ~5 to 10 hours) contraction. In the context of the developed mathematical model, one may assume that the initial tumor cell mass fraction is below the equilibrium cell mass fraction and infer that this initial chemical potential unbalance drives the system toward the equilibrium configuration, *t* = 0 hours, which is beyond the scope of the present work. Following contraction, the growth of the cellular cylindroid is observed to be first exponential before becoming linear after about 2 days. As depicted in Fig. 5A, numerical calculations are in excellent agreement with the experimental data.

Despite the partial release of confinement, which could a priori allow the multicellular cylindroid to grow freely along the main axis of the tube without undergoing any compaction, we counterintuitively observe that cell mass fraction drastically increases since it even exceeds the equilibrium value *t* = 50 hours.

Figure 5A shows that a good agreement is found between numerical calculations and experimental data until *t* = 50 hours. Significant deviation is however observed in the later stages: The pressure value derived from the experimental data tends to level off, suggesting that cell proliferation rate and diffusion velocity are balanced in a quasi-steady state, while the value of the pressure derived from the mathematical model continues to increase. This discrepancy may be due to the fact that the deformability of the alginate tubes was not explicitly considered in the mathematical model (as explained in the “The mathematical model” section). As shown in Fig. 5E, the pressure felt by the cells depends on the position of the cells along the *z* axis (Fig. 5C). Closer inspection shows that pressure exceeds the threshold value *24*).

At first sight, the concentration of oxygen predicted by the model shows a reasonably good oxygenation of tumor cells, with an oxygen mass fraction between 7.5 × 10^{−6} and 9.0 × 10^{−6} for the whole time window that we explored (see Fig. 5D). Note, however, that some pieces of the cell cylindroids do experience an oxygen mass fraction slightly lower than the threshold value

## DISCUSSION

Growth of cell aggregates in hollow tubes and spheres was performed with the same cell lines. Yet, the main features for the growth kinetics were significantly different. Qualitatively, we may argue that the differences in the confinement conditions are the main determinant. In both experimental configurations, layered structures with different growth kinetics have been identified once the aggregate reaches a threshold size and become confined in the capsules. To push further our understanding of the underlying mechanisms governing cell aggregate growth under mechanical stresses, we performed a thorough and quantitative analysis by developing a mathematical model that relies on a limited set of variables, namely, the cell mass fraction, the oxygen mass fraction, and the internal pressure. Using the same computational model for both experimental conditions, we could reproduce the two growth kinetics and gain a predictive value on cell aggregate behavior within the capsules. These results are promising for investigating dynamics of small tumor cell aggregates in configuration difficult to create or monitor experimentally. This in silico approach also allowed to exclude the hypothesis of an oxygen-induced inhibition in spherical capsules and quantitatively assesses the impact of hypoxia in tubes. We could distinguish two confinement regimes: A strong confinement effect in spherical capsules, for which confluence is sudden and results in a homogeneous pressure rise above the inhibitory pressure value, and a weak confinement in tubes, where pressure continuously increases in the more central zone (i.e., older part) of the cylindrical aggregate, whereas tips (i.e., the younger part) remain proliferative. This distinction is of primary importance since, as a direct consequence of the confinement effect, a switch in the phenotype of peripheral cells has been clearly identified in the case of alginate spherical capsules (*11*): After prolonged confinement, most of peripheral cells become highly motile with a characteristic migration velocity of the order of 5 μm/hour. On the contrary, the cylindrical configuration does not display these features. Although further investigation is required to unveil the mechanism of phenotype switch, the differences in the confinement regime might be one of the causes.

If tumor cell growth can be predicted from the developed mathematical formulation in controlled environments, some additional developments must be achieved to use the model for more complex configurations, such as in vivo simulations. In our model, surface tension of tumor cells plays a leading role, controlling pressure inhibition and cell motion. The value of surface tension identified by comparing numerical results with experimental data is in the low range of commonly used values (see Table 1). This may suggest that surface tension is varying during the course of the experiment: It could be linked to the mechanical state of cells, as proposed experimentally for other cell properties (*25*). Upon pressure increase, the aggregate internal energy rises, and the system becomes largely out of equilibrium. Biophysical mechanisms could mitigate this condition by decreasing aggregate internal pressure, thanks to cell surface tension reduction. A quantitative study must be performed to clearly identify these processes and implement them in our mathematical formulation. Also, mobility of cells could be investigated more deeply to better characterize cell accumulation effect. The merging of cell aggregates, both in constrained and unconstrained environments, may be studied at the light of the developed mathematical model to evaluate reference values for cell mobility. Oxygen measurements inside the multicellular spheroids and cylindroids would also be useful to critically test the numerical simulations, although they are not easy to implement in a noninvasive manner. Last, an extension of both the experimental configuration and numerical model could be used to produce and analyze vascularized cellular aggregates constructed around one central engineered blood vessel [as developed in (*13*)].

## MATERIALS AND METHODS

### Cell lines and cell culture

We used wild-type mouse colon carcinoma CT26 cells (purchased from the American Tissue Culture Collection, ATCC CRL-2638). Cells were maintained in Dulbecco’s modified Eagle’s medium (DMEM; Invitrogen) supplemented with 10% fetal bovine serum (Invitrogen) and antibiotics [streptomycin (100 μg ml^{−1}) and penicillin (100 U ml^{−1}) (Gibco BRL)] in a humidified atmosphere containing 5% CO_{2} at 37°C with medium changed every 2 days. Cells were grown as subconfluent monolayers to prepare the CS used for encapsulation.

### Cell encapsulation

The coextrusion device was fabricated using a 3D printer (Envisiontec Micro Plus HI-Res, Earketyp 3D, France) using the computer-assisted design file provided in (*26*) and following the same procedure. A Teflon tip or a hydrophobic forged glass capillary was glued to the end of the 3D-printed device as a nozzle for spherical capsule production. As described in (*13*), this extra tapered nozzle was unnecessary for tube generation. Sterilization was performed instead by rinsing the chip with Biocidal ZF from Biovalley. The three fluid phases [CS, culture medium (CM), and alginate solution (AS)] were loaded into syringes (10MDR-LL-GT SGE, Analytical Science) mounted on computer-controlled pumps (neMESYS low pressure module V2). AS was prepared by dissolving 2.5% (w/v) sodium alginate (FMC, Protanal LF200S) in water and by adding 0.5 mM SDS surfactant (VWR international). The solution was filtered at 1 μm (Pall Life Science) and stored at 4°C.

The CM phase was either a 300 mM d-sorbitol (Merck) solution or DMEM. The CS phase was obtained by detaching cells from the walls of the culture flask with 0.5% trypsin-EDTA (Invitrogen) and by resuspending them in 300 mM sorbitol solution at an approximate concentration of 2 × 10^{6} cells/ml. Most experiments were performed using the following injection flow rates: (i) for spherical capsules: AS (30 ml hour^{−1}), CM (20 ml hour^{−1}), and CS (20 ml hour^{−1}), and (ii) for tubular capsules: AS (2 ml hour^{−1}), CM (1 ml hour^{−1}), and CS (1 ml hour^{−1}). Once formed, spherical and cylindrical capsules were placed inside an incubator (37°C, 5% CO_{2}, ~100% relative humidity) for about 24 hours and then transferred to an observation chamber placed on the stage of an optical microscope.

### Image analysis

Phase-contrast images of the spherical capsule were analyzed using a custom-made gradient-based edge detection algorithm implemented in MATLAB (MathWorks) [see (*11*) for details]. To determine the temporal evolution of the cylindroids length, we applied a variance filter followed by morphomath treatment and median filter to all images, and we derived the length *L*(*t*) by measuring the Feret diameter (*27*). This parameter is well adapted to measure sizes along one specific direction. It relies on the principle of the caliper. Registration of cylindroids was performed using homemade macros in Fiji (ImageJ). The mean diameter through time was measured on a registered cylindroid by edge detection and reslice of the central section in time. Cylindroid profiles for deformation curves were detected manually.

### Pressure derivation from spherical capsule dilation

Capsule deformation was monitored using video microscopy, and the pressure exerted by the MCTS was calculated using the formalism of thick-walled internally pressurized spherical vessels. We assumed that the alginate gel was isotropic and that strains were small (i.e., <10%). However, if the condition *t/R* < <1 is not fulfilled, the assumption of constant tangential stress across the thickness of the vessel does not hold. In the general case of Poisson’s ratio ν ≠ 1/2, we shall recall (*28*) the expressions for the radial and hoop stresses*R*_{in} *≤ R ≤ R*_{out}.

The radial displacement *u*(*R*) is obtained from Hooke’s law

Collecting these results, we arrive at

If the material is incompressible, we may simplify this equation and write it for the particular case of, for instance, *R = R*_{in}

Last, from volume conservation of the shell, we have

Using this equation, we separate the two time variables *R*_{in}(*t*) and *R*_{out}(*t*) and write the pressure *P*(*t*) as a function of *R*_{in}(*t*). Experimentally, we thus need to measure only the initial outer and inner radii and track the evolution either the inner radius of the capsule

### Pressure derivation from tubular capsule dilation

In the tubular capsule configuration, stress is recovered numerically to avoid restrictive hypothesis of the analytical model. The deformation of the capsule is monitored experimentally at several times along the tube main axis (see fig. S6). The difference between the initial tube radius and its temporal evolution permits to obtain the radial relative displacement of the alginate. For each recorded time, this displacement is imposed as BC in a finite element (FE) model of the cylindrical capsule implemented in Cast3M (FE code of the French Atomic Energy Commission, http://www-cast3m.cea.fr/). Once results are computed, the stress induced by tumor aggregate is obtained by extracting the normal stress at the inner bound of the cylindrical capsule FE model.

### Derivation of the mathematical model

*Mixture definition and conservation equations*. When the spherical/tubular capsule is formed, the CS (in the inner tube) and the culture medium (CM, in the intermediate interspace) come in contact and form a mixture. This mixture has a liquid nature, and its physical properties are correlated with the presence and relative abundance of the dissolved species. Two main dominant species can be identified in the mixture: the tumor cell species (*t*) and the medium fluid species (*m*). Once built, the spherical/cylindrical capsule is placed in a culture well where certain conditions (e.g., medium properness, nutrients) are prescribed to assure growth of cell aggregates. Between these conditions, oxygenation of the medium in the culture well has a capital role.

Oxygen is in the external culture medium but can also diffuse across alginate walls to reach the intracapsular space allowing for cells oxygenation. Hypoxia is one of most critical factors affecting tumor cell proliferation rate (*4*, *22*). Hence, oxygen is accounted as the nutrient species affecting tumor cell growth and metabolism and explicitly considered in the computational model as species of the mixture. Its mass fraction is indicated with ω* _{n}*. In addition to oxygen, other diluted chemical species are present in the mixture, but they are not considered explicitly.

If a properly representative elementary volume of mixture is defined, each point of the intracapsular space can be characterized by a certain mass fraction of *t*, *m*, and *n*: ω* _{t}*, ω

*, and ω*

_{m}*respectively (Fig. 2). Obviously, this constraint must be respected*

_{n}**v**, defined as the weighted sum of species velocities (

*14*)

The deviations of species with respect to mixture velocity are called diffusion velocities and read

Combining equations (Eqs. 10 to 12), the following condition holds

These definitions allow to properly describe the motion of cells and other species within the mixture. In particular, tumor cells may move due to advective transport (related to mixture velocity **v**) and diffusive transport (related to its own diffusive velocity, **u*** _{t}*). So the spatial form of the mass conservation equation of the

*t*species reads

*r*is an exchange of mass to account for chemical-biological events (e.g., cell division) inducing mass transfer from other species of the mixture to

_{t}*t*species. Two analogous equations govern mass conservation of the medium species,

*m*, and oxygen species,

*n*

*r*and

_{m}*r*are source or sink terms to account for mass exchange between species (e.g.,

_{n}*r*is the mass of oxygen consumed by cells due to metabolism and growth).

_{n}Mass created in one species (for example, in the *t* species due to cell division) has to come from other species, so the following constraint condition holds

Summing mass conservation equations of all species of the mixture (*t*, *m*, *n*, and those nonexplicitly considered) and accounting for constraint Eqs. 10, 13, and 17 give the mass conservation equation of the mixture as

The mixture behavior has also to obey to its momentum conservation equation, which reads**t** is the mixture stress tensor, and **b** are volumetric forces.

*Model hypotheses*. The following two hypotheses can be reasonably assumed in the computational model to simplify the set of governing equations:

1) The local density of the mixture, ρ, varies negligibly (so it is considered as a constant in the computational physics-based model).

2) Mixture velocity, **v**, is very low, so it is neglected in the equations.

Thanks to these hypotheses, Eq. 18 is satisfied. Furthermore, Eqs. 14, 16, and 19 can be rewritten as follows**t** = − *p***1**, with *p* being the mixture pressure. Equations 20 to 22 allow describing the behavior of the mixture. However, Eqs. 20 to 22 constitute a nonsolvable system of PDEs because the number of system unknowns (species mass fractions, diffusive velocities, etc.) is larger than that of equations, so the system has an infinite number of solutions. That is, the physical model is ill-posed, and additional hypothesis/relationships are needed to properly capture specificity of mixture behavior.

### Mixture potential and diffusion of tumor cell species

A diffuse interface approach was used to track interface between the cell aggregate and the medium. An interface’s thickness and the composition profile through the interface are set by the competition between individual motion and reciprocal attraction/repulsion of cells. The resulting free energy density, * _{t}* and its spatial derivatives (

*15*). Expanding about the homogenous state after symmetry considerations, the simplest model for the free energy density consists of two parts: a first gradient energy term, which accounts for the excess of free energy due to the inhomogeneous distribution of mass fraction in the interfacial region, and a bulk energy term, which depends on the assumed double-well potential (Eq. 2). We use here the form proposed in (

*20*)

The chemical potential, defined as the variation (rate of change) of free energy in the overall capsule volume * _{t}*, reads

The constant α in the previous equation is a normalization term, which depends on the assumed form of Eq. 2. It can be calculated from the 1D equilibrium analytical profile, imposing that the excess of the free energy density per unit surface area equals the fluid-fluid interfacial tension (*20*). According to the expression of the assumed bulk free energy density, calculated α reads

Diffusion of tumor cells is driven by the mixture chemical potential, μ, and its local distribution over the domain*M* is the tumor cell mobility. Active cell movement driven by biological and/or chemical signals is a nondominant mechanism in the analyzed cases, so it is not accounted in the mathematical model.

*Diffusion of oxygen*. As oxygen is a diluted species, its diffusion is assumed to follow Fick’s law* _{m}* in the considered case) allows for diffusion, a first approximation for the effective diffusion is provided from Rayleigh (

*29*)

^{−9}m

^{2}/s) (

*30*).

Oxygen diffuses also across capsule alginate walls. As alginate walls are not modeled explicitly to account for diffusion across them, a convective-type condition is assumed at the inner walls for oxygen normal inflow, *q**h*_{alg} is an equivalent convection coefficient, *t* is the thickness of the alginate wall, and ^{−9} m^{2}/s (*30*)]. In the two considered cases with *t* of the order of 30 μm, *h*_{alg} ≈ 4.5 × 10^{−5}. Oxygen mass fraction in the external culture medium is calculated with Henry’s law providing

*Intraphase exchange of mass terms: r _{t} and r_{n}*. Cell division is modeled as a reaction term (

*r*in Eq.

_{t}*20*) calibrated on replication rate measured experimentally. We hypothesize mitosis requires a certain amount of oxygen as well as medium species. Mechanical pressure within the mixture is also taken into account and assumed as inhibitor of cell division process (

*4*). Hence, depending on the pressure in the mixture and oxygen mass fraction, either cells have optimal rate replication or growth rate is partially or fully inhibited (fig. S2). To account for dependence on oxygen and pressure, the following equation is here proposed

*is a coefficient identified experimentally,*

_{t}*H*

_{ωn}and

*H*are two regularized Heaviside functions, which depend respectively on ω

_{p}*and*

_{n}*p*(mixture pressure) and vary between 0 and 1.

Each one of such regularized Heaviside functions depends on two parameters, *H* having the following form

The product of *H*_{ωn}(1 − *H _{p}*) is qualitatively depicted in fig. S2. The parameter

*H*= 0) and

_{p}*H*

_{ωn}= 1). A similar form has the function proposed for oxygen consumption

_{ng}allows accounting for consumption related to cell division, while the parameter γ

_{n0}serves to model oxygen consumption due to cell metabolism. Note that regularized Heaviside function and associated parameters are the same of Eq. 30. Hence, Eq. 32 can also be rewritten in the following equivalent form

*Interphase volumetric forces*. As indicated in the main text, chemical potential (*1*) times the gradient of tumor cells mass fraction allows to account for pressure difference between medium and cell aggregate induced by surface tension (*20*). Being the interface aggregate-medium diffuse, this is considered via a body forces in Eq. 22, which reads

*Final set of equations and computational aspects*. Introducing the defined constitutive relationships, a solvable final set of equations is obtained* _{t}*, chemical potential μ, nutrient mass fraction ω

*, and the mixture pressure*

_{n}*p*are the primary variable of the computational model.

The equations are discretized in space by the FE method and in time by finite differences and solved numerically with FEniCS (*31*). Table 1 reports material and geometrical parameters used in the numerical simulation.

To test the predictive potential of the mathematical model, material parameters have been firstly identified by inverse modeling of the tubular capsule case and then used for simulation of the spherical one. Uncertainties quantification (*32*, *33*) could be performed in a successive study to better assess uncertainties in the CCT protocol and their impact on the output of the mathematical model.

Convergence of the numerical solution has been analyzed to set the adequate FE size and time step (spatial and temporal discretization of model equations). Homogenous FE meshes have been generated and used for the numerical simulations (linear triangular elements having a mean characteristic size of the order of 0.5 μm). With fine meshing only needed in the interface zone, a local remeshing technique (*34*, *35*) would surely help in reducing computational costs as perspective improvement of the developed computational framework.

## SUPPLEMENTARY MATERIALS

Supplementary material for this article is available at http://advances.sciencemag.org/cgi/content/full/6/13/eaaz7130/DC1

Movie S1. Phase-contrast time-lapse video of a growing multicellular spheroid in a spherical capsule.

Movie S2. Phase-contrast time-lapse video of growing multicellular cylindroids in a cylindrical capsule.

Movie S3. Numerical simulations of aggregate growth in spherical capsule.

Movie S4. Numerical simulations of aggregate growth in tubular capsule.

Fig. S1. Images of tumor aggregate growth in alginate capsule.

Fig. S2. Pressure and oxygen mass fraction tumor cell growth dependence.

Fig. S3. Numerical diffusive flow in alginate capsule.

Fig. S4. Influence of hypoxia on volume evolution for spherical and tubular capsule.

Fig. S5. Numerical radial plots of tubular capsule primary variables.

Fig. S6. Tube radius as a function of distance from the cylindroid center along the main axis.

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is **not** for commercial advantage and provided the original work is properly cited.

## REFERENCES AND NOTES

**Acknowledgments:**P.N. and K.A. thank L. Munn, D. Drasdo, and J. Prost for the fruitful discussions about cylindroid growth.

**Funding:**This work was supported by the ANR VascTubes (ANR-15-CE18-0019) and, in part, by the ANR MecaTiss (ANR-17-CE30-0006-03).

**Author contributions:**V.L.M, P.N., and G.S. wrote the manuscript. B.G. and K.A. performed the experiments. B.G., K.A., and P.N. contributed to the experimental design and analysis. V.L., G.S., and H.B. contributed to the development of the mathematical model and performed numerical simulations.

**Competing interests:**The authors declare that they have no competing interest.

**Data and materials availability:**All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. Additional data related to this paper may be requested from the authors.

- Copyright © 2020 The Authors, some rights reserved; exclusive licensee American Association for the Advancement of Science. No claim to original U.S. Government Works. Distributed under a Creative Commons Attribution NonCommercial License 4.0 (CC BY-NC).