Poroelastic osmoregulation of living cell volume

Summary Cells maintain their volume through fine intracellular osmolarity regulation. Osmotic challenges drive fluid into or out of cells causing swelling or shrinkage, respectively. The dynamics of cell volume changes depending on the rheology of the cellular constituents and on how fast the fluid permeates through the membrane and cytoplasm. We investigated whether and how poroelasticity can describe volume dynamics in response to osmotic shocks. We exposed cells to osmotic perturbations and used defocusing epifluorescence microscopy on membrane-attached fluorescent nanospheres to track volume dynamics with high spatiotemporal resolution. We found that a poroelastic model that considers both geometrical and pressurization rates captures fluid-cytoskeleton interactions, which are rate-limiting factors in controlling volume changes at short timescales. Linking cellular responses to osmotic shocks and cell mechanics through poroelasticity can predict the cell state in health, disease, or in response to novel therapeutics.

Cell height changes can be finely captured by defocusing microscopy Water permeation and cellular deformability regulate dynamics of cell volume changes Poroelasticity describes the dynamics of cell volume changes The response of cell to hypo or hyperosmotic shocks are modeled by poroelasticity

INTRODUCTION
Cell volume regulation involves key homeostatic mechanisms that act in concert to maintain a well-defined intracellular ionic environment when external osmotic challenges arise (Chamberlin and Strange, 1989;Strange, 2004). For that, passive water and ion flow through membrane channels, together with active ion exchange via ATP-driven co-transporters, help to maintain cell volume in health and disease states (McManus et al., 1995;Cadart et al., 2019). Several mathematical models have been proposed to describe such volume adaptations (Herná ndez and Cristina, 1998;Lim et al., 2006;Solenov et al., 2011;Sachs and Sivaselvan, 2015;Cadart et al., 2019). The true complexity of volumetric responses displayed by cells exposed to external osmotic perturbations is not explained by modeling cells as simple water-filled receptacles composed of a lipid-rich plasma membrane (Fels et al., 2009). Instead, the mechanics of the cytoskeleton-a key determinant of cellular behavior, phenotype, and shape (Cosgrove, 1993;Ofek et al., 2009;Jiang and Sun, 2013)-is intimately linked with the volume of mammalian cells (Ferná ndez and Pullarkat, 2010). Additionally, cell shape and volume are regulated by geometrical constraints, adhesive and mechanical interactions with extracellular environment (Ferná ndez and Pullarkat, 2010;Ambriz et al., 2018), the rheological properties of cytoskeleton (Ferná ndez and Pullarkat, 2010), the spatial distribution of water flux across the cell membrane, and hydraulic pressure acting at the cell surface (Moeendarbary et al., 2013;Han et al., 2020;Tao et al., 2017). In the context of tumor biology, for example, it has been shown that cancer cells positioned on the outer layers of tumor spheroids display larger volumes and are softer, conferring a more invasive phenotype compared to cells located within the tumor core (Han et al., 2020). Analogous to a fluid-filled sponge, a solid skeleton of the cell, consisting of cytoskeletal proteins, organelles, and macromolecules, is permeated by intracellular fluid. The dynamics of these fluid movements within the gelatinous cytoplasm can control the rate of adaptive cell volume changes under both hypoosmotic and hyperosmotic conditions (Strange, 2004;Moeendarbary et al., 2013).
Under the poroelastic framework, cells are treated as fluid-filled spongy materials, providing a more representative depiction of their inherent hydrogel-like composition (Fels et al., 2009;Moeendarbary et al., 2013;Gautier et al., 2015;Moeendarbary and Charras, 2015). Hydrogels are cross-linked polymer networks permeated by water molecules and their response to mechanical forces depends on the hydrodynamic interaction between the solid phase (e.g. in cells, the cytoskeleton) and the fluid phase (the cytosol). Indeed, the velocity of fluid redistribution through the pores of the solid meshwork and the deformability of the solid meshwork control the rate at which hydrogels can expand or shrink (Nam et al., 2004;Mitchison et al., 2008;Persson and Halle, 2008;Moeendarbary et al., 2013;Esteki et al., 2020). In addition to the hydrogel-like passive mechanical behavior, cells actively remodel their cytoskeleton at minute timescales, potentially influencing the dynamics of shape change (Li et al., 2018;Garenne et al., 2020;Rizzelli et al., 2020). On the one hand, fluid flows into a cell due to hypoosmotic conditions lead to cell swelling and increase the intracellular hydraulic pressure (Sardini et al., 2003;Wilson and Mongin, 2018). Consequently, the cell swelling may lead to changes in gene transcription via Ca 2+ regulated mechanotransduction signaling pathways and relocalization of transcription factors (Okada et al., 2001;Nilius et al., 2003) as the cytoskeleton is connected to protein complexes that regulate the opening and closing of mechanosensitive ion channels (e.g. TRPV4). Swelling can also cause membrane blebbing and increase the metastatic potential of cancer cells (Khan et al., 2019). On the other hand, hyperosmotic stress and cellular shrinkage may decrease intracellular hydraulic pressure and lead to pathologies induced by either cytoplasmic stiffening (Prystopiuk et al., 2018;Han et al., 2020), DNA damage and apoptosis (Burg et al., 2007), oxidative stress and inflammation (Vaziri and Rodríguez-Iturbe, 2006), or cell-cycle arrest, and cellular senescence (Myers et al., 2007). Indeed, as cells age, the density and cross-linking of the cytoskeletal network (i.e. F-actin filaments) increases leading to a decrease in cytoskeletal pore size, slower fluid diffusion within old cells (Phillip et al., 2015), and cellular stiffening (Wilson and Mongin, 2018). As such, a poroelastic theory that mechanistically captures cell volume changes in response to osmotic changes can be used to make predictions regarding the biological significance of cell volume changes in different cell types (Mitchison et al., 2008;Moeendarbary et al., 2013;Sachs and Sivaselvan, 2015).
Pioneered by Terzaghi for the 1D settling of water-saturated soils (Terzaghi, 1925), and later generalized and extended in 3D by Biot (Biot, 1941(Biot, , 1956, poroelasticity has been applied to study rheological behavior of hydrogels and more recently fluid flows and pressurization in cells (Hu et al., 2010;Chan et al., 2012;Moeendarbary et al., 2013;Esteki et al., 2020). Passive intracellular pressurization often arises from solute concentration differences between the extracellular and intracellular spaces, which are separated by a semipermeable cell membrane. In the absence of active ion transporters, the cell membrane allows solutions to equilibrate, that is, water permeates toward the solution of higher ionic concentration causing an extra pressure on the membrane defined as the osmotic (or turgor) pressure (P) ( Van't Hoff, 1966;Liu et al., 2019aLiu et al., , 2019b. Application of hydrostatic pressure (p) difference across the cell membrane can also drive water movement. Therefore, the net pressure inside the cell (P in = p i -P i ) and the net pressure outside the cell (P out = p o -P o ) define the total effective pressure (P eff. = P out -P in ) across the cell membrane which drives the flow of water in or out of the cell. A positive effective pressure difference (P eff. >0) causes fluid to flow into the cell and swelling. Conversely, a negative effective pressure difference (P eff. <0) leads to shrinkage. The swelling or shrinkage involves the water flowing into or out of the cell by crossing the lipid cell membrane. While swelling or shrinkage depends on the water permeation that primarily relies on the membrane permeability as well as the hydraulic resistance of the cytoplasm, the rheological properties of the cytoplasm critically modulate the rate of the volume changes (Moeendarbary and Charras, 2015;Liu et al., 2019aLiu et al., , 2019bMalandrino and Moeendarbary, 2019). Poroelasticity, as a coarse-grained bottom-up framework (Bausch and Kroy, 2006), can account for these factors and thus provides a mechanistic framework for understanding the swelling/ shrinkage kinetics.
In recent swelling experiments, the mechanical behavior of several hydrogels and tissues has been characterized (Chandran and Barocas, 2004;Yoon et al., 2010;Bouklas and Huang, 2012;Gupta and Shivakumar, 2012;Bouklas et al., 2015;Wong et al., 2015;Castro et al., 2018;Tanska et al., 2020) through determining the poroelastic parameters by analytical solutions and investigating their effects on the transient and steady-state responses (Buckley et al., 2009;Yoon et al., 2010;Bouklas and Huang, 2012;Gupta and Shivakumar, 2012;Li et al., 2013;Bouklas et al., 2015;Caccavo et al., 2018;Yu et al., 2018). In these studies, it is assumed that initial pressurization ramps instantaneously (i.e. rise time t r = 0), while an instantaneous pressure rise assumption is mostly unrealistic under experimental circumstances and may lead to incorrect estimation of the swelling/deswelling poroelastic parameters during a transient osmotic stimulus (Javanmardi et al., 2021). Moreover, active processes at the cell membrane (i.e. cellular regulatory volume processes via ion channels) and osmotic and hydrostatic challenges that cells experience in different physiological and disease states may create gradual effective pressure differences, that is, with a finite rise time. It has been concluded that neglecting the effects of finite rise times (t r ) during the analysis of the indentation of poroelastic materials may lead to unrealistic estimations of poroelastic parameters (Oyen, 2008;Galli and Oyen, 2009;Hu and Suo, 2012;Mollaeian et al., 2017;Mollaeian et al., 2018aMollaeian et al., , 2018b  iScience Article 2019a, 2019b; Esteki et al., 2020). However, the influence of finite pressure rise-time associated with osmotic and hydrostatic challenges on the realistic evaluation of swelling/deswelling dynamics has not been considered yet.
Here, we employed a combination of cell culture experiments, high-resolution nanoparticle tracking, analytical solutions, and finite element (FE) simulations to investigate the poroelastic swelling/deswelling behavior of adherent cells in response to transient pressure differences caused by osmotic challenges. We measured the dynamics of HeLa cell swelling and deswelling by defocusing microscopy and ran a series of analytical calculations and FE simulations to identify critical poroelastic parameters, including geometrical and rise times effects. In our modeling framework, considering hydraulic (p o À p i ) and osmotic (P o À P i ) pressures, we defined the effective pressure difference as P eff: = ðp o À p i Þ-ðP o ÀP i Þ , which is a net driving pressure for cellular volumetric deformations. This phenomenological approach has the advantage of being directly implementable in the FE simulations. Considering a realistic 3D cell geometry, we investigated the poroelastic responses of cells under finite pressure rise-time. Our results indicate that dynamics of cell volume changes at short timescales can accurately be described by the poroelasticity framework and our estimated cellular poroelastic parameters are consistent with previously reported values.

Experimental results
The temporal evolution of the volume of HeLa cells in response to osmotic shocks was measured by acquiring fast xyzt spinning disk confocal stacks. Figure 1A shows an example of orthogonal profiles of a HeLa cell before and after application of a hyperosmotic shock. The decrease in cell volume resulted in an increase in fluorescence intensity indicating the increase in GFP concentration due to water expelling. The 3D reconstruction of cell profile from z-stacks before and after application of the shock displayed a substantial cell height decrease ( Figure 1A, bottom projections) suggesting that probing the cell height is a good indicator of the volume changes. We then measured cell volume changes in response to changes in osmolarity for up to 30 min to ascertain that volume changes persisted long enough for our AFM indentation measurements ( Figure 1B). We verified that photobleaching due to the acquisition of confocal stacks over an extended period did not artefactually affect our measurement of cell volume (control, Figure 1B). To ensure a stable volume increase in hypoosmotic conditions over the timescale necessary for AFM experiments ($30 min), during measurements, we added regulatory volume decrease inhibitors to the medium (Lang et al., 1998). We measured an unstable cell volume increases up to 10% upon addition of water for hypoosmotic treatment while a stable cell volume increase of 22 G 2% was achieved upon addition of water + N + D. Conversely, upon addition of 100 and 250 mM sucrose, cell volume decreased to 12 G 3% and 21 G 6% of the original volume, respectively.
We maximized the rate of application of the osmotic shocks by exchanging the media of different osmolality immediately (over less than a few seconds). Our results indicate that most volume changes occur within the first minute following the exchange of the media ( Figure 1B). Furthermore, employing spinning disk confocal z-stacks and image processing algorithms, we estimated errors of at least 20% associated with capturing cell height (i.e. 400 nm is the best resolution we could achieve in z direction using confocal microscopy). Also, the maximum rate at which we could acquire the stacks covering a whole cell was $1 min. Taken together, despite spinning disk confocal imaging providing a good estimation of volume changes over long timescales (tens of minutes), it is not possible to capture the full transitory dynamics of cell height changes (occurring at millisecond to second timescales) with this technique.
To address the limiting factors (i.e. lack of spatial and temporal resolutions) associated with typical confocal 3D imaging and capture the short timescale dynamics (shaded area in Figure 1B), we employed epifluorescence defocusing microscopy (see STAR Methods). We exposed HeLa cells to hyperosmotic or hypoosmotic medium at t r $0 while tracking the vertical displacement of fluorescent nanobeads as a function of time and analyzing the defocused 2D images taken at 100 ms intervals ( Figures 1D and 2A). The changes in the thickness of the cells d at the bead positions (located at specific heights) as a function of time are plotted in Figures 2B and 2C for swelling and deswelling conditions. In swelling experiments, most vertical displacements were larger for greater cell heights while this trend was weaker in deswelling cases. Furthermore, the timescale associated with deswelling was smaller and most curves reached a stable plateau after 20s compared to 50s for swelling curves that reach a less stable plateau. These observations suggest the inherent differences in the dynamics of cell swelling and deswelling that mainly originate from different ll

OPEN ACCESS
iScience 24, 103482, December 17, 2021 3 iScience Article mechanical responses of the cell membrane, cytoskeleton, and other cellular structures (nucleus, organelles, macromolecules) to stretch or compression (Isenberg et al., 2003;Chao et al., 2018). When normalizing the vertical displacements d ( Figures 2B and 2C) to their respective maximum final vertical displacement d N ( Figures S3A and S3B), we noticed distinct dynamics with varied responses of different positions Inset shows the averaged radial intensity curves normalized to the maximum intensity in their center. (F) Corresponding calibration curve demonstrating the linear relationship between the radial distance of outer ring and the height (z-position). The radial position of the ring (peak intensity) is linearly related to the distance between the imaging plane and the nanobead focal plane and therefore the z-position of the nanobeads can be tracked with very high accuracy ($10 nm). iScience Article on the cell surface and also at short and long timescales. Interestingly, the normalized displacement curves did not collapse on each other and therefore, we asked whether this observation is, in principle, valid when considering the cell to be an ideal poroelastic material.
We also investigated the effects of hypoosmotic or hyperosmotic shocks on the cytostructural organization of the HeLa cells and particularly F-actin which is the main determinant of the cell mechanical properties. The Factin cytoskeleton did not reorganize in response to hyperosmotic treatment ( Figure 2D) or hypoosmotic shock ( Figure 2E), as the density and morphology of actin filaments were not noticeably altered. These observations, together with the contribution of intrinsic viscoelasticity of filaments (arising from cytoskeletal reorganization) to relaxation times at long timescales (over 60s), suggest that considering a purely elastic solid skeleton at short timescales (less than 60s) is a well-justified consideration in our poroelastic model.
To measure the elastic modulus required in the estimation of poroelastic parameters in osmotic experiments, we conducted AFM indentation tests (see STAR Methods). Additionally, these indentation tests provided an independent way of measuring poroelastic diffusion constant (Hu et al., 2011;Chan et al., 2012;Esteki et al., 2020) and showed how changes in cell volume influence the poroelastic parameters.
In the case of swelling, addition of water increases cell volume and significantly increases (55%) the poroelastic diffusion constant and decreases (23%) cellular elastic modulus ( Figures 2F and 2G).

Poroelastic analysis of simple 2D geometries
We first considered a rectangular geometry and linear poroelastic FE model ( Figure 3A and Table S1) to investigate swelling/deswelling in response to osmotic perturbations that were modeled by imposing an effective pressure P eff. boundary condition that linearly reaches its final value over the rise time t r ( Figure 3B).

OPEN ACCESS
iScience 24, 103482, December 17, 2021 5 iScience Article Figure 3. Effects of spatial location, rise time, effective pressure for constrained swelling of a simple geometry investigated by FE linear poroelastic simulations (A) Idealised swelling of a thin layer of a poroelastic material attached from one side to a substrate and initially at equilibrium P eff : .$ P out -P in = 0. Axisymmetric FE model of a thin disk with thickness H and length L was built to simulate effective pressure-driven swelling considering P eff : > 0 / P in < P out . Vertical displacements of specific nodes located at different positions P1, P2, and P3 were studied. (B) Linear increase of effective pressure from 0 to P eff : = 0.05 kPa over a rise time t r = 0.1s. (C) The swelling type deformation of the hydrogel due to application of pressure difference and inflowing of water. Dash and solid lines indicate nondeformed and deformed material boundaries respectively. The fluid velocity field is distributed non-uniformly at boundaries and within the domain and induces a non-uniform swelling behavior across the geometry.
(D) Plot of vertical d over time considering the ramp in pressure indicated in (B).
(E) Normalization of (D) considering d dN for y axis and plots of 1D finite thickness and self-similar solutions. (F) Normalization of (D) considering t = 4Dt pH 2 for x axis and d dN for y axis. (G) Application of ramp with two different P eff : and two different t r .

OPEN ACCESS
The rise time can be interpreted as the time of pressurization during deswelling, or depressurization during swelling. The consideration of rise time in our simulation is a new attempt to phenomenologically capture the effects of active membrane processes that make the cell membrane a time-dependent barrier to the transport of molecules and water in or out of the cell (Yamamoto et al., 2014;Yang and Hinner, 2015;Li et al., 2016;Emami et al., 2018;Yang et al., 2019). Furthermore, cells are not exposed to stepwise, instantaneous osmotic perturbations in vivo (Haskew-Layton et al., 2008;Pilizota and Shaevitz, 2012), and in our in vitro experiments, it takes at least a few seconds before cells are fully exposed to the media with different osmolarities despite our efforts to exchange the media immediately.
We tested the spatial heterogeneity of swelling/deswelling dynamics as the experimental vertical displacement curves differed largely for different positions ( Figures 2B and 2C). We considered a set of poroelastic properties in the range of our AFM cellular measurements (E = 1kPa, D = 10mm 2 s À1 ) and effective pressure of 50 Pa which induces deformations at the scale of the experimental osmotic perturbations. Figure 3C shows swelling deformations and local fluid velocities predicted from the linear poroelastic FE model considering a short rise time of t r = 0.1s ( Figure 3B) and the constrained (no-slip condition for the bottom surface) boundary condition (see Figure S4A for constrained deswelling and Figure S6A for free swelling). Different local fluid velocities (Figures 3C and S4B) were observed at different positions of the rectangular geometry (P1, P2, and P3) leading to spatially different vertical displacements curves (d vs t) for both swelling/deswelling conditions ( Figures 3D and S4C). Also switching from no-slip to sliding boundary condition led to significantly different vertical displacement dynamics ( Figure S6). To better compare constrained and free swelling, the velocity field and deformations for different boundary conditions are shown in Figure S8. Analysis of deformations and fluid profiles, suggests that a positive pressure gradient causes fluid entrance and swelling deformation ( Figures S8A and S8C), while the negative pressure gradient induces expulsion of fluid and deswelling/shrinkage ( Figures S8B and S8D).
Similar to the experimental observations ( Figures S3A and S3B), normalization of vertical displacements with respect to their respective maximum final vertical displacement, d N , did not lead to collapse of curves ( Figures 3E and S4D) particularly at short time scales (t < 1 s). This demonstrates that different points at the cell surface can show unique poroelastic dynamics depending on their location relative to the boundaries. In Figure 3E, we also plotted the 1D finite thickness (Equation 25) and self-similar (Equation 26) analytical solutions considering stepwise ramp (considering the same poroelastic properties as in our FE simulations and H = 5.5 mm). The self-similar solution only captures the dynamics of deformations at very short time scales (t<<0.2s), while the finite thickness solution roughly follows the dynamics for the P1 position, which is the point that experiences the most symmetrical condition and the least boundary effects. Furthermore, normalization of time by t = 4Dt pH 2 based on Equation 26 and previous reports (Yoon et al., 2010; Bouklas and Huang, 2012)) did not improve overlaying of the curves ( Figure 3F). We next studied the sensitivity of the model predictions to t r and P eff considering three different scenarios: (i) P eff. = 0.03 kPa and t r = 0.6s, (ii) P eff. = 0.05 kPa and t r = 1s and (iii) P eff. = 0.05 kPa and t r = 0.6s ( Figure 3G). Figures 3H and S4G shows the d-t curves (at position P1) for swelling and deswelling cases and suggest that for a fixed maximum pressure (P eff. = 0.05 kPa) but different rise times (t r = 1s or 0.6s), the shape of d-t curves is different particularly at short-time scales. This is mainly due to different fluid redistribution in the early phases of pressure relaxation. At longer time scales, the relaxation curves tend to overlay each other and plateau to a stable value of d N at time t N . In swelling, a stable value of d N $0.4 mm was reached in t N $4 s while for deswelling, it took t N $3.5 s for d N $0.4 mm ( Figures 3H and S4G). Interestingly, further normalization of time using 4Dt pH 2 illustrates that curves with equal values for both P eff : dN and t r fully overlay onto each other ( Figures 3I and S4H), indicating the unique dependence of poroelastic dynamics to both time and length scales. Similar results were obtained by employing a porohyperelastic model (Figure S5), confirming the suitability of porohyperelasticity for application to more complex settings such as those we investigated in the next section. Moreover, we additionally examined the effects of non-linear increase in the P eff. (t r ) and confirmed that the height dynamics did not change ( Figure S9) compared to the simple linear ramp function, which is assumed throughout the rest of the simulations.  Figures S4D, S4E, S4H, and S4I) highlight the noncompliance of analytical solutions with FE simulations, primarily due to different geometrical and boundary conditions. In fact, both analytical solutions (i.e. self-similar and finite thickness 1D models) assume simple 1D boundary with parameters that depend only on one direction while in our FE model flow and pressure distributions are multidirectional. Moreover, consideration of instantaneous application of pressure in analytical solutions versus finite rise time in FE simulations is another major reason for different solutions. The timescale of pressurization has been ignored in several recent studies (Rice and Cleary, 1976;Yoon et al., 2010;Bouklas and Huang, 2012;Li et al., 2013;Bouklas et al., 2015;Tanska et al., 2020) while here, it appears to be of fundamental relevance. Indeed, during the deformation of a poroelastic material, the key mechanism for fluid pressure changes relies on the ability of the interstitial fluid to redistribute through the pores of the solid phase. The distinctive feature of the poroelastic mechanical behavior is that the length scale l and the timescale t of pressurization or relaxation are correlated through t p $ l 2 =D, where D is the poroelastic diffusion coefficient (Hu et al., 2010;Yoon et al., 2010;Bouklas and Huang, 2012;Kalcioglu et al., 2012;Esteki et al., 2020). The presence of a characteristic length scale l is thus a unique feature of poroelastic mechanical responses. In our study, the length scale for both the pressurization and the subsequent fluid relaxation is the material thickness l $ H. Here, we propose a limiting case when this length scale and the poroelastic diffusion constant are such that t r and t p $ H 2 =D are of the same orders, that is, when pressurization and volume changes occur on similar timescales. On the contrary, our model reconciles with previous analyses (Hu et al., 2010;Yoon et al., 2010;Bouklas and Huang, 2012;Kalcioglu et al., 2012;Esteki et al., 2020) in when t r is significantly shorter than t p , t r ( t p $ H 2 =D. To examine the effects of different permeabilities of membrane and cytoplasm, we simulated the membrane by considering an additional top thin layer (500nm thickness) with different diffusion coefficient compared to the base layer ( Figure S10A). Indeed, the overall fluid redistribution can be captured by considering an apparent diffusion coefficient D that relies on both the diffusion coefficient of the membrane (D 2 ) and the diffusion coefficient of the intracellular space (D 1 ). No-slip condition between the two layers was assumed and other conditions were set similar to the single-layer FE model. We considered 2 sets of poroelastic properties (D 1 = 10mm 2 s À1 , D 2 = 10, 1, 0.1mm 2 s À1 and D 1 = 1mm 2 s À1 , D 2 = 1, 0.1, 0.01mm 2 s À1 ) and an effective pressure increase of 50 Pa. Figures S12C and S12E show vertical displacements curves predicted by the linear poroelastic FE model considering a short rise time of t r = 0.1s at position P1. Varying membrane diffusion coefficient led to different vertical displacement curves (d vs t). Interestingly, we found a correspondence between the use of different t r in the one-layer FE model (t r = 0.1, 0.6, 6s in Figures 3D and 3H) and different membrane diffusion coefficients ( Figures S10C, and S10E). These results suggest that considering different diffusion coefficients for the cell membrane changes the shape of d-t curves (Figures S10C and S10E) in a similar fashion as when changing rise time in the single-layer FE model. Therefore, considering a finite rise time in a single layer model may provide a simple way of mimicking the membrane effects.

Toward more realistic analysis for cellular setting
Our osmotic perturbation experiments indicate that cells undergo large deformations (strains over 50%) that may fall outside of the linear elasticity theory implemented in FE simulations. Also, our porohyperelastic simulations for a simple geometry ( Figure S4) suggest trends that were consistent with the linear poroelastic results (Figure 3). Therefore, we considered a porohyperelastic formulation (Table S2) for an axisymmetric geometry of a typical cell adhered to a substrate and refined the model by employing an unstructured mesh ( Figure 4A). Osmotic perturbations were simulated by imposing a linear ramp (P eff. = 50 Pa) over a fast rise time t r = 0.01 s ( Figure 4B), leading to water influx and a complex velocity field as illustrated in Figure 4C (and Figure S11B for deswelling). The vertical displacement vs time curves (d-t) in specific cellular locations (P1, P2, and P3 with respective heights H = 5.5, 4, 2.5 mm) are shown in Figure 4D. Normalization of vertical displacement using d dN did not show any meaningful overlapping of the curves. However, normalization of both vertical displacement and time (using t = 4Dt pH 2 ) resulted in some fair overlapping of curves for positions P1 and P2 ( Figure 4F). Furthermore, while self-similar (Equation 26) and 1D finite thickness (Equation 25) solutions did not fit the overall trend of FE curves, the self-similar solution shows a good fit for curves at short time scales ( Figure S7). Interestingly, for positions with larger heights (e.g. P1 and P2) the self-similar solution appropriately fit the d dN FE curves up to 1s (Figures S7A and S7B) while for P3 (that has a significantly lower height), the self-similar solution fit the FE curve only up to 0.5s ( Figure S7C). Consistent with our observations in the rectangular geometry (Figure 3), the differences in the dynamics of different locations is due to geometrical boundary effects: in particular, P1 lays on the symmetry axis, and relatively far from the bottom boundaries, while P3 is the farthest of the three points and closest to the bottom.  Figure 4H. Normalization of vertical displacement using d dN led to the overlap of the curves only when t r <<1s and 1D analytical solutions could fit these curves at short time scales. However, for t r >>1s (i.e. t r = 3 and 10s), the shape of the curves significantly differs, yet reaching a plateau of d N $1.65 mm at t N $11s ( Figures 4H and 4I), suggesting that at the steady-state the poroelastic relaxation responses converge to similar values. To consider the effect of transiently changing osmotic pressure in realistic geometries, and obtain the corresponding values of poroelastic parameters, we modeled the living cell as an isotropic porohyperelastic continuum (as described in the previous section), and extracted the material parameters by fitting experimental curves in Figures 2B and 2C with FE simulations. Considering the swelling/deswelling curves, fundamentally we can only extract three parameters (P eff. , D and t r ). Therefore, we used the experimental values of elastic modulus estimated from AFM indentation tests (average of swelled and control, Figure 2G, E = 0.975kPa, for swelling and average of shrunk and control, Figure 2G, E = 1.15kPa, for deswelling) and ran optimization procedures to fit the swelling/deswelling experimental curves ( Figures 2B and 2C) with FE simulations and estimate effective pressure (P eff. ), rise time (t r ) and poroelastic diffusion constant (D) for the swelling and deswelling of the HeLa cells. For simplicity, we considered a fixed Poisson's ratio of y = 0.3 as suggested in recent studies (Moeendarbary et al., 2013;Wu et al., 2018;Mokbel et al., 2020). To avoid the influence of active remodeling and intrinsic viscoelasticity of the cytoskeleton, we only considered the first 40 s of swelling/deswelling curves in fitting procedures. The optimization procedure involved initial running of the simulations to fit the experimental curves and find the first approximated values for Peff. and D, considering t r = 0. Then a step increase of 0.1s and small variations in P eff. and D were allowed to find the best parameters that fit the experimental curves with a minimal error.
Representative results of fitting for two swelling cases of two different cells (with nanobeads positioned at heights of 2.9 and 3.9 mm) and two deswelling cells (with nanobeads positioned at heights of 2.5 and 4.6 mm) are shown in Figures 5A and 5B. We found the following optimum effective pressures, diffusion constants and rise times for swelling: P eff. = 0.6 kPa, D = 2 mm 2 s À1 and t r = 4.5 s for H = 2.9 mm, and P eff. = 0.8 kPa, D = 2.4 mm 2 s À1 and t r = 4 s for H = 3.9 mm ( Figure 5A). For deswelling, we kept the same geometry but forced a negative effective pressure and found the following optimum parameters: P eff. = -0.7 kPa, D = 0.8 mm 2 s À1 and t r = 5 s for H = 2.5 mm, and P eff. = -0.8 kPa, D = 3 mm 2 .s À1 and t r = 2.5 s for H = 4.6 mm ( Figure 5B). Finally, running the fitting for all experimental curves, we estimated the average P eff. = 0.93 G 0.3 kPa, D ave . = 3.35 G 1.36 mm 2 .s À1 and t r = 3 G 1.47 s for swelling and P eff. = -0.72 G 0.21 kPa, D ave . = 1.88 G 0.92 mm 2 s À1 , t r = 3.92 G 0.92 s for deswelling ( Figures 5C, 5D, and 5E). iScience Article In our simulations, we only considered the influence of the cell geometrical features and therefore the abovementioned coefficients can only be regarded as approximations to the bulk, average poroelastic properties of a living cell. The structural heterogeneity and complex nature of the cell can be the key reason that different experimental settings report different values for cellular mechanical and rheological properties. Indeed, cellular structures such as the cell membrane, the cortex, the nucleus, and the granular filamentous interior, with different permeability and elasticity characteristics, contribute to setting the overall swelling/shrinking cell dynamics. Figures 2F and 2G show that cells exposed to sucrose and after deswelling exhibited no significant changes in elastic modulus or diffusion constant. Nonetheless, the mild increase in cell elastic modulus ($10%), although not statistically significant, is comparable in magnitude to the stiffening of cytoplasmic material after water efflux (Guo et al., 2017) and indicates that cells relax less rapidly and become stiffer with decreasing fluid fraction. Consistent with previous AFM indentation results by us and others, we measured the diffusion coefficients in the range from D = 10 to 20 mm 2 .s À1 (Charras et al., 2005(Charras et al., , 2009Moeendarbary et al., 2013) and elastic modulus of E $ 1 kPa (Lim et al., 2006;Moeendarbary et al., 2013;Gautier et al., 2015). Also, increases in cell volume resulted in a significant increase in poroelastic diffusion constant and a significant decrease in cellular elastic modulus. Our indentation simulations on the final swelled/shrunk state of a simple geometry (see STAR Methods) showed that D increased while E decreased for the swelled state, and reversely D decreased and E increased for the shrunk state (Figures S12 and S13) which are consistent with our experimental observations.

AFM indentation test results in
Within the timescale of our experiments, and separating swelling and shrinkage phenomena, a unique set of poroelastic parameters reproduced most of the features of fast volume changes in HeLa cells. This is important as it provides the ability to better understand the early adaptive physiological responses of living cells to external osmotic challenges. However, the diffusion coefficients estimated from experimental swelling/deswelling curves are around one order of magnitude smaller than those obtained from AFM indentation tests. Several reasons could contribute to this discrepancy. Firstly, in AFM tests, a thin layer of actin cortex located beneath the plasma membrane mostly contributes to the measured poroelastic parameters. During the osmotic perturbation, however, water infiltration through the whole cell, including the cytoplasm and specific regions of the cell such as cell nucleus, regulate the dynamics of swelling/deswelling (Potma et al., 2001). Furthermore, while in AFM tests the pressurized fluid redistributes only horizontally in small sections of the membrane because of the compression of the impermeable spherical indenter, the whole-cell membrane significantly controls the osmoregulated volume changes in swelling/deswelling experiments. Finally, the timescale of ramp (pressurization) and relaxation phases are different across the two methodologies. In the AFM experiments, we used t r $500ms while, in osmotic tests, it takes a certain amount of time for the osmotic perturbations to act on cell surfaces (t r in the order of a few seconds) and delays associated with media mixing could be another reason for the discrepancy in the estimation of diffusion coefficients.
Comparing deswelling and swelling experiments, the consistently faster poroelastic diffusion in deswelling might be due to the irreversible dynamics of membrane water channels (aquaporins) (Mitchison et al., 2008;Ibata et al., 2011) and to the variable responses of the cytoskeleton (and indeed whole-cell components) and to tensional and compressional forces created during hyperosmotic shrinkage and hypoosmotic swelling (Sehy et al., 2002;Mitchison et al., 2008). The imaging of actin cytoskeleton (which is the main determinant of cellular mechanical properties) before and after short-time of application of osmotic perturbations ( Figures 2D and 2E), suggests that while cytoskeleton can be considered elastic at very short timescales (as remodeling and reorganization of actin cytoskeleton occurs at long time scales), slight variations of cytoplasmic pore size and crowding can effectively influence the poroelastic properties at short timescales. Although the intrinsic viscoelastic timescale of the cell cytoskeleton may fall within a range close to poroelastic timescales, the early osmoregulation of the cell volume is predominantly regulated by fluid-solid interactions, that is, the poroelasticity.
Accurate predictions of cell volume changes in response to pressure challenges can help to quantify the mechanical behavior of living cells and contributes to the understanding and diagnosis of various pathologies (Liu et al., 2019a(Liu et al., , 2019b iScience Article organ-on-chip technologies (Borsali and Pecora, 2008;Feksa et al., 2018;Whelan et al., 2021). In the present work, we used the finite element method to simulate the osmotic-driven swelling and deswelling kinetics of the cytoplasm of a living cell. In cell biology, it is important to be able to model and predict the response of different cell types to osmotic perturbations using models that incorporate realistic time frames for physiological phenomena. In the context of tumor biology, for example, being able to predict the response of cancer cells to chemotherapies that alter their internal hydrostatic pressure may help to explain their effectiveness at treating tumors, or conversely, may reveal how tumor cells develop resistance to a particular drug (Micalet et al., 2021). Our model closely recapitulates the volumetric responses of HeLa cells with rapid redistribution of fluid due to hypo and hyperosmotic stimuli. We confirm that the poroelastic model, with minimal parameters, explains the swelling/deswelling of living cells in vitro and found that, in both simplified and cell-like geometries, the effect of duration of initial stimulation, that is, rise time, is critical. We also found that including more realistic geometries obtained from confocal reconstructions images of the cell improves uncertainties arising from spatial inhomogeneities and is simply sufficient for predicting the cell mechanical responses and volume regulation without the need for additional complex parameters.

Limitations of study
The main limitation of our work is the consideration of isotropy and constant values for poroelastic parameters. However, both simulations (Figures S12 andS13) and experiments ( Figures 2F and 2G) on swelled and shrunken cell states suggest that the effective elastic modulus and diffusion constants progressively change during volume changes. Interestingly, these changes are in the opposite direction for swelling compared to deswelling and this may explain the experimentally observed differences between the dynamics of swelling and deswelling. Also, additional features, such as the viscoelasticity of the solid-phase, were not considered in our model. In swelling experiments, most vertical displacements were larger for greater heights while this trend was weaker in deswelling cases. Therefore, we speculate that in the deswelling experiments, the height dynamics may strongly be affected by the position of the nucleus and the presence of the (stiffer) nucleus may have a stronger impact during deswelling compared to swelling. Notably, in this study we did not consider the presence of the nucleus and, therefore, the inclusion of the nucleus and its properties would be an interesting addition in future studies.

STAR+METHODS
Detailed methods are provided in the online version of this paper and include the following:

Materials availability
This study did not generate new reagents.

Data and code availability
Original/source data in this study is available upon request. All original code is available in this papers supplemental information. Any additional information required to reanalyze the data reported in this study is available from the lead contact upon request.

Cell volume measurements
Confocal stacks of cells expressing cytoplasmic GFP were acquired using a spinning disk confocal microscope (Yokogawa CSU-22, Japan, 100x oil immersion Olympus objective lens, NA = 1.4, piezo-electric zdrive, NanoscanZ, Prior Scientific, Rockland) to measure changes in cell volume in response to osmotic shock. Exposure time and laser intensity were set to minimize photobleaching and stacks (40 images separated by 0.2 mm) were acquired for a total of 30 mins with the shortest achievable time interval ($2 min).
Using Matlab (Mathworks Inc) Image Processing Toolbox, the background noise of stack images was removed, and images were smoothed and binarised. A series of erosion and dilatation operations were performed on binarised stacks to create a contiguous cell volume image devoid of isolated pixels (Figure 1A). The cell volume at each time step was calculated from the sum of nonzero pixels in each stack multiplied by the volume of a voxel. Prior to application of osmotic shock, four stacks were captured and then cell volume was followed for a further 25 minutes ( Figure 1B) after the change in osmolarity. For hyperosmotic shocks, cells were exposed to medium to which 100mM and 250mM sucrose had been added. For hypoosmotic conditions, we diluted culture medium using water and observed a transient swelling, signifying that regulatory volume decrease mechanisms acted on the time-scale of the experiment. Therefore, cells were treated with regulatory volume decrease inhibitors (Lang et al., 1998) NPPB and DCPIB (N+D on the graph) to ensure stable volume increases over the longer time-scale necessary for AFM measurements ($30 min). The measured volume at each time point was normalized to the initial cell volume at t=0 s (Figure 1B). Additionally, to investigate effects of osmotic perturbations on cellular structures, a confocal laser scanning microscope (FV1000, Olympus) was used to image F-actin using the live reporter construct Lifeact-Ruby ( Figures 2D and 2E) and the nucleus using the cell-permeant nuclear marker Hoechst 34332 (Merck-Biosciences).

Attachment of nanobeads to the cell surface and measurement of cell height
HeLa cells were plated onto 50 mm glass bottomed Petri dishes (Fluorodish, World Precision Instruments, Milton Keynes, UK). Yellow-Green carboxylate-modified fluorescent nanobeads with a diameter of 500 nm (FluoSpheres, Molecular Probes, Invitrogen) were coated with collagen-I following the manufacturer's protocol. To attach nanobeads to the cell membrane prior to experiments, Life-Act ruby HeLa cells were incubated for $30 min with a diluted solution containing the fluorescent collagen coated nanobeads. A Piezoelectric z-stage (NanoscanZ, Prior, Scientific, Rockland, MA) was used to calibrate defocused images and to estimate the height of the attached nanobeads ( Figures 1C and 1D). The z-stage was fixed on top of the optical fluorescence microscope stage (IX-71, Olympus) and was piloted using mManager (Micromanager, Palo-Alto, CA) to step up/down while taking epifluorescence images. In order to find the relationship between the z-positions of the bead relative to the plane of focus and the radius of the outer rings formed in the defocused images, we acquired a series of images of the nanobead for different distances from the coverslip with step size of 100 nm. The initial height of the bead (vertical distance between the bottom of cell attached to the coverslip to the bead located on cell membrane) is a critical parameter for estimating the poroelastic diffusion constant using either analytical solutions or FE simulations. The position of the bead relative to the coverslip was estimated by acquiring fluorescence z-stack images of the cells and beads attached to cell with a step size of 200 nm over a distances of $10 mm. Starting a few microns above the cell and ending a few microns below the bottom of the glass coverslip, the z-planes that showed the stress fibres at the bottom of the cell in focus were considered as the bottom of cell (the reference zero height), and the plane in which the bead was in focus was considered as the location of bead and therefore the initial height of the bead was estimated ( Figure 1C).

Principle of defocusing microscopy
Through tracking of fluorescent nanobeads attached to the cell membrane ( Figure 1C), we monitored the response of cells under different hyperosmotic or hypoosmotic perturbations and investigated the swelling/deswelling kinetics of the cytoplasm. Imaging away from the nanobead plane, the bead appears as a set of concentric rings in the image plane ( Figure 1D) and the shape and size of the rings depends on the diffraction pattern of the bead and point spread function of the imaging system (Rosenbluth et al., 2008). Off-plane defocused images of fluorescent nanobeads ( Figure 1D) enabled precise measurements of time-dependent cellular deformations in the z-direction with $10 nm resolution (Perlman and et al., 2013).Furthermore, the defocused images of the nanobead were taken at 100 ms intervals providing an excellent temporal tracking resolution. The spatio-temporal resolution of defocusing technique is significantly higher than spinning disk confocal microscopy enabling precise investigation of cell volume changes at short timescales. We focused on a certain plane that shows a slightly defocused image of the bead and while acquiring images (at 100 ms intervals), either a hypoosmotic or a hyperosmotic shock were applied to drive water in (swelling) or out (deswelling) of the cell, respectively. After application of such osmotic pressure gradients, changes in the positions of the membrane-bound fluorescent nanobeads over time were determined by monitoring the evolution of the radius of the ring in the defocused fluorescent images (Figures 1D, 1E, and 1F).
Image processing to determine changes in cell height from radial projection A 2D defocused image of a bead is circularly symmetric, thus a 1D radius function was used to detect the radius of the outer ring in each defocused image. Because there were several nanobeads attached to each cell, the rough central position of the bead of interest was first selected visually and then the precise coordinates of the centroid of the chosen bead ðx c ; y c Þ were estimated by fitting a 2D Gaussian function to the image. An arbitrary annulus interval D was set to generate a radial distance vector r i = iD; i = 0.N.
By knowing the precise center of the image, the distance of each pixel in the image from the centroid was estimated. According to the radial distance of each pixel d, there were a total of P i pixels lying between radii r iÀ1 and r i + 1 . For each pixel within the annulus between radii r iÀ1 and r i + 1 the averaged radial intensity vector I R ðr i Þ was calculated (see Figure 1E) by splitting the intensity of each pixel I in this annulus linearly according to its radial position d: To find an estimate with sub-pixel resolution for the radial position of the outer rings, the averaged radial intensity curve was normalized to its maximum intensity as shown in the inset of Figure 1E. The intersection of the curve with the line y = 0.05 provided a good estimate for the radial position where the intensity of the outer ring diminishes to zero. The changes in the height of cells versus time were obtained after hypo/hyperosmotic shock were normalized with respect to the final change in cell height. The curves were then used for fitting with analytical solutions and for comparison to FE simulations.

Atomic force microscopy measurements
To estimate the cellular poroelastic parameters including the elastic modulus and poroelastic diffusion constant, we employed atomic force microscopy (AFM) indentation tests (JPK Nanowizard-I, JPK, Instruments, Germany) on central regions of the cells adhered on glass bottom petri dishes. The AFM cantilever (MLCT, Bruker, interfaced with 15 mm latex microspheres, Invitrogen) was pressed into the cells with a fast approach velocity of 10 mm s À1 to excite poroelastic modes of relaxation. Once a maximum force of 5 nN was reached, the cantilever was kept at constant height (via z-closed loop feedback implemented on the JPK Nanowizard) for 5s to enable acquisition of the force-relaxation curves. The force-indentation and force-relaxation curves were analysed to extract the elastic modulus and poroelastic diffusion coefficient (see (Moeendarbary et al., 2013) for detailed methodology and curve fitting procedures).

One-dimensional analytical solution
To study the poroelastic behaviors of living cells, we first evaluated the applicability of one dimensional analytical solution to investigate the swelling/shrinkage kinetics of simple 2D geometries under the framework of linear poroelasticity. We considered the swelling/deswelling kinetics of a linearly isotropic poroelastic material fully immersed in a solvent and bound to an impermeable substrate on its bottom surface ( Figures 3A and S1). Consider a thin layer of a poroelastic material attached from one side ( Figure S1) to a fixed substrate and immersed in a solvent for a long enough time to reach an equilibrium, that is, a state of homogeneity where C 0 and m 0 are the initial concentration and chemical potential of solvent respectively. In the present study, the initial state is assumed to be isotropically swollen from the dry state ( The free energy can be written as a function of the strain tensor in the linear case of an isotropic gel: Combining equations (Equations 12 and 13) results in equation (Equation 2) replacing pore pressure by the term ðm À m 0 Þ=U: Navier and diffusion equations similar to the Biot poroelasticity equations can then be derived using the above equations: These equations can be used to define the dynamics of gel swelling/shrinkage in response to a change in extracellular osmolarity and to extract diffusion coeffcient D. Here we consider the situation where all the components of the strain tensor except ε zz = vuz vz are zero so that the poroelastic material is only allowed to deform in one direction ( Figure S1). In this case the non-zero components of the stress tensor reduce to: 8 > > > < > > > : Combining equations (Equations 6 and 17) and using the fact that ε zz = q, the diffusion equation for the pore pressure is obtained: This equation is an inhomogeneous diffusion equation that can be solved by specifying a stress condition.
In the following the one dimensional time-dependent settlement of a poroelastic material subjected to certain types of boundary conditions is studied. Consider a thin layer of a poroelastic material (a gel) attached from one side ( Figure S1) to a fix substrate and immersed in a solvent for a long enough time to reach an equilibrium homogeneous state with C 0 and m 0 the initial concentration and chemical potential of solvent respectively. When a solvent with different chemical potential of m is suddenly exposed to the gel, the chemical potential gradient m À m 0 drives the solvent to flow into/out of the gel (Yoon et al., 2010). The displacement and chemical potential fields can be considered one dimensional in z if the lateral dimensions of the gel are much larger than its thickness (L x ;L y [H) then there is a negligible amount of flow from the edges and the solvent flows mostly from top surface. The gel is constrained in x-y so ε xx = ε yy = 0. The gel can freely swell and deswell from the top surface (s zz = 0) (Yoon et al., 2010;Bouklas and Huang, 2012) and thus the equations (Equations 17 and 19) can be written in terms of chemical potential: 8 > > > < > > > : iScience Article effects of osmotic changes in the extracellular medium. To investigate the spatial and rise time effects, we first idealised the cell geometry by approximating it to a thin cylindrical disk, with a diameter of 30 mm and a thickness of 5.5 mm, immersed in a solvent sufficiently long to enable it to reach an equilibrium effective pressure with the same value as the outside cell pressure (P eff. = P out -P in = 0). The net pressure inside (P in ) and outside (P out ) of the cell may arise from different combinations of osmotic and hydrostatic pressures. Both free-sliding (free swelling) and non-slip conditions (constrained swelling) for the bottom surface were studied while a no flow condition was always imposed at the bottom surface ( Figure 3A). Later, we considered a more realistic geometry and modeled the cell as an axisymmetric elliptical cap rigidly fixed to a substrate ( Figure 4A). The finite element model of both the idealized and realistic cell was constructed in ABAQUS (version 2018). The geometry and boundary conditions of each discretized model are shown in Figures 3A and 4A, respectively. Given the large mechanical deformation range observed in cell biology (Kang et al., 2008;Lin et al., 2008;Ding et al., 2013;Vargas-Pinto et al., 2013;Nguyen et al., 2016;Zhang and Yang, 2017;Castro et al., 2018), and particularly during our osmotic perturbation experiments, we employed nonlinear-geometry (large displacement formulation) and unstructured mesh options in our FE simulations, and also considered a neo-Hookean isotropic porohyperelastic model (large strain formulation) (Lim et al., 2006;Moeendarbary et al., 2013;Vargas-Pinto et al., 2013;Nguyen et al., 2016;Castro et al., 2018;). The initial shape of HeLa cells was estimated considering 3D stacks of confocal images of a typical cell ( Figure 1A). Averaged elastic modulus, estimated from AFM experiments on the swelled and deswelled cells, were considered in the FE simulations to fit the experimental swelling/deswelling curves and obtain P eff. , D and t r . Our linear isotropic homogeneous assumptions provided a good approximation to capture the dynamics of cell volume changes using a small number of independent parameters.
The implicit ABAQUS/standard finite element solver was applied with automatic time steps of Newtontype iterations. A typical swelling/deswelling simulation included a ramp phase (0 < t < t r ) where the internal pressure was linearly changed up to a certain value over the rise time t r resulting in gradual displacements and change of height at the top surface d(t). The difference between internal and external pressure induces pressure gradients across the cell boundaries and causes fluid to flow from outside to inside, resulting in bulk volumetric swelling/deswelling. For swelling P eff. > 0 / P out > P in was considered, and for deswelling simulations P eff. < 0 / P out < P in was considered. The poroelastic material domain was discretised using the quadratic quadrilateral CPE8P element and the sensitivity of the FE simulations to domain size and mesh element numbers were checked: a mesh convergence study was performed in the ramp phase by decreasing the mesh size until the estimated maximum displacement d N , at the top surface in the finer mesh, differed by less than 0.1% of the d N in the coarser mesh. In the FE solver, a tolerance parameter for maximum pore pressure change (in our case pore pressure being equal to effective pressure as we do not explicitly model osmotic pressure) per increment was used. Because the pore pressure and fully relaxed maximum displacement d N are the main effective output parameters of these simulations, the iteration convergence was achieved by decreasing the pore pressure rate change until the error in the maximum displacement d N was minimized ( Figure S2). For completeness, we also compared linear poroelasticity and neo-Hookean porohyperelasticity. As for the cell-related porohyperelastic parameters (C 10 , D 1 , D, h), they were taken from previous studies (Moeendarbary et al., 2013;Valero et al., 2016) and are summarized in Tables S1 and S2.
To computationally evaluate effects of cell volume changes on the poroelastic parameters (i.e. E, D), we conducted indentation simulations on central regions of a simple constrained geometry as shown in Figure 3. An infinitely rigid indenter of size R = 4.5mm was pressed into the simulated samples (control, swelled and shrunk cases) with a fast approach velocity. The force-indentation and force-relaxation curves were analyzed to extract the elastic modulus and poroelastic diffusion coefficient (see (Esteki et al., 2020) for detailed methodology and (Hu et al., 2010) for curve fitting procedures).