Modeling the mechanobioelectricity of cell clusters

We propose a continuum finite strain theory for the interplay between the bioelectricity and the poromechanics of a cell cluster. Specifically, we refer to a cluster of closely packed cells, whose mechanics is governed by a polymer network of cytoskeletal filaments joined by anchoring junctions, modeled through compressible hyperelasticity. The cluster is saturated with a solution of water and ions. We account for water and ion transport in the intercellular spaces, between cells through gap junctions, and across cell membranes through aquaporins and ion channels. Water fluxes result from the contributions due to osmosis, electro-osmosis, and water pressure, while ion fluxes encompass electro-diffusive and convective terms. We consider both the cases of permeable and impermeable cluster boundary, the latter simulating the presence of sealing tight junctions. We solve the coupled governing equations for a one-dimensional axisymmetric benchmark through finite elements, thus determining the spatiotemporal evolution of the intracellular and extracellular ion concentrations, setting the membrane potential, and water concentrations, establishing the cluster deformation. When suitably complemented with genetic, biochemical, and growth dynamics, we expect this model to become a useful instrument for investigating specific aspects of developmental mechanobioelectricity.


Introduction
Recent endeavors have documented that, alongside genetic and biochemical cues, bioelectrical and mechanical signaling is important for development (McCaig et al. 2005;Mammoto and Ingber 2010).
In particular, bioelectricity deals with the study of the ion redistribution within a network of non-excitable cells and its environment, modulating the membrane potential ). The latter qualifies as both a key readout and regulator of several developmental processes, such as proliferation and differentiation (Sundelacruz et al. 2009), at the single cell level, and symmetry breaking (Levin et al. 2002), wound healing and regeneration (Levin 2007), and cancer progression (Chernet and Levin 2013), at the tissue level.
In addition to experimentation through ion channels manipulation, deciphering the bioelectrical dynamics requires ad hoc simulators. The BioElectric Tissue Simulation Engine (BETSE), a finite volume code proposed by Pietak and Levin (2016), allows one to predict the spatiotemporal evolution of the ion concentrations and membrane potential within a cluster of closely packed cells, in response to a perturbation of the bioelectric state. Later, BETSE has been augmented to consider the interplay between genetic, biochemical, and bioelectrical dynamics, so as to explain aspects of planaria regeneration (Pietak and Levin 2017) and developmental brain damage and rescue in frogs (Pai et al. 2018).
Notably, as argued by Silver and Nelson (2018), bioelectrical and mechanical cues are expected to affect each other. Indeed, on the one hand, the membrane potential and the osmotic pressure are strictly related. On the other hand, several ion channels are mechanosensitive, that is, they respond electrically to changes in the membrane mechanics (Martinac 2004).
Then, Silver et al. (2020) have proved that mechanotransduction may effectively direct the establishment of membrane potential gradients within a tissue. In particular, they show that connexin hemichannels, which are mechanosensitive, preferentially open in the peripheral regions of mammary epithelial tissues, characterized by a higher endogenous mechanical stress, thus leading to local depolarization. This, in turn, is responsible for transcriptional changes that promote cell proliferation.
In order to numerically address the coupling between mechanical and bioelectrical signaling, BETSE has been endowed with a solid mechanics module in Leronni et al. (2020), leading to mecBETSE. Specifically, such a module allows one to compute the cluster deformation due to the osmotic pressure gradients determined by the bioelectrical activity and to account for mechanosensitive channels.
However, mecBETSE is limited to small deformations, thus hampering developmental applications. Moreover, it does not account for the water flow triggered by osmotic pressure gradients, which, according to poromechanics frameworks (Coussy 2004), largely employed for cells (Moeendarbary et al. 2013), mechanistically establishes the deformation field.
So as to overcome the previous limits, we propose a continuum finite strain theory coupling the bioelectricity and the poromechanics of cell clusters. In this theory, the bioelectrical response is governed by mass balances for the intracellular and extracellular concentrations of mobile ions and by Gauss laws for the intracellular and extracellular electric potentials. The poromechanical response is determined by mass balances for the intracellular and extracellular concentrations of water molecules and by a momentum balance for the displacement field of the solid network of cytoskeletal filaments and anchoring junctions.
After introducing the object of modeling, that is, the cluster with its main constituents involved in the mechanobioelectrical response, in Sect. 2, we systematically derive the theory in Sect. 3, starting from first principles. In Sect. 4, with reference to a 1D axisymmetric benchmark, we discuss the finite element solution of the proposed model obtained through the commercial software Comsol Multiphysics ® . Finally, in Sect. 5 we draw the conclusions of our study and outline possible future developments.

Modeling object
We refer to a cluster of closely packed animal cells, as sketched in Fig. 1. Each cell is endowed with a cytoskeleton of actin filaments, microtubules, and intermediate filaments. Similar to epithelia, actin and intermediate filaments of neighboring cells are mechanically joined by adherens junctions and desmosomes, respectively (Alberts 1983). Importantly, there exists a thin space between cells (Tsukita et al. 2001), constituting a porosity network referred to as the extracellular (EC) space, within which water and ions can freely flow. At the cluster boundary, we assume that perfectly sealing tight junctions (TJs) (Tsukita et al. 2001) may be either absent or present, respectively allowing or preventing the water and ion exchange between the EC space and the bath surrounding the cluster. Water and ions can also flow directly between neighboring cells through gap junctions (GJs) (Goodenough and Paul 2009;Gao et al. 2011). We refer to the porosity network formed by the cytoplasm and the GJs as the intracellular (IC) space. Finally, the transmembrane water and ion transport, that is, the transport between the IC and EC spaces, is allowed by aquaporins (Agre 2006) and ion channels (Hille 1984), respectively. In a cluster of plant cells, anchoring junctions are replaced by rigid cell walls, and plasmodesmata play the role of GJs (Alberts 1983).

Modeling framework
In the following, we develop a continuum Lagrangian finite strain theory addressing the mechanobioelectricity of the cell cluster described in Sect. 2. By relying on mixture theory (Ateshian 2007), we assume that the solid network and the IC and EC spaces coexist within the same material point, such that transmembrane fluxes should be regarded as local fluxes.

Balance equations
The momentum balance, to be solved for the displacement field , describes the mechanics of the solid network of cytoskeletal filaments and anchoring junctions. Under the quasi-static approximation, and in the absence of bulk loads, it reduces to the mechanical equilibrium where Div is the material divergence and is the nominal stress tensor, such that (Div ) i = P iJ ∕ X J , in which denotes the material position vector and the small and capital case subscripts indicate the spatial and material coordinates, respectively.
The mass balances for the IC and EC water concentrations read where 0 and e 0 are the initial IC and EC porosities, C w and C e w are the IC and EC molar concentrations of water per unit reference volume of the IC and EC spaces, w and e w are the IC and EC nominal molar fluxes of water, J m w is the transmembrane nominal molar flux of water, positive if water moves from the IC to the EC space, A c is the reference cell membrane area, and V c is the reference cell volume. The symbol ̇ indicates time derivative, that is, Ċ w = C w ∕ t . Since the cluster is constituted by closely packed cells, e 0 ≪ 0 . The terms on the right-hand sides are self-balancing, that is, a source of water for the IC space is a sink for the EC space, or vice versa.
Similarly, the mass balances for the IC and EC concentrations of mobile ion i read where C i and C e i are the IC and EC molar concentrations of ion i per unit reference volume of the IC and EC spaces, i and e i are the IC and EC nominal molar fluxes of ion i , and J m i is the transmembrane nominal molar flux of ion i. For the sake of clarity, w and i , e w and e i , and J m w and J m i represent the water and ion fluxes between cells through GJs, in the interconnected intercellular spaces, and across cell membranes through aquaporins and ion channels, respectively.
Finally, the Gauss laws for the IC and EC electric potentials and e read (1) Div = , where and e are the IC and EC nominal electric displacements, F is the Faraday constant, and z i is the valency of ion i . The terms on the right-hand sides represent the IC and EC nominal free charges, and account for both mobile and fixed ions.

Boundary and initial conditions
In our benchmark, we consider that the cluster is traction-free, that is, we supplement Eq. (1) with the boundary condition where is the outward unit normal to the reference boundary.
The initial conditions for the mass balances (2) and (3) read We assume that water and ions can be exchanged with the bath surrounding the cluster through the EC space only, provided that TJs are absent. Therefore, we equip Eqs. (2a) and (3a) with the boundary conditions where the symbol ⋅ denotes the inner product, such that ⋅ = J I N I . As for Eqs. (2b) and (3b), in the absence of TJs we impose chemical equilibrium at the boundary, that is, where e w and e i are the EC chemical potentials of water and ion i , coinciding with those of the bath ̄w and ̄i , in turn supposed to be equal to the initial EC ones. Instead, in the presence of sealing TJs we impose the boundary conditions We assign to Eq. (4a) the boundary condition in which D m is the transmembrane electric displacement at the boundary, which is defined in Sect. 3.6. As for Eq. (4b), in the absence of TJs we impose electrical equilibrium, that is, where ̄ is the electric potential of the bath, assumed to be zero. Instead, in the presence of TJs, we impose with D tj denoting the electric displacement across the TJs, which is defined in Sect. 3.6.
While in the present work we focus on endogenous mechanobioelectricity, our framework may be also adopted to investigate the effect of an externally applied mechanical load or electric field, or of the exposure to an hypotonic or hypertonic environment, by enforcing appropriate boundary conditions.

Thermodynamic restrictions
We follow the approach of Gurtin et al. (2010) for coupled problems of mechanics and species transport, suitably augmented to account for the electric charge of ions. Therefore, under isothermal conditions, the energy balance reads where U is the nominal internal energy density, = + ∇ is the deformation gradient (with denoting the secondorder identity tensor and ∇ denoting the material gradient, such that (∇ ) iJ = u i ∕ X J ), w and i are the IC chemical potentials of water and ion i, are the IC and EC nominal electric fields, and are the IC and EC electrochemical potentials of ion i.
Upon combining Eq. (13) with the second law of thermodynamics and introducing the nominal Helmholtz free energy density W, we obtain the free energy imbalance We assume that W is a function of the independent variables , C w , C i , C e w , C e i , , and e , such that Eq. (16) becomes By relying on the Coleman-Noll procedure, we obtain the following constitutive prescriptions: Consequently, Eq. (17) reduces to the dissipation inequality We remark that the terms in the second line are local dissipation terms, arising from the exchange of water and ions across the cell membrane within the same material point. Since water and ions share the same IC and EC spaces, we assume that each IC or EC flux is affected by the chemical potential gradient of water and by the electrochemical potential gradients of all mobile ions, that is, where the constitutive operators can be collected into the symmetric (Onsager 1931) mobility matrices in which the off-diagonal entries account for so-called crossdiffusion (Vanag and Epstein 2009). Conversely, aquaporins and ion channels are specific for water and ion transport. Therefore, the transmembrane fluxes of water and ion i only depend on the difference between the EC and IC chemical potentials of water and electrochemical potentials of ion i , respectively: We assume the mobility matrices M and M e to be positive definite and the mobility coefficients M m w and M m i to be positive, in order to fulfill Eq. (19), as detailed in Sect. 3.7.

Free energy density
We choose the following additive decomposition for the free energy density: where W mec accounts for the elasticity of the solid network, W mix and W e mix account for the mixing of water and ions in .
the IC and EC spaces, and W pol and W e pol account for the dielectric polarization of the IC and EC spaces.
We adopt for W mec the compressible Neo-Hookean model proposed by Simo and Pister (1984): where and G are the first and second Lamé parameters of the solid network, = T is the right Cauchy-Green deformation tensor, and is the Jacobian, that is, the volume ratio.
We assume the IC and EC solutions of water and ions to be ideal, such that W mix and W e mix are purely entropic and read (Ateshian 2007) where R is the gas constant and T is the absolute temperature. We note that W mix and W e mix account for both mobile and fixed ions, also the latter being part of the IC and EC solutions. We further hypothesize that the IC and EC solutions are dilute, that is,  (Hong et al. 2010) in which 0 is the vacuum permittivity and r is the relative permittivity of the IC and EC solutions, assumed to coincide with that of water given their diluteness.

Constraint on the volume ratio
We assume that the solid network, water, and ions are incompressible, such that the volume ratio of Eq. (25) is inextricably related to the redistribution of water and ions, namely where v w and v i are the molar volumes of water and ion i . In the limit of dilute IC and EC solutions, Eq. (29) reduces to implying that to be replaced in Eq. (6a). Moreover, are the IC and EC molar concentrations of ion i per unit current volume of the IC and EC spaces. In order to impose the constraint (30), we modify the free energy density (23) as (Hong et al. 2010) where p w is a Lagrange multiplier field assuming the role of water pressure. For later developments, we rearrange Eq. (30) for the IC water concentration: Notably, this operation removes C w from the list of the independent variables, in favor of the independent variable p w introduced by Eq. (33).

Conservative constitutive laws
We obtain the nominal stress by combining Eqs. (18a), (24), and (33): where ⊗ denotes the tensor product, such that ( ⊗ ) IJ = D I D J . The stresses w , pol , and e pol could be regarded as active stresses (or eigenstresses), to be balanced by mec through equilibrium (1). The corresponding Cauchy stress is where = T is the left Cauchy-Green deformation tensor and = J −1 and e = J −1 e are the IC and EC current electric displacements. We define the pressure as adopting the convention that each contribution to p is positive if compressive.
By using Eqs. (18b), (27), and (33), we obtain the following IC and EC chemical potentials of water: Similarly, the electric displacements at the boundary across the cell membranes and the TJs in Eqs. (10) and (12)

Dissipative constitutive laws
We choose the following form for the IC and EC mobility matrices of Eqs. (21): in which D w and D e w are the water diffusivities in the IC and EC spaces, while D i and D e i are the diffusivities of ion i in the IC and EC water. As for the transmembrane mobility coefficients of Eqs.
In the perspective of mixture theory, one would obtain the same expressions by relying on the individual momentum balances for water and ions in the IC and EC spaces separately, and assuming that the friction between the different ion species and between the ions and the solid network in either the IC or EC space is negligible, given the diluteness of the solutions (Huyghe and Janssen 1997). Combining Eqs. (38), (42), and (48) provides the following IC and EC fluxes: where the index i refers to mobile ions, whereas the index j refers to fixed ions.
In particular, the first terms in Eqs. (49a) and (49c) account for the water flux down its pressure gradient; they arise because of the coupling of water transport with mechanics. The last terms in Eqs. (49b) and (49d) account for migration, that is, the ion transport in an electric field; they result from the coupling of ion transport with electrostatics. The second terms in Eqs. (49a) and (49c) describe the electro-osmosis of water with mobile ions, while the first terms in Eqs. (49b) and (49d) are associated with the convection of ions with water; electro-osmosis and convection originate from the coupling of water and ion transport. In the absence of ions ( fluxes (49a) and (49c) reduce to Darcy-like fluxes, as in standard poromechanics (Coussy 2004). For immobile water ( D w = D e w = 0 ), the ion fluxes (49b) and (49d) reduce to standard Nernst-Planck fluxes, describing the electro-diffusion of ions (Rubinstein 1990 where is the membrane potential. Equation (50a) accounts for the transmembrane osmosis through aquaporins, whereas Eq. (50b) accounts for the transmembrane electro-diffusion of ions through ion channels, historically addressed through the Goldman-Hodgkin-Katz flux equation (Hille 1984).
Finally, we note that, in light of Eq. (34), substituting Eqs. (49a) and (50a) into Eq. (2a) provides an equation to be solved for the water pressure p w .

One-dimensional axisymmetric benchmark
As a representative benchmark, we consider a circular cell cluster, of reference radius R cl , whose innermost circular region, of reference radius R cl ∕2 and denoted as in , is characterized by a transmembrane diffusivity to sodium D m Na + ten times larger than the surrounding annular region, denoted as out , simulating an overexpression of sodium channels. Given the axial symmetry of the problem, the results only depend on the radial coordinate R. We assume plane stress conditions in Eq. (1), and that each cell is circular in the reference configuration, such that, in Eqs. (2) and (3), A c ∕V c = 2∕R c , with R c denoting the reference cell radius.
We derive the governing equations for this 1D axisymmetric problem in "Appendix A" and detail their finite element implementation in Comsol Multiphysics ® in "Appendix B". After listing the model parameters in Sect. 4.1, we first present the results of the simulation in the absence of both GJs and TJs, in Sect. 4.2. Then, we introduce and comment on the role of GJs in Sect. 4.3. Finally, we further account for TJs in Sect. 4.4.

Parameters
The model parameters are listed in Table 1. We refer to an average animal cell of reference radius R c = 5 m and membrane thickness T m = 5 nm . The cluster reference radius is R cl = 500 m , which is much larger than R c , that is, the characteristic size of a material point, thus ensuring the validity of our continuum formulation. We assume that cells are separated by a reference intercellular space of about 30 nm in size (Pietak and Levin 2016), and that approximately the 70% of the cluster is occupied by fluid. Correspondingly, we obtain an estimate of 0 = 0.695 and e 0 = 0.005 for the initial IC and EC porosities.
The simulations are conducted at body temperature T = 310 K . Accordingly, r = 80 and v w = 18 cm 3 /mol are reliable estimates for the relative permittivity and molar volume of water. We employ m r = 3 for the relative permittivity of the cell membrane (Gramse et al. 2013). The thickness of a TJ complex is about T tj = 500 nm (Tsukita et al. 2001), and we use tj r = 30 for its relative permittivity, which is an average between those of bulk water and proteins inside (Li et al. 2013).
We With reference to a typical mammalian cell, the more abundant ions involved in bioelectricity are sodium, potassium, and chloride. We adopt the following initial IC and EC concentrations (Alberts 1983): C 0 Na + = 10 mol/m 3 , IC Na + , K + , and Cl − diffusivity (with GJs) D i 10 −12 m 2 ∕s IC water diffusivity (with GJs) D w 10 −9 10 −10 ÷ 10 −8 m 2 ∕s C e,0 Na + = 145 mol/m 3 , C 0 K + = 140 mol/m 3 , C e,0 K + = 5 mol/m 3 , C 0 Cl − = 10 mol/m 3 , and C e,0 Cl − = 110 mol/m 3 . We also consider a fixed generic monovalent anion, whose IC and EC concentrations are uniform and constant and equal to C A − = 140 mol/m 3 and C e A − = 40 mol/m 3 . In the IC space, A − is intended to represent negatively charged proteins, nucleic acids, and other cellular constituents. Notably, C A − and C e A − ensure the initial electroneutrality in both the IC and EC spaces (that is, and also the equality of the initial IC and EC osmotic concentrations (that is, C 0 = C e,0 = 300 mol/m 3 ). By using Eq. (31), we obtain C 0 w = C e,0 w = 1∕v w ≈ 5.6 × 10 4 mol/m 3 . Therefore, the IC and EC solutions are actually dilute; indeed, C 0 ∕C 0 w = C e,0 ∕C e,0 w ≈ 0.5%. We employ the following transmembrane ion diffusivities: D m Na + = 10 −18 m 2 ∕s , D m K + = 5 × 10 −18 m 2 ∕s , and D m Cl − = 5 × 10 −17 m 2 ∕s . These are on the order of those reported in Pietak and Levin (2016), but account for the fact that the permeability of artificial lipid bilayers to Na + , K + , and Cl − is not the same (Alberts 1983). As anticipated, in in we set instead D m Na + = 10 −17 m 2 ∕s . We consider that the transmembrane water diffusivity is D m w = 10 −8 m 2 ∕s , that is, ten order of magnitudes larger than D m Na + , as documented in Alberts (1983) with reference to artificial lipid bilayers. In Sect. 4.2, we further analyze the case of smaller D m w , simulating an underexpression of aquaporins.
We follow Pietak and Levin (2016) and assume D e i = 10 −9 m 2 ∕s for the diffusivity of all ions in EC water. We set D e w = 10 −7 m 2 ∕s for the EC water diffusivity, as approximately obtained through the Kozeny-Carman equation (Coussy 2004). Given the uncertainty in this parameter, in Sect. 4.2 we also explore how the response changes by increasing or decreasing D e w of one order of magnitude. In the presence of GJs, we adopt D i = 10 −12 m 2 ∕s for the diffusivity of all ions in IC water. In particular, D i ≤ 10 −14 m 2 ∕s should be excluded, as it has no impact on the behavior of the cluster. Finally, we adopt D w = 10 −9 m 2 ∕s for the IC water diffusivity and, in Sect. 4.3, we further explore how the cluster behavior is affected by variations of D w of one order of magnitude.

Results in the absence of gap and tight junctions
We first assume that GJs are either absent or closed, such that D w = D i = 0 in Eqs. (49a) and (49b). Therefore, the mass balances (2a) and (3a) reduce to ordinary differential equations. Moreover, we assume that TJs are absent, such that boundary conditions (8) and (11) hold. In Fig. 2, we represent the relevant bioelectrical and mechanical fields as a function of R at different times.  Fig. 2 a Relative IC osmotic concentration (C − C 0 )∕C 0 , b relative EC osmotic concentration (C e − C e,0 )∕C e,0 , c IC electric potential , d EC electric potential e , e relative IC water concentration (C w − C 0 w )∕C 0 w , f relative EC water concentration (C e w − C e,0 w )∕C e,0 w , g water pressure p w , h Jacobian J, i areal Jacobian J a , and j radial displacement u as a function of R∕R cl at different times The large D m Na + in in leads to a prominent influx of Na + from the EC to the IC space, that is, down its concentration gradient. Correspondingly, the IC osmotic concentration C rapidly increases. We register a little increase in C in out as well, which is essentially due to the large C A − compared to C e A − , as explained by the Gibbs-Donnan effect (Overbeek 1956). While C presents a steep gradient at R = R cl/2 , due to the lack of GJs connecting in and out , the EC osmotic concentration C e is smoother, because of the interconnection of the intercellular spaces. Moreover, while C e initially decreases with time, it then increases as Na + is transported from the outside to the inside of the cluster.
The redistribution of ions establishes a negative IC electric potential , as of the beginning of the simulation. This is again explained by the Gibbs-Donnan effect. In particular, in is depolarized, that is, at a higher , with respect to out , and the depolarization increases over time due to the influx of Na + . The EC electric potential e remains rather small, such that the membrane potential m nearly corresponds to . In particular, the value of about −60 mV , registered in out , is representative of the resting m associated with the adopted initial ion concentrations and transmembrane diffusivities, which can be estimated through the Goldman-Hodgkin-Katz voltage equation (Hille 1984).
As Na + is transported from the EC to the IC space through ion channels in in , water follows by osmosis through aquaporins. Correspondingly, the IC water concentration C w increases with time. In the EC space, as Na + ions enter the cluster to cope with the request for Na + in in , they drag water molecules by electro-osmosis. Therefore, the EC water concentration C e w increases after initially decreasing, similar to C e .
As water enters in , the water pressure p w increases therein and is equilibrated by the mechanical stress. The IC and EC electrostatic pressures p pol and p e pol , not represented here, are both irrelevant, being orders of magnitude lower than p w . The increase in C w in in is also accompanied by an increase in the Jacobian J. We remark that, given the smaller variation of C e w compared to C w except for the initial transient, and, mostly, the close cell packing, implying 0 ≫ e 0 , the contribution of the variation of C e w to J, as described by Eq. (30), is negligible. The areal Jacobian J a , given by the product of the radial and circumferential deformation gradient components (or stretches) F rR and F , indicates in-plane expansion everywhere, larger in in . By comparing J and J a , we infer that there occurs an out-of-plane expansion in in and a little out-of-plane compression in out . Finally, the radial displacement u increases from R = 0 to R = R cl ∕2 , and then decreases.  Fig. 3 a Relative IC osmotic concentration (C − C 0 )∕C 0 , b relative EC osmotic concentration (C e − C e,0 )∕C e,0 , c IC electric potential , d EC electric potential e , e relative IC water concentration (C w − C 0 w )∕C 0 w , f relative EC water concentration (C e w − C e,0 w )∕C e,0 w , g water pressure p w , h Jacobian J, i areal Jacobian J a , and j radial displacement u as a function of R∕R cl at t = 600 s for different EC water diffusivities D e w 1 3 In Fig. 3, we explore how the same fields are affected by variations of the EC water diffusivity D e w . By increasing D e w , water is transported more rapidly from the outside to the inside of the cluster through the EC space, and consequently from the EC to the IC space across cell membranes. Therefore, fixed the time, C w , C e w , p w , J, J a , and u increase. Interestingly, C and C e increase as well, since the convective contribution to ion transport grows with D e w . The change in the ion redistribution also impacts on m , though mildly. Finally, we note that, by decreasing D e w to 10 −8 m 2 ∕s , the demand for water in in can be barely sustained, such that, for a given time, C w decreases from R = R cl ∕2 to R = 0 ; furthermore, C w also decreases from R = R cl to R = R cl ∕2 , suggesting that water is transported from the IC to the EC space in out , and then from out to in through the EC space, resulting in an in-plane shrinkage of out .
In Fig. 4, we consider the effect of the transmembrane water diffusivity D m w , modulated by the density and open fraction of aquaporins. Increasing D m w above 10 −8 m 2 ∕s or decreasing it up to 10 −12 m 2 ∕s does not affect the results. By further decreasing D m w to 10 −14 m 2 ∕s , the transmembrane water transport is hampered, such that C w strongly decreases and C e w increases. In turn, this determines a decrease in p w , J, J a , and u. Importantly, while varying D e w strongly impacted on the ion redistribution, changing D m w mildly affects it.

Introducing gap junctions
We now investigate on the role of GJs on the mechanobioelectricity of the cluster. As reported in Sect. 4.1, we adopt D i = 10 −12 m 2 ∕s uniformly for all ions and D w = 10 −9 m 2 ∕s. In Fig. 5, we compare the relevant fields at the end of the simulation with and without GJs. If GJs are present, the Na + ions entering the IC space in in flow down their IC electrochemical potential gradient toward out . Therefore, accounting for GJs smooths out the steep gradient of C at R = R cl ∕2 , thus leading to a reduction in C in in and to an increase in C in out . The different ion redistribution in the IC space also influences , with a lesser depolarization occurring in in and a larger one characterizing out . Similarly, the water entering the IC space in in flows toward out through GJs, mainly dragged by ions through electro-osmosis. Predictably, the EC fields are almost no affected by GJs. Notably, the IC water redistribution in the presence of GJs leads to a decrease in J a in in and to an increase in J a in out ; correspondingly, u(R cl ∕2) diminishes, but u(R cl ) remains equal. Without GJs With GJs Fig. 5 a Relative IC osmotic concentration (C − C 0 )∕C 0 , b relative EC osmotic concentration (C e − C e,0 )∕C e,0 , c IC electric potential , d EC electric potential e , e relative IC water concentration (C w − C 0 w )∕C 0 w , f relative EC water concentration (C e w − C e,0 w )∕C e,0 w , g water pressure p w , h Jacobian J, i areal Jacobian J a , and j radial displacement u as a function of R∕R cl at t = 600 s without and with gap junctions  Fig. 6 a Relative IC osmotic concentration (C − C 0 )∕C 0 , b relative EC osmotic concentration (C e − C e,0 )∕C e,0 , c IC electric potential , d EC electric potential e , e relative IC water concentration (C w − C 0 w )∕C 0 w , f relative EC water concentration (C e w − C e,0 w )∕C e,0 w , g water pressure p w , h Jacobian J, i areal Jacobian J a , and j radial displacement u as a function of R∕R cl at t = 600 s for different IC water diffusivities D w 1 3 In Fig. 6, we compare the responses for different values of D w . By increasing D w , more water is transported from in to out , such that the difference in C w between in and out reduces, along with that in J a . However, again, u(R cl ) remains the same. The difference in C between in and out reduces as well, which confirms the relevance of ion transport by convection as D w is risen. Decreasing D w below 10 −10 m 2 ∕s does not affect further the results.
In Fig. 7, we examine the cluster behavior by varying its Young modulus E. Increasing it up to one order of magnitude does not affect the results, except for p w , which grows proportionally to E. Therefore, we conclude that, for suitably small values of E, proper of animal cells, the response of the cluster to the imposed bioelectrical perturbation is independent of E. More specifically, the ion redistribution triggers the water redistribution, which establishes the cluster deformation. However, for larger values of E, which may be proper of plant cells endowed with stiff walls, the mechanics affects the water redistribution as well. Indeed, the accumulation of water in the IC space of in determines the build-up of a large water (turgor) pressure gradient ∇p w , forcing water to flow back toward out , through both GJs and the EC space according to our model [see Eqs. (49a) and (49c)]. Consequently, C w and C e w both diminish, along with J a and u. Though at a lesser extent, C and C e are affected as well, while m practically remains unaltered, given the similar reductions in both and e .
In Fig. 8, we investigate the cluster behavior in the presence of GJs until the steady state, reached at about 24 h . The IC fields monotonically increase with time, both at R = 0 and, albeit slower, as ions and water flow outward through GJs, at R = R cl . At the steady state, the IC fields attain the uniform values (C − C 0 )∕C 0 = (C w − C 0 w )∕C 0 w ≈ 2.5 and = 0 . The EC fields at R = 0 rapidly decrease in the first 2 min and then slowly increase. At the steady state, the EC fields C e − C e,0 , C e w − C e,0 w , and e attain uniform zero values. Therefore, the EC space is undeformed at the steady state. Notably, at the steady state ]; moreover, we could show that c i = c e i = C e,0 i ∀i . The volume ratio J behaves similar to C w , as C e w is not significant for J in the absence of TJs. While initially u(R cl ∕2) ≈ u(R cl ) , they progressively diverge and, at the steady state, u(R cl ) = 2u(R cl ∕2) ≈ 0.4 R cl . To conclude, this simulation reveals that, in the absence of TJs and for a sufficiently compliant cluster [such that p w is irrelevant in Eq. (40)] devoid of ion pumps, at the steady state the current IC and EC ion concentrations and the IC and EC electric potentials become equal to those of the bath surrounding the cluster, in turn coinciding with the initial EC ones. This is accompanied by large cluster deformations, exclusively attributable to the deformation of the IC space.  Fig. 7 a Relative IC osmotic concentration (C − C 0 )∕C 0 , b relative EC osmotic concentration (C e − C e,0 )∕C e,0 , c IC electric potential , d EC electric potential e , e relative IC water concentration (C w − C 0 w )∕C 0 w , f relative EC water concentration (C e w − C e,0 w )∕C e,0 w , g water pressure p w , h Jacobian J, i areal Jacobian J a , and j radial displacement u as a function of R∕R cl at t = 600 s for different Young moduli E, in the presence of gap junctions

Introducing tight junctions
In this section, we comment on the cluster response in the presence of TJs, that is, by considering that the EC space cannot exchange neither water nor ions with the bath surrounding the cluster. Boundary conditions (9) and (12) now hold. As in Sect. 4.3, we also account for GJs. We display the results of the simulation in Fig. 9, by focusing on a relatively short time interval of 30 s.
As for the previous case, the large D m Na + in in leads to a rapid inflow of Na + from the EC to the IC space. However, in the presence of TJs, the ions of the outside bath cannot replace those lost by the EC space. Therefore, given that 0 ≫ e 0 , the increase in C with time is limited, while the decrease in C e is more pronounced and nearly uniform with R.
While e remained nearly null everywhere in the absence of TJs, such that m practically coincided with , here both and e contribute to m , being comparable in magnitude. We observe that, initially, m is very close to the value registered in the absence of TJs. We further note that the positive e at R = R cl corresponds to the transepithelial potential established by TJs (Nuccitelli 2003).
Following Na + , water molecules pass from the EC to the IC space by osmosis through aquaporins, leading to an increase in C w and to a decrease in C e w . As reported for C and C e , given the impermeability of the boundary to water and the large difference between 0 and e 0 , C w little increases, while C e w strongly decreases. Given the limited water redistribution, the mechanical fields are much smaller in magnitude than in the absence of TJs. Furthermore, we highlight that J and consequently p w are negative in out , meaning that the decrease in the EC volume is larger than the increase in the IC volume. Indeed, we remark that, in the presence of TJs, the great disparity between |C w − C 0 w | and |C e w − C e,0 w | makes both contributions important for the estimation of J through Eq. (30). The radial displacement u increases from R = 0 to R = R cl ∕2 , though remaining very small, and then decreases becoming nearly zero at R = R cl .
In Fig. 10, we display the time evolution of the relevant fields until the steady state. After quickly increasing in the first 30 s , C(R = 0) slowly decreases to a steady state value ≈ 1.005 C 0 , reached at about t = 60 min . Similarly, C(R cl ) increases quite rapidly in the first 30 s , but then keeps increasing, though slower, until the same steady state value of C(0). Both C e (0) and C e (R cl ) decrease to the same steady state value ≈ 0.3 C e,0 at about t = 30 s . Therefore, we conclude that Na + ions electro-diffuse from the EC to the IC space of in in the first 30 s and then flow from in to out through GJs until the steady state, when C becomes uniform within the whole cluster.  Fig. 9 a Relative IC osmotic concentration (C − C 0 )∕C 0 , b relative EC osmotic concentration (C e − C e,0 )∕C e,0 , c IC electric potential , d EC electric potential e , e relative IC water concentration (C w − C 0 w )∕C 0 w , f relative EC water concentration (C e w − C e,0 w )∕C e,0 w , g water pressure p w , h Jacobian J, i areal Jacobian J a , and j radial displacement u as a function of R∕R cl at different times, in the presence of gap and tight junctions R = R cl /2 Fig. 10 a Relative IC osmotic concentration (C − C 0 )∕C 0 , b relative EC osmotic concentration (C e − C e,0 )∕C e,0 , c IC electric potential , d EC electric potential e , e relative IC water concentration (C w − C 0 w )∕C 0 w , f relative EC water concentration (C e w − C e,0 w )∕C e,0 w , g water pressure p w , h Jacobian J, and i areal Jacobian J a at R = 0 and R = R cl , and j radial displacement u at R = R cl ∕2 and R = R cl , as a function of time t in the presence of gap and tight junctions The evolution of C w and C e w at R = 0 and R = R cl in the first 30 s is analogous to that observed for C and C e . However, between approximately t = 30 s and t = 3 min , C w (0) increases, while C w (R cl ) decreases. This suggests that, in this time interval, some water flows from out to in , either directly through GJs or by passing through the EC space. After 3 min , water starts flowing back from in to out , until both C w (0) and C w (R cl ) reach the same steady state value ≈ 1.005 C 0 w . Notably, c = C 0 = C e,0 = c e at the steady state; moreover, we could show that c i = c e i ≈ C 0 i ∀i. The Jacobian J and the areal Jacobian J a increase at R = 0 and decrease at R = R cl until t = 3 min and then asymptotically tend to one. In particular, J a is equal to the out-ofplane stretch J∕J a . The radial displacement u at R = R cl ∕2 increases until t = 3 min and then goes to zero.
To conclude, in the presence of TJs and in the absence of ion pumps, at the steady state the current IC and EC ion concentrations are equal and close to the initial IC values, and both the IC and EC spaces are electroneutral. Moreover, within the same material point, the volume increase in the IC space balances the volume decrease in the EC space, such that the cluster is globally undeformed.

Concluding remarks
We have proposed a continuum finite strain theory for the coupling of electrostatics, ion transport, water transport, and mechanics of a closely packed cell cluster.
Specifically, we have regarded the cluster as the superposition of a solid network of cytoskeletal filaments and anchoring junctions and intracellular (IC) and extracellular (EC) solutions of water and ions. We have described the mechanics of the solid network through compressible hyperelasticity. Given the diluteness of the IC and EC solutions and the incompressibility of all the constituents, volumetric deformations are established by water redistribution only. We have obtained the IC and EC fluxes, the first being allowed by gap junctions, by considering cross-diffusing effects. Correspondingly, the IC and EC water fluxes result from the contributions of water pressure and electroosmosis, while the IC and EC ion fluxes are due to electrodiffusion and convection. We have further accounted for transmembrane osmosis and ion electro-diffusion through aquaporins and ion channels, respectively.
We have tested our model to an in-plane circular cluster whose central region in presents an overexpression of sodium channels. The model correctly predicts the accumulation of ions and, consequently, of water, in the IC space of in , and the resulting depolarization and in-plane expansion. The presence of gap junctions smooths out the steep gradients in all the relevant fields otherwise existing at the boundary of in . The contribution of the pressure to the water flux becomes relevant in stiff plant cell clusters, while in deformable animal cell clusters the water flow is exclusively dictated by osmotic phenomena. In the absence of tight junctions, the ion and water redistribution may be severe, leading to large deformations; moreover, the membrane potential and the volumetric deformation are essentially established by the IC fields only, as the EC space remains nearly electroneutral and undeformed, except for the very initial transient. Differently, in the presence of cluster-sealing tight junctions, the ion and water redistribution is much more limited, resulting in small deformations, and both the IC and EC spaces contribute to setting the membrane potential and the volumetric deformation.
The model may be quite straightforwardly complemented with the addition of (i) the active transmembrane ion transport through ion pumps, (ii) the dependency of the diffusivity of ion channels and gap junctions on the membrane potential and tension (so as to account for voltage gating and mechanosensitivity), and (iii) genetic and biochemical dynamics required for specific applications, as already addressed in the literature Levin 2016, 2017;Leronni et al. 2020).
A major and cumbersome advancement would be the inclusion of growth into the model, toward exploring the interplay between mechanical and bioelectrical dynamics in development and regeneration. This could be achieved by introducing an inelastic (growth) deformation gradient, multiplying the elastic contribution and being modulated by the membrane potential, so as to relate growth and depolarization (Sundelacruz et al. 2009;Ambrosi et al. 2019;Silver et al. 2020).
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creat iveco mmons .org/licen ses/by/4.0/.

Governing equations for the one-dimensional axisymmetric benchmark
In a 1D axisymmetric space dimension, equilibrium (1) reduces to where P � rR = P rR ∕ R , and P rR and P are the radial and circumferential nominal stresses, given by [see Eq. (35)] where with and F zZ denoting the radial, circumferential, and out-ofplane deformation gradient components (or stretches), and u (52) being the radial displacement. Under plane stress conditions, the out-of-plane nominal stress P zZ is zero: We equip Eq. (52) with the following boundary conditions [see Eq. (5)], with that in R = 0 ensuing from symmetry: Mass balances (2) and (3) rR 2 0 r JF zZ D 2 + (D e ) 2 = 0 .