Fluid pumping and active flexoelectricity can promote lumen nucleation in cell assemblies

TL;DR

This paper explores how fluid pumping and active flexoelectricity influence lumen nucleation in cell assemblies, highlighting the coupling of mechanical, electrical, and hydraulic factors in tissue development.

Contribution

It introduces a continuum two-fluid model incorporating fluid pumping and active electric effects to explain lumen formation in tissues.

Findings

Fluid pumping promotes lumen nucleation.

Flexoelectricity significantly influences lumen stability.

Multiple long-term tissue states are possible.

Abstract

We discuss the physical mechanisms that promote or suppress the nucleation of a fluid-filled lumen inside a cell assembly or a tissue. We discuss lumen formation in a continuum theory of tissue material properties in which the tissue is described as a two-fluid system to account for its permeation by the interstitial fluid, and we include fluid pumping as well as active electric effects. Considering a spherical geometry and a polarized tissue, our work shows that fluid pumping and tissue flexoelectricity play a crucial role in lumen formation. We furthermore explore the large variety of long-time states that are accessible for the cell aggregate and its lumen. Our work reveals a role of the coupling of mechanical, electrical and hydraulic phenomena in tissue lumen formation.

Tables

Table 1. Table 1: Experimental values and references (left columns) and estimated values (right columns) of the phenomenological parameters of the model appearing in the constitutive equations.

Parameters	Exp. values	Parameters	Estimations
$\eta$ Forgacs et al. (1998)	${10}^4$ Pa$\cdot$s	$\overline{\kappa }$	${10}^3$ A$\cdot$s$/$m3
$\overline{\eta }$ Montel et al. (2011)	${10}^9$ Pa$\cdot$s	${\lambda }_1$	${10}^8$ N$/$m3
${\gamma }_{1,2}$ Forgacs et al. (1998)	${10}^{-3}$ N/m	${\lambda }_2,{\lambda }_3$	${10}^6$ N$/$m${}_{}{}^2/$V
${\kappa }^{-1}$ Netti et al. (2000)	${10}^{-13}$ m2/Pa/s	${\lambda }_4$	${10}^3$ N$/$m2
${\Pi }_{1,2}^{\mathrm{ext}}$ Prather et al. (1968); *brace1977	${10}^3$ Pa	${\Lambda }_1$	${10}^0$ A$/$m2
$v_{1,2}$ Delarue et al. (2014)	${10}^{-10}$ m/s	${\Lambda }_2,{\Lambda }_3$	${10}^{-2}$ A$/$V/m
$P_h^c$ Montel et al. (2011)	${10}^3$ Pa	${\Lambda }_4$	${10}^{-5}$ A$/$m
${\zeta }_0$ Delarue et al. (2014)	${10}^3$ Pa	${\nu }_1$	${10}^1$ N/m/V
${\zeta }_1$ Delarue et al. (2014)	$-{10}^{-1}$	${\nu }_2$	${10}^8$ Pa$\cdot$s/m
		${\nu }_3$	${10}^0$ N/m/V
		${\nu }_4$	${10}^8$ Pa$\cdot$s/m
		${\Lambda }_1^f$, ${\Lambda }_2^f$	${10}^{-11}$ m/Pa/s
		$J_{p,1}$, $J_{p,2}$	${10}^{-10}$ m/s
		$1-\varphi$	${10}^{-2}$
		${\nu }_0$	$1$

Table 2. Table 2: Dimensionless parameter values used for plotting the figures.

Figure	Parameters values
	${\delta }_1$	${\delta }_2$	$\chi$	${\hat{\gamma }}_0$	${\hat{\gamma }}_1$	${\hat{\gamma }}_2$	$\hat{\lambda }$	${\hat{v}}_1$	${\hat{v}}_2$
3(a)	1	-5	1	0	-0.01	0.02	0	0.1	0.5
3(b)	-1	-0.83	1	0	-${\mathrm{1.3\hspace{0.17em}10}}^{-4}$	${\mathrm{9.6\hspace{0.17em}10}}^{-4}$	0	${\mathrm{8.3\hspace{0.17em}10}}^{-3}$	${\mathrm{8.3\hspace{0.17em}10}}^{-3}$
3(c)	1	-0.1	1	0	-0.05	0.15	-5	0.1	0.1
3(d)	-1	33	1	0	-0.03	0.037	0	0.33	-0.13
3(e)	1	-1	1	0	-0.3	0.4	-4	0.1	0.1
3(f)	1	-0.9	10	0	0.4	0.4	0	0.01	0.3
3(g)	1	0.1	1	0	-${\mathrm{4\hspace{0.17em}10}}^{-3}$	${\mathrm{6\hspace{0.17em}10}}^{-3}$	0	-0.05	0.05
3(h)	1	1	1	0	-0.4	0.6	0	-2	-1
5(b)	1	-1	1	0	-0.218	0.318	-2	0.1	0.1
5(c)	1	-1	1	0	-0.28	0.38	-2	0.1	0.1
5(d)	1	-1	1	0	-1.43	1.53	-2	0.1	0.1

Equations

$$\frac{1}{r^2}\frac{\partial }{\partial r}\left(r^2{\tilde{\sigma }}_{rr}^{\mathrm{c}}\right)-\frac{{\tilde{\sigma }}_{\theta \theta }^{\mathrm{c}}+{\tilde{\sigma }}_{\phi \phi }^{\mathrm{c}}}{r}+\frac{\partial {\sigma }^{\mathrm{c}}}{\partial r}=\frac{\partial P^{\mathrm{f}}}{\partial r},$$

$${\sigma }^{\mathrm{c}}+P_{\mathrm{h}}^{\mathrm{c}}=\bar{\eta }{v}_{\gamma \gamma }^{\mathrm{c}}-{\nu }_0{\tilde{\sigma }}_{\alpha \beta }^{\mathrm{c}}{q}_{\alpha \beta }-{\nu }_1{p}_{\alpha }{E}_{\alpha }-{\nu }_2{p}_{\alpha }(v_{\alpha }^{\mathrm{c}}-v_{\alpha }^{\mathrm{f}}),$$

$${\tilde{\sigma }}_{\alpha \beta }^{\mathrm{c}}=2\eta {\tilde{v}}_{\alpha \beta }^{\mathrm{c}}+\zeta q_{\alpha \beta }-{\nu }_3[E_{\alpha }{p}_{\beta }{]}_{\mathrm{st}}-{\nu }_4[v_{\alpha }{p}_{\beta }{]}_{\mathrm{st}}.$$

$$\zeta ={\zeta }_0-{\zeta }_1({\sigma }^{\mathrm{c}}+P_{\mathrm{h}}^{\mathrm{c}}),$$

$$f_{\alpha }=-\kappa (v_{\alpha }^{\mathrm{c}}-v_{\alpha }^{\mathrm{f}})+{\lambda }_1{p}_{\alpha }+{\lambda }_2{E}_{\alpha }+{\lambda }_3{q}_{\alpha \beta }{E}_{\beta }+{\lambda }_4{\partial }_{\beta }{q}_{\alpha \beta }.$$

$$j_{\alpha }=-\bar{\kappa }(v_{\alpha }^{\mathrm{c}}-v_{\alpha }^{\mathrm{f}})+{\Lambda }_1{p}_{\alpha }+{\Lambda }_2{E}_{\alpha }+{\Lambda }_3{q}_{\alpha \beta }{E}_{\beta }+{\Lambda }_4{\partial }_{\beta }{q}_{\alpha \beta },$$

$$-{\sigma }_{rr}^{\mathrm{c}}(R_1)+P^{\mathrm{f}}(R_1)$$

$$-{\sigma }_{rr}^{\mathrm{c}}(R_2)+P^{\mathrm{f}}(R_2)$$

$$v_1^{\mathrm{f},\mathrm{ext}}-\mathrm{d}{R}_1/\mathrm{d}t$$

$$v_2^{\mathrm{f},\mathrm{ext}}-\mathrm{d}{R}_2/\mathrm{d}t$$

$$v_r^{\mathrm{c}}(R_1)=\mathrm{d}{R}_1/\mathrm{d}t+v_1,$$

$$v_r^{\mathrm{c}}(R_2)=\mathrm{d}{R}_2/\mathrm{d}t-v_2,$$

$$\begin{array}{rl}{P}_{1,2}^{\mathrm{eff}}& ={\Pi }_{1,2}^{\mathrm{ext}}-P_{\mathrm{h}}^{\mathrm{c}}-\frac{J_{\mathrm{p},1,2}}{{\Lambda }_{1,2}^{\mathrm{f}}}\\ & -\frac{2}{3}\left({\zeta }_0{\nu }_0+{\lambda }_4+\frac{{\Lambda }_4\lambda +(3{\nu }_1/2-2{\nu }_0{\nu }_3){\Lambda }_1}{\Lambda }\right),\end{array}$$

$${\gamma }_{1,2}^{\mathrm{app}}={\gamma }_{1,2}\mp 4{\nu }_3{\Lambda }_4/\Lambda ,{\gamma }_0^{\mathrm{app}}=\left(3{\nu }_1/2-2(2+{\nu }_0){\nu }_3\right)({\Lambda }_4/\Lambda ).$$

$${\lambda }^{\mathrm{eff}}={\lambda }_1-{\Lambda }_1\lambda /\Lambda ,$$

$$\frac{\mathrm{d}{r}_1}{\mathrm{d}\hat{t}}$$

$$\frac{\mathrm{d}{r}_2}{\mathrm{d}\hat{t}}$$

$${\hat{\gamma }}_{0,1,2}=\frac{{\gamma }_{0,1,2}^{\mathrm{app}}}{\bar{\eta }{\Lambda }_1^{\mathrm{f}}|P_1^{\mathrm{eff}}|},\hat{\lambda }=\frac{{\lambda }^{\mathrm{eff}}\bar{\eta }{\Lambda }_1^{\mathrm{f}}}{|P_1^{\mathrm{eff}}|},{\delta }_{1,2}=\frac{P_{1,2}^{\mathrm{eff}}}{|P_1^{\mathrm{eff}}|},$$

$$\mathrm{d}{r}_1/\mathrm{d}\hat{t}=f(r_2)-2{\hat{\gamma }}_1/r_1,$$

$$r_1^c=2{\hat{\gamma }}_1/f(r_2),$$

$$4{\nu }_3{\Lambda }_4/\Lambda >{\gamma }_1.$$

$${\eta }_{1}v(r)\left(\frac{1}{r^2}+\frac{1}{L_0^2}+\frac{1}{L_{1}r}\right)+{\eta }_2{v}^{\prime }(r)\left(\frac{1}{L_2}-\frac{1}{r}\right)-{\eta }_3{v}^{\prime \prime }(r)={\tilde{\lambda }}_1+\frac{2\tilde{\zeta }}{r}+\frac{2\tilde{\gamma }}{r^2},$$

$$L_0^{-1}$$

$$L_2^{-1}$$

$${\eta }_1$$

$${\eta }_2$$

$${\eta }_3$$

$$\lambda ={\lambda }_2+\frac{2}{3}{\lambda }_3,\Lambda ={\Lambda }_2+\frac{2}{3}{\Lambda }_3,{\kappa }^{\mathrm{eff}}=\kappa -\bar{\kappa }\frac{\lambda }{\Lambda },$$

$$\tilde{\gamma }=\left(\frac{4}{3}(2+\nu ){\nu }_3-(1+4{\zeta }_1/3){\nu }_1\right)\frac{{\Lambda }_4}{\Lambda },\tilde{\zeta }={\zeta }_0+{\tilde{\lambda }}_4+(2{\nu }_3-{\zeta }_1{\nu }_1)\frac{{\Lambda }_1}{\Lambda },$$

$$\alpha =1-\frac{2}{3}{\nu }_0{\zeta }_1,{\tilde{\lambda }}_1=\alpha \left({\lambda }_1-{\Lambda }_1\frac{\lambda }{\Lambda }\right),{\tilde{\lambda }}_4=\alpha \left({\lambda }_4-{\Lambda }_4\frac{\lambda }{\Lambda }\right).$$

$$-{\eta }_{1}v(r)/r^2+{\eta }_2{v}^{\prime }(r)/r+{\eta }_3{v}^{\prime \prime }(r)+{\tilde{\lambda }}_1+\frac{2\tilde{\zeta }}{r}+\frac{2\tilde{\gamma }}{r^2}=0.$$

$$v(r)=A_1{r}^{{\beta }_1}+A_2{r}^{{\beta }_2}+k_0+k_{1}r+k_2{r}^2,$$

$$k_0=\frac{2\tilde{\gamma }}{{\eta }_1},k_1=\frac{2\tilde{\zeta }}{{\eta }_1-{\eta }_2},k_2=\frac{{\tilde{\lambda }}_1}{{\eta }_1-2({\eta }_2+{\eta }_3)},$$

$${\beta }_{1,2}$$

$$\begin{array}{rl}& \frac{\mathrm{d}{r}_1}{\mathrm{d}\hat{t}}+\frac{3({\hat{v}}_1+\mathrm{d}{r}_1/\mathrm{d}\hat{t})r_1^2}{r_2^3-r_1^3}+\frac{3({\hat{v}}_2-\mathrm{d}{r}_2/\mathrm{d}\hat{t})r_2^2}{r_2^3-r_1^3}={\delta }_1\\ & -\frac{2}{r_1}\left({\hat{\gamma }}_1-{\hat{\gamma }}_0\frac{r_1(r_1+r_2)}{r_1^2+r_2^2+r_1{r}_2}\right)+\hat{\lambda }(r_2-r_1)\frac{r_1^2+2r_1{r}_2+3r_2^2}{4(r_1^2+r_2^2+r_1{r}_2)},\end{array}$$

Full text

Introduction

A fundamental problem in biology is to understand the collective organization of many cells that can give rise to complex structures and morphologies. Such phenomena can be studied either in living embryos or developing organisms, but also in vitro for example in organoid systems that recapitulate morphogenetic processes Lancaster and Knoblich (2014); Simunovic and Brivanlou (2017) or by studying even simpler cell assemblies. In these systems, it is often observed that liquid-filled cavities or lumens appear within cell assemblies Sigurbjörnsdóttir et al. (2014): in respiratory, circulatory, and secretory organs, it is typically an interconnected network of tubular lumens that forms Andrew and Ewald (2010). Alternatively, spherical lumens can also form, such as cysts, acini, alveoli or follicles in mammalian epithelial organs O’Brien et al. (2002); Ferrari et al. (2008). Strikingly, this ability of cell assemblies to self-organize and form internal fluid cavities is maintained in simpler model systems such as organoids and even in multicellular spheroids formed by a single cell type. Examples include Madin–Darby canine kidney cells and mammary epithelial cells (MCF-10A) that are observed to form polarized spherical aggregates with liquid-filled lumen Martín-Belmonte et al. (2008); Debnath et al. (2002).

The formation of a well-delimited cavity surrounded by a cohesive cellular structure has been observed to rely on various mechanisms Lubarsky and Krasnow (2003); Andrew and Ewald (2010); Martín-Belmonte et al. (2008); Debnath et al. (2002). Programmed cell death induced at the structure center for instance leads to lumen formation by a process known as cavitation Lubarsky and Krasnow (2003); Sigurbjörnsdóttir et al. (2014). In addition, lumenogenesis must also rely on cells ability to transport water and ions in a collective fashion to open fluid-filled cavities. This capacity of cells to pump fluid has for instance been quantified in experiments on rabbit corneal epithelia Maurice (1972); Sánchez et al. (2002). The ion pumps that are necessary to generate fluid flows also produce ion flows that can lead to the build up of an electric field across the tissue. Strong experimental evidence indeed supports the presence of a voltage difference across many tissues Cereijido (1978); Josephson (1979); Hay and Geddes (1985). More generally, control of cell proliferation by electric mechanisms has received an increasing attention Blackiston et al. (2009); Cervera et al. (2018); Levin and Martyniuk (2018), and is for instance suspected to play a role in zebrafish fin growth control Daane et al. (2018).

In this paper, we use a continuum theory of radially polarized cell spheroids to reveal key physical mechanisms underlying lumen formation. This coarse-grained approach is especially suited to study the combined effects of fluid permeation, electric fields and currents, as well as mechanical stresses stemming from cell division and death Ranft et al. (2012); Sarkar et al. (2019). We show that lumen formation is an active nucleation problem, governed by tissue growth, fluid pumping and active electric effects. In particular, we discover a surprising role of tissue flexoelectricity in lumen nucleation. Flexoelectricity was first observed as a bending of the nematic order of liquid crystals when an external electric field is applied Meyer (1969); De Gennes and Prost (1974). This also implies that electric fields are generated when liquid crystal orientational order is bent. In tissues, flexoelectricity describes the emergence of electric fields when cell polarity orientation is bent or deformed. We also discuss the state diagram of lumen formation as a function of key parameters and explore the interplay of lumen and spheroid growth at long time.

The structure of this manuscript is as follows. We first present the geometry, boundary conditions and material properties of a permeated tissue in the presence of electric fields. Having derived the dynamical equations for the inner and outer radii of the spheroid, we then focus on lumen nucleation and we show how lumen formation at early time is influenced by pumping and active electric effects. Finally, we explore in Sec. III the long-time states of the spheroid and its lumen.

I Constitutive equations of a permeated tissue in the presence of electric fields

I.1 Tissue geometry and notation

Following Ref. Sarkar et al. (2019), we adopt in this paper a coarse-grained, hydrodynamic description of tissues to study the formation of lumen in a spherical aggregate of cells. In the simplified model we consider here, the tissue is permeated by the interstitial fluid and described in a two-fluid framework. The cell active pumping as well as the electric field and current are introduced in the constitutive equations that describe the tissue properties.

Although spheroids and lumens are found with a large variety of shapes, we consider here a simplified geometry with a spherical symmetry, such that the dynamical quantities only vary along the radial direction $\bm{e}_{r}$. The spheroid has an outer radius $R_{2}$ and encloses a fluid-filled lumen of radius $R_{1}<R_{2}$, as illustrated in Fig. 1. The spherical aggregate is surrounded on the outside (for $r>R_{2}$) by a fluid with a hydrostatic pressure $P^{\mathrm{ext}}_{2}$, and containing osmolites that enter neither the tissue nor the interstitial fluid, such that there exists an osmotic pressure $\Pi^{\mathrm{ext}}_{2}$. We define similarly the hydrostatic and osmotic pressures $P^{\mathrm{ext}}_{1}$ and $\Pi^{\mathrm{ext}}_{1}$ inside the lumen ($r<R_{1}$). The precise boundary conditions involving in particular tissue permeation at the inner and outer surfaces as well as growth rates at the boundaries will be specified later, and we now focus on the bulk equations describing the cell aggregate.

I.2 Permeation of the spheroid by the interstitial fluid

To account for fluid transport and the formation of a fluid-filled cavity inside the spheroid, we decompose the total tissue stress as $\sigma_{\alpha\beta}=\sigma^{\rm c}_{\alpha\beta}+\sigma^{\rm f}_{\alpha\beta}$, where $\sigma^{\rm c}_{\alpha\beta}$ is the stress associated with the cells and $\sigma^{\rm f}_{\alpha\beta}$ is the stress associated with the interstitial fluid. For simplicity, we moreover consider that the anisotropic stress in the interstitial fluid vanishes over length scales large compared to that of the cells. In this case the interstitial fluid flow is driven by pressure gradient and the fluid stress is simply $\sigma^{\rm f}_{\alpha\beta}=-P^{\rm f}\delta_{\alpha\beta}$ Ranft et al. (2012). Neglecting inertia and in the absence of external forces, the force balance in spherical coordinates reads:

[TABLE]

where we have decomposed the cell stress as the sum of an isotropic and an anisotropic part: $\sigma^{\rm c}_{\alpha\beta}=\sigma^{\rm c}\delta_{\alpha\beta}+\tilde{\sigma}^{\rm c}_{\alpha\beta}$, and where $\tilde{\sigma}^{\rm c}_{\theta\theta}=\tilde{\sigma}^{\rm c}_{\varphi\varphi}=-\tilde{\sigma}^{\rm c}_{rr}/2$ because of the spherical symmetry of the cell aggregate.

I.3 Isotropic and anisotropic cell stress in a permeated cell aggregate

We now introduce the constitutive equations for a polar tissue permeated by a fluid and subject to an electric field due to ion transport. Cells are considered polar, i.e., they exhibit a structural anisotropy that can be characterized by a unit polarity vector ${\bm{p}}$. This introduces a nematic order parameter for the cells given by $q_{\alpha\beta}=p_{\alpha}p_{\beta}-\delta_{\alpha\beta}/3$. Notice that we consider in the following that cells are polarized along the radial direction, such that $\bm{p}=\bm{e}_{r}$.

Even though tissues are elastic at shorter time scales, cell division and apoptosis enable them to release stress at longer time scales and to become effectively fluid-like Ranft et al. (2010, 2012). This fluidization of the cells is made possible because cells can probe and sense the local stress and react to it, either by dividing or starting apoptosis, or by generating active stress (for instance in the cytoskeleton) as a result of energy consumption at the molecular level.

An expansion based on symmetry near the homeostatic pressure $P^{\rm c}_{\rm h}$ – defined as the pressure at which cell death exactly compensates cell division in the isotropic state Basan et al. (2009) – allows us to derive a constitutive equation for the isotropic cell stress $\sigma^{\rm c}$ (see Ref. Sarkar et al. (2019) for more details). In the quasistatic limit, the elastic stress has relaxed due to cell division and death, and the constitutive equation for the isotropic part reads111Notice that we have added the term $\nu_{2}$ in Eq. (2) and the term $\nu_{4}$ in Eq. (3) that were not present in Ref. Sarkar et al. (2019) although they are allowed by symmetry.:

[TABLE]

where we have introduced the cell and fluid velocities, $v^{\rm c}_{\alpha}$ and $v^{\rm f}_{\alpha}$. We have also defined the cell strain rate tensor $v^{\rm c}_{\alpha\beta}=(\partial_{\alpha}v^{\rm c}_{\beta}+\partial_{\beta}v^{\rm c}_{\alpha})/2$ and the summation over repeated indices is implied. We have introduced $\bar{\eta}$ the effective bulk viscosity of the tissue due to the response of cell growth to stress, $\nu_{0}$ is a dimensionless coefficient that takes into account the possible dependence of the homeostatic pressure on the anisotropic part of the stress, $\nu_{1}$ characterizes the influence of the electric field on the homeostatic pressure, and $\nu_{2}$ is a coefficient accounting for the effects of the relative motion of the cells and the interstitial fluid to homeostatic pressure Sarkar et al. (2019).

A similar expansion for the traceless anisotropic part of the stress tensor reads:

[TABLE]

where $\tilde{v}^{\rm c}_{\alpha\beta}$ is the traceless part of the cell strain rate tensor, and we have moreover defined the symmetric traceless part of the projection of the cell polarity on the electric field: $[E_{\alpha}p_{\beta}]_{\rm st}\equiv p_{\alpha}E_{\beta}+p_{\beta}E_{\alpha}-(2/3)p_{\gamma}E_{\gamma}\delta_{\alpha\beta}$ and similarly for the projection of the cell polarity on the velocity difference: $[v_{\alpha}p_{\beta}]_{\rm st}\equiv p_{\alpha}(v^{\rm c}_{\beta}-v^{\rm f}_{\beta})+p_{\beta}(v^{\rm c}_{\alpha}-v^{\rm f}_{\alpha})-(2/3)p_{\gamma}(v^{\rm c}_{\gamma}-v^{\rm f}_{\gamma})\delta_{\alpha\beta}$. We introduce $\eta$ the isotropic shear viscosity of the tissue, and we ignore the fourth rank tensor nature of the viscosity. The coefficient $\nu_{3}$ describes the coupling of the electric field to the anisotropic cell stress and $\nu_{4}$ represents the magnitude of the coupling induced by the interstitial fluid flow through the anisotropic cells. The magnitude of the active anisotropic cell stress $\zeta$ can include both a collective component arising from cell division and death, and a contribution from each cell due to the activity of their cytoskeleton Delarue et al. (2014). This active stress can be regulated by the cells and we therefore consider that it depends on the local pressure at linear order as:

[TABLE]

where $\zeta_{0,1}$ are assumed to be constant. Notice that $\zeta_{1}$ is a dimensionless parameter.

I.4 Fluid permeation and electric currents

To close our system of equations, we moreover write a constitutive equation for the momentum exchange $f_{\alpha}=\partial_{\beta}\sigma^{\rm f}_{\alpha\beta}$ between interstitial fluid and cells Ranft et al. (2012). This term, which is balanced by the interstitial fluid pressure gradient $-\partial_{r}P^{\rm f}$ in our geometry, can be expressed as:

[TABLE]

The coefficient $\kappa$ describes the friction due to the (relative) flow of the interstitial fluid in the nanometric cleft between cells, and corresponds to Darcy’s law Darcy (1856) in the description of porous media. The second term on the right-hand side represents the active pumping of fluid by the cells, with $\lambda_{1}$ the active pumping coefficient. The third and fourth terms, proportional to $\lambda_{2}$ and $\lambda_{3}$ respectively, represent the isotropic and anisotropic parts of the force density generated by the electric field. The last term proportional to $\lambda_{4}$ characterizes the sensitivity of the pumping to the bending of the tissue Ramaswamy et al. (2000).

Similarly, we can also write a constitutive equation for the electric current density $j_{\alpha}$:

[TABLE]

where $\bar{\kappa}$ is the coefficient that characterizes the current due to the (relative) flow of ions between cells as a consequence of a reverse electroosmotic effect Kirby (2013). The coefficient $\Lambda_{1}$ characterizes the contribution of ion pumping to the electric current, while $\Lambda_{2}$ and $\Lambda_{3}$ are the isotropic and anisotropic part of the electric conductivity tensor. The coefficient $\Lambda_{4}$ is an active flexoelectric coefficient, indicating that a spatially nonuniform cell polarity orientation is obtained in response to an electric field222Notice that both $\lambda_{4}$ and $\Lambda_{4}$ were already introduced in Sarkar et al. (2019) but they had a vanishing contribution.. This term plays in particular a crucial role in lumen nucleation as we explain in the following.

Finally, assuming that cells and interstitial fluid have the same mass density, and in the limit of an incompressible tissue, that we consider in the following, mass conservation can be rewritten in such a way that the total volume flux is divergence-free Ranft et al. (2012). Considering that there is no fluid flow inside the lumen, the incompressibility yields a relation between the cell velocity and the fluid velocity inside the tissue: $v_{r}^{\rm f}=-\frac{\phi}{1-\phi}v^{c}_{r}$ where we have introduced the cell volume fraction $\phi$, that we assume to be constant in our model. Similarly, the charge conservation in the quasistatic limit $\partial_{\alpha}j_{\alpha}=0$ can be integrated directly in the absence of external current and yields $j_{r}=0$ throughout the tissue.

I.5 Boundary conditions

The spheroid is surrounded by an external fluid both inside (in the lumen), and outside. This fluid exerts a hydrostatic pressure on the tissue that must be balanced by the tissue surface tension and by the total normal stress at the boundaries:

[TABLE]

where we have introduced the inner and outer tissue surface tensions $\gamma_{1}$ and $\gamma_{2}$. Fluid exchange between the spheroid and the outside is driven by osmotic conditions:

[TABLE]

Here, $\Lambda^{\rm f}_{i}$ is the permeability of the interface to water flow. The fluxes $J_{{\rm p},i}=-\Lambda^{\rm f}_{i}(\Pi^{\mathrm{ext},0}_{i}-\Pi^{\mathrm{int},0}_{i})$, with $\Pi^{\mathrm{ext},0}_{i}$ and $\Pi^{\mathrm{int},0}_{i}$ denoting respectively the outside and inside osmotic pressures of osmolites that can be exchanged between external fluid and tissue, can be nonzero as a result of active pumps and transporters that maintain an osmotic pressure difference, and act effectively as water pumps. We have also defined $v^{\rm f,ext}_{1,2}$, the external flows imposed at the inner and outer boundaries. Consistently with the assumption made earlier that there is no flow inside the lumen, we consider in the following that $v^{\rm f,ext}_{1}=v^{\rm f,ext}_{2}=0$.

The normal velocity of the cells at the boundaries has to match the growth of the spheroid radii. An increased cell proliferation in a thin surface layer has been observed in growing spheroids Delarue et al. (2014); Montel et al. (2011); Delarue et al. (2013). Thus, for the sake of generality, we allow for a thin surface layer of cells, both facing outside and to the lumen (see Fig. 1), to have a growth rate that differs from the bulk. The cell velocity boundary conditions then read:

[TABLE]

where $v_{i}=\delta k_{i}n^{\rm c}_{i}e/n^{\rm c}$ with $e$ the thickness of the boundary layers, $n^{\rm c}_{i}$ and $\delta k_{i}$ the cell number density and the cell growth rate in the surface layers, respectively.

II Lumen nucleation in a spherical cell aggregate

II.1 Equations for the dynamics of the spheroid and its lumen

In the previous section we have introduced the bulk equations that describe the properties of the tissue and can be integrated to obtain the cell velocity profile. The values of the different phenomenological parameters that we have defined can be obtained using experimental data and order-of-magnitude estimates (see Table 1, App. C, and Ref. Sarkar et al. (2019)). In particular, for a spheroid which typical radius is of the order of the hundred of micrometers $10^{-6}\lesssim r\lesssim 10^{-3}$ m, our estimates indicate that effects that are relevant at length scales larger than experimentally accessible ones can be neglected to obtain a simpler velocity profile (see App. A for details).

The bulk cell velocity profile together with the boundary conditions introduced in the previous section then allow us to obtain the dynamics of the inner $R_{1}(t)$ and outer $R_{2}(t)$ radii of the spheroid in the quasistatic limit. We obtain two coupled nonlinear differential equations for the spheroid radii, Eqs. (35) and (36) (see App. A). In these equations, six effective parameters are introduced. Two effective pressures:

[TABLE]

where $\Lambda=\Lambda_{2}+2\Lambda_{3}/3$ is an effective conductivity and where $\lambda=\lambda_{2}+2\lambda_{3}/3$. The effective pressure $P^{\rm eff}_{1}$ (resp. $P^{\rm eff}_{2}$) can be seen as a modification of the homeostatic pressure $P^{\rm c}_{\rm h}$ by the external osmotic pressure $\Pi_{1}^{\rm ext}$ (resp. $\Pi_{2}^{\rm ext}$), the pumping flux $J_{{\rm p},1}$ (resp. $J_{{\rm p},2}$) and electric and active contributions. In particular, if all other quantities are kept constant, we observe that a positive effective pressure $P^{\rm eff}_{1}>0$ indicates a stress on the spheroid inner boundary leading to a shrinkage of the tissue, similar to the increase in cell apoptosis due to a pressure larger than the homeostatic pressure in simpler settings Ranft et al. (2010, 2012). Three apparent tension parameters are also introduced:

[TABLE]

The apparent tension $\gamma^{\rm app}_{1}$ (resp. $\gamma^{\rm app}_{2}$) of the inner (resp. outer) surface is a modification of the tissue surface tension $\gamma_{1}$ (resp. $\gamma_{2}$) stemming from the flexoelectric term proportional to $\Lambda_{4}$ – which generates an electric field as a result of curvature – and the field-induced anisotropic stress characterized by $\nu_{3}$. This modification can in principle lead to a negative apparent surface tension for the inner or outer part of the spheroid. A negative inner apparent surface tension ($\gamma^{\rm app}_{1}<0$) will be shown in the following to enhance lumen formation and allow for spontaneous nucleation. The apparent tension parameter $\gamma^{\rm app}_{0}$ is also due to flexoelectricity. It enters the dynamical equations of the inner and outer radii of the spheroid in a similar way as the apparent surface tensions (see Eqs. (35) and (36) in App. A). Finally, an effective pumping coefficient is introduced:

[TABLE]

which is a combination of the active pumping term $\lambda_{1}$ and of an electric contribution $\Lambda_{1}\lambda/\Lambda$ to the pumping due to electroosmotic effects. This term comes as a prefactor of the spheroid thickness $R_{2}-R_{1}$, indicating that the whole tissue acts as a pump.

II.2 Lumen nucleation: a competition between pumping and electric effects

We now discuss lumen formation in spherical cell aggregates. Assuming that the radius of the lumen is small at the early stage of its formation, we can cast the lumen early dynamics into the form of a nucleation problem. In particular, we highlight in the following the crucial role in lumen formation of pumping and of the active flexoelectrity which contributes to the apparent surface tension $\gamma^{\rm app}_{1}$. In the small lumen limit $R_{1}\ll R_{2}$, the equations describing the radii dynamics partially decouple and we obtain the following dimensionless equations:

[TABLE]

where we have introduced dimensionless radii: $r_{i}(\hat{t})=R_{i}(t)/R_{0}$ with $R_{0}=\Lambda^{\rm f}_{1}\bar{\eta}$ and a dimensionless time $\hat{t}=t/\tau_{0}$ with $\tau_{0}=\bar{\eta}/|P^{\rm eff}_{1}|$. We have also introduced the dimensionless parameters:

[TABLE]

and $\chi\!=\!\Lambda^{\rm f}_{1}/\Lambda^{\rm f}_{2}$, $a\!=\!3\hat{v}_{2}-2\hat{\gamma}_{0}$, $b\!=\!3\hat{v}_{2}-2(\hat{\gamma}_{2}+\hat{\gamma}_{0})$, $\hat{v}_{2}\!=\!v_{2}/(\Lambda^{\rm f}_{1}|P^{\rm eff}_{1}|)$. Note that the inner surface growth velocity $v_{1}$ does not contribute to lumen nucleation. Indeed, lumen growth is mainly fed by inward fluid flow, while an increased cell division at the inner surface comes as a lower order effect.

To allow for an analogy with nucleation of a droplet in a fluid, the equation for the dynamics of $r_{1}$ can be rewritten as:

[TABLE]

with $f(r_{2})=\delta_{1}-a/r_{2}+(b-\delta_{2}r_{2}+\hat{\lambda}r_{2}^{2}/4)/(r_{2}(3+\chi r_{2}))+3\hat{\lambda}r_{2}/4$. In this form, we can make an analogy with the nucleation of a droplet in a fluid: lumen nucleation is driven by a competition between the bulk contribution (first term in Eq. (19)) and the surface term (second term in Eq. (19)). To continue the analogy with droplet nucleation, we also introduce a lumen critical radius:

[TABLE]

which is the radius above which a lumen starts growing. Notice that in this nonequilibrium, active system, the lumen critical radius depends on the value of the outer radius $r_{2}$, and therefore on time. Moreover, the (dimensionless) apparent surface tension $\hat{\gamma}_{1}$ can be negative as the result of active flexoelectricity, and in this case a lumen can open spontaneously starting even from a vanishingly small radius.

Let us now be more specific and consider the case where $r_{2}$ is kept fixed, and the volume contribution $f(r_{2})$ is positive. The phase portrait of the system for this case is plotted on the left panel of Fig. 2. Let us first emphasize that the condition $f(r_{2})>0$ can always be satisfied if the effective pumping term $\hat{\lambda}$ is positive (inward pumping) and if $r_{2}$ is sufficiently large: a signature that the aggregate acts collectively for providing fluid to the lumen. According to Eq. (15), a positive effective pumping can be achieved with a positive active fluid pumping coefficient $\lambda_{1}$, indicating an inward pumping, or with a negative ion pumping coefficient $\Lambda_{1}$, indicating an inward pumping of ions, or with both. Because of this pumping, the larger the spheroid (that is, $r_{2}$), the more cells contribute to the inward flow, and therefore the smaller the critical radius $r_{1}^{c}$. This effect is so dramatic that for an arbitrarily small early lumen $r_{1}(t=0)=\varepsilon$, one can always find a large enough spheroid such that $\varepsilon>r_{1}^{c}$ and thus such that the lumen starts growing. Notice also that a positive effective pressure ($P^{\rm eff}_{1}>0$, which implies $\delta_{1}>0$), which is a modification of the homeostatic pressure, indicates an unfavorable environment for cells in the center and is therefore favorable for lumen formation.

Moreover, when the apparent surface tension is positive ($\gamma^{\rm app}_{1}>0$, which implies $\hat{\gamma}_{1}>0$), the usual droplet nucleation picture is preserved in the sense that the system must first perform work against the apparent surface tension (nucleation barrier) to effectively form a growing lumen. However, the value of $\gamma^{\rm app}_{1}=\gamma_{1}-4\nu_{3}\Lambda_{4}/\Lambda$ is a competition between the tissue surface tension $\gamma_{1}$ and the active flexoelectric contribution $4\nu_{3}\Lambda_{4}/\Lambda$, which effectively lowers the tissue surface tension. As $\hat{\gamma}_{1}$ is decreased – which can for instance be achieved by reducing the tissue surface tension at the boundary with the lumen $\gamma_{1}$, or by increasing the tissue flexoelectric coefficient $\Lambda_{4}$ – the cost for nucleation is lowered (black dots on the left panel of Fig. 2) until it eventually vanishes for $\gamma^{\rm app}_{1}\leq 0$ (which implies $\hat{\gamma}_{1}\leq 0$). As a result, the nucleation barrier vanishes as the apparent surface tension changes sign.

The case $f(r_{2})<0$ is also illuminating: in this case, the usual droplet nucleation picture would indicate that the droplet phase is not stable since it does not lower the system energy. The phase portrait for this scenario is plotted on the right panel of Fig. 2: when the apparent surface tension is positive ($\gamma^{\rm app}_{1}>0$, which implies $\hat{\gamma}_{1}>0$), there is no fixed point and the system is always driven to the lumenless state, as one would expect in the droplet nucleation picture. However, as soon as the apparent surface tension becomes negative, a new attractive fixed point exists at a finite radius, meaning that a small lumen forms spontaneously. In this case, the negative surface tension term drives lumen formation even if a lumen is disfavored by the volume term.

Whether the volume term favors a lumen state or not, we find that a lumen spontaneously nucleates whenever its apparent surface tension becomes negative ($\gamma^{\rm app}_{1}<0$), that is when

[TABLE]

This happens when the flexoelectric effects (proportional to $\Lambda_{4}$) overcome the usual tissue surface tension.

For fixed values of $r_{2}$, we have seen that lumen formation can be cast in the form of a nucleation problem. In fact, one can even introduce the (dimensionless) volume $V_{1}=4\pi r_{1}^{3}/3$ and the function $\Psi(r_{1})=4\pi r_{1}^{2}\hat{\gamma}_{1}-(4/3)\pi r_{1}^{3}f(r_{2})$ such that one has $\mathrm{d}V_{1}/\mathrm{d}\hat{t}=-\partial_{r_{1}}\Psi(r_{1})$ and the lumen radius is obtained by minimizing $\Psi$ at fixed $r_{2}$. The final picture is however more subtle because the outer radius $r_{2}$ is time dependent and evolves according to Eq. (17). The curl of $(\mathrm{d}r_{1}/\mathrm{d}\hat{t},\mathrm{d}r_{2}/\mathrm{d}\hat{t})$ does not vanish and the dynamics given by equations (16) and (17) thus does not result from the gradient of an effective potential. Indeed, the nonpotential nature of the dynamics allows for thickness oscillations of the spheroid, as we discuss in the next section.

III Spheroid and lumen dynamics at long time

III.1 Spheroid long-time states

After lumen nucleation, the spheroid and its lumen follow a dynamics that depends on the parameters of the model. The spheroid and its lumen keep evolving according to Eqs. (35) and (36), and the diversity of the long-time scenarii is illustrated in Fig. 3, where the time evolution of the spheroid radius $R_{2}$ and the lumen radius $R_{1}$ are plotted for parameter values illustrating the different regimes.

In the case where lumen nucleation is favored, which can be achieved for instance with a negative apparent inner surface tension $\hat{\gamma}$, corresponding to conditions where the flexoelectric coupling overcomes the inner tissue surface tension, several long-time fates are possible for the cell aggregate. The spheroid and its lumen can grow (see Fig. 3(a)): this is achieved in this example because a positive inner effective pressure $\delta_{1}$ favors lumen growth and a negative outer effective pressure $\delta_{2}$ favors the outer radius growth. Switching the inner effective pressure $\delta_{1}$ to negative values, the lumen can reach a steady state while the spheroid grows (Fig. 3(b)). Alternatively, both lumen and spheroid can reach steady states (Figs. 3(c) and 3(d)). This is for instance seen in the case of an outward effective pumping of fluid ($\hat{\lambda}<0$) which limits the spheroid growth. Under conditions where lumen formation is favored by flexoelectricity while the pumping is directed outwards, the spheroid can even undergo thickness oscillations if the time scales associated with these mechanisms are sufficiently different (Fig. 3(e)).

When the lumen is not favored, it shrinks and eventually vanishes (Fig. 3(f)). In this example the coefficients describing active pumping and flexoelectricy are zero, and the apparent inner surface tension is positive. Figure 3(g) shows an example of a growing spheroid with growing lumen where the lumen radius approaches the outer radius. This could correspond to the formation of a single layer spheroid. Finally, Fig. 3(h) provides an example of spheroid where the outer layer shrinks faster than the lumen, leading to the disappearance of the spheroid when $R_{1}=R_{2}$. This behavior results from a positive outer effective pressure, $\delta_{2}>0$. Similar behavior can also occur when the lumen grows faster than the outer layer, for example if both inner and outer surface growth velocities $\hat{v}_{1,2}$ are negative.

III.2 State diagram of the spheroid dynamics

We now discuss a typical state diagram to illustrate parameter regions where different behaviors displayed in Fig. 3 can be found. In the state diagram shown in Fig. 4, we vary the dimensionless apparent inner surface tension $\hat{\gamma}_{1}$ and the dimensionless effective pumping coefficient $\hat{\lambda}$, while the other parameters are held fixed. We consider the case where the effective pressure $P^{\rm eff}_{1}$ is fixed. The fixed parameters are chosen as follows. We consider a negative $\delta_{2}$ and positive $\delta_{1}$, promoting growth of the outer radius $R_{2}$ and of the lumen $R_{1}$, respectively (see Eqs. (35) and (36)). The values of the surface layer growth velocities $\hat{v}_{1,2}$ are chosen positive to match experimental results Delarue et al. (2014); Montel et al. (2011); Delarue et al. (2013). Their precise values play a minor role in the lumen dynamics, which is mostly controlled by the fluid influx and outflux across the surface of the cellular aggregate. Inner and outer surface permeabilities are chosen to be equal such that $\chi=1$. The apparent tension $\hat{\gamma}_{0}$ is set to 0 as it plays a marginal role in lumen nucleation (see Eqs. (16) and (17)). Finally we choose the sum of the inner and outer tissue surface tensions, such that $\hat{\gamma}_{1}+\hat{\gamma}_{2}=0.1$. These dimensionless values correspond to physical quantities given in Table 1.

For negative apparent inner surface tensions $\hat{\gamma}_{1}<0$, which corresponds to the left half plane of the state diagram, a lumen always forms, see Eq. (16). If in addition the effective pumping $\hat{\lambda}$ is positive (inward pumping), lumen and spheroid grow and there is no steady state. If instead, $\hat{\lambda}$ is negative (outward pumping), steady state spheroids can occur but also oscillations of the spheroid in certain parameter ranges. Onset of spheroid oscillations corresponds to the appearance of a stable limit cycle in the dynamics and is delimited in the diagram by the solid black curve between the yellow and red regions in the lower left quadrant. As discussed in App. B, the system undergoes a Hopf bifurcation when going from the steady-state region to the oscillation region. For positive apparent inner surface tensions $\hat{\gamma}_{1}>0$, lumen tends to shrink with one exception. If in addition $\hat{\lambda}$ is negative, lumen always shrink while the spheroid can be stationary or growing. Finally, if both $\hat{\gamma}_{1}$ and $\hat{\lambda}$ are positive, they have opposite effects on lumen size, and the outcome depends on the initial radii $R_{1}$ and $R_{2}$ and on parameters values. If the initial radius of the lumen is smaller than the critical radius $R_{1}^{c}$ (which is a function of $R_{2}$), the lumen shrinks, while if this initial radius is larger than this critical value, the spheroid and the lumen grows. For sufficiently large outer radius $R_{2}$, a lumen always grows.

The values of $\hat{\gamma}_{1}$ and $\hat{\lambda}$ depend on the physical coefficients describing fluid and ion pumping, flexoelectricity and surface tension. In the case where the effective pressure $P^{\rm eff}_{1}$ is kept constant, the apparent inner surface tension $\hat{\gamma}_{1}$ depends linearly on the tissue inner surface tension $\gamma_{1}$ and on the flexoelectric coupling $\Lambda_{4}$ (see Eq. (14)). Similarly, the effective pumping coefficient $\hat{\lambda}$ depends linearly on the fluid pumping coefficient $\lambda_{1}$ and on the ion pumping coefficient $\Lambda_{1}$ (see Eq. (15)). For instance, an increase in the flexoelectric coupling $\Lambda_{4}$ can result in a sign change from positive to negative of the apparent inner surface tension $\hat{\gamma}_{1}$ and thus promote lumen nucleation. Similarly, changing the direction of fluid or ion pumping from inwards to outwards, which corresponds to a sign change of the coefficients $\lambda_{1}$ or $\Lambda_{1}$, respectively, can result in a sign change of the effective pumping coefficient $\hat{\lambda}$. In the case of positive apparent inner surface tension, this sign change leads to lumen suppression. Conversely in the case of negative apparent inner surface tension, this leads to the existence of a steady state with finite radii or to thickness oscillations.

Notice that spheroid oscillations have been observed both in vitro Ruiz-Herrero et al. (2017) and in vivo Fütterer et al. (2003) in the case of monolayers, but to the best of our knowledge this phenomenon has not been observed in thick spheroids. For monolayered spheroids, oscillations are usually explained by a cycle of growth of the spheroid that builds a stress on the shell, which finally burst open and therefore shrinks due to the outward fluid flow. This is followed by a healing of the hole and the whole process repeats Fütterer et al. (2003); Ruiz-Herrero et al. (2017). For the thick spheroids we consider here, such a bursting process cannot take place, but interestingly oscillatory regimes can still be predicted.

In the state diagram shown in Fig. 4, we find most of the behaviors shown as examples in Fig. 3. However the behavior of a steady-state lumen radius while the spheroid grows indefinitely (Fig. 3(b)) requires a different choice of parameters. Indeed, a necessary condition to find this specific behavior is $\delta_{1}<0$ (together with $\hat{\lambda}=0$, $\delta_{2}>0$ and $\hat{\gamma}_{1}<0$).

IV Conclusion

Our theoretical work has shown that the formation of a fluid-filled lumen in cell assemblies is governed by nucleation equations that resemble those describing the nucleation of a droplet in a fluid. The nucleation of a lumen however depends on additional effects that are fundamentally active: tissue response to mechanical stress, tissue fluid and ion pumping and tissue active flexoelectricity. In the present context, flexoelectricity describe the ability of a polar tissue to generate an electric current when the polarity axis is splayed. This effect could also be observed for nonpolar cells, provided they have an axis of anisotropy. One might expect that such effects are small compared to those generated by uniform polarity, but our findings here show that they can be significant in the nucleation process of lumen. One possible mechanism for generating a flexoelectric current results from the wedge shape of cells in a splayed tissue, leading to a difference in ion pumping on basal and apical sides and thus in electric current generation.

In addition to exhibiting the role of the coupling between mechanical, hydraulic and electrical mechanisms in lumen formation, we have also used our model to explore the role of this coupling in the long-time dynamics of the aggregate and its lumen. In particular, tissue active flexoelectricity, associated with the mechano-electric response of the tissue, is revealed to play a crucial role in early lumen formation as it generates a bulk term in Eq. (14) that acts as a surface tension. Our order-of-magnitude estimations (see Table 1) indicate that this effect is significant, and the flexoelectric contribution could overcome the tissue surface tension, leading to a negative apparent surface tension. Such a negative apparent inner surface tension guarantees the nucleation of a lumen. Note that the negative contribution to the apparent surface tension stems from a bulk stress proportional to $1/R$ rather than from a genuine surface tension. Therefore negative apparent surface tension should not lead to surface shape instabilities.

Similarly, we observe that the fluid pumping is also influenced by electric effects: the effective pumping term that is defined in Eq. (15) contains a contribution that stems directly from active cell pumping, but also an additional one which can be understood as an electroosmotic contribution to the active pumping. Electroosmotic flows are generated when an electric field is applied to a fluid close to a charged surface Marbach and Bocquet (2019) and recent studies suggest that these electroosmotic flows could be for instance dominant in the corneal fluid transport Sánchez et al. (2002); Fischbarg et al. (2017). In our analysis, an inward effective pumping ensures lumen nucleation if the spheroid is large enough. We have moreover shown that active tissue flexoelectricity and pumping can work in concert, in which case lumen formation is facilitated (or prohibited if both effects tend to close the lumen). However, if fluid pumping is directed outwards while flexoelectric effects tend to open a lumen, steady states of lumen and spheroid are observed. This antagonism between flexoelectric and pumping effects can also give rise to electrohydraulic oscillations, which are radically different from the oscillations observed for spherical cell monolayers that rely on a burst and healing mechanism Fütterer et al. (2003); Ruiz-Herrero et al. (2017); Chan et al. (2019)

Because our framework is based on symmetry considerations and does not rely on specific cell-based mechanisms, we expect our results to give a robust qualitative picture of lumen formation in cell assemblies. Below the cell scale, other approaches based on similar principles become relevant Dasgupta et al. (2018). Lumen nucleation in the mouse embryo has for instance been observed to rely on the nucleation and coarsening of multiple micrometer-sized lumens Dumortier et al. (2019). Once these micrometer-sized lumens have fused, our approach should capture its further evolution. Lumen formation in a cell aggregate can also lead to the formation of a monolayer spheroid Andrew and Ewald (2010); Sigurbjörnsdóttir et al. (2014); Datta et al. (2011). Such monolayers arise in our analysis in certain regimes, see Figs. 3(d) and 3(g) for instance.

Finally, the role of mechanical stress on tissue morphogenesis, the importance of electric effects and fluid transport in cells and tissues are well-known facts. However, we have highlighted in this paper the interplay between these effects and how they need to be combined to understand lumen formation in cell assemblies. The potential importance of such an interplay is suggested for instance in a recent work on the zebrafish fin regeneration Daane et al. (2018), which shows the importance of potassium channels in growth phenotypes. This suggests that electric effects may couple to growth process and size control of tissues Blackiston et al. (2009); Cervera et al. (2018); Levin and Martyniuk (2018). Such observations therefore open the door to studies of tissue morphogenesis where hydraulic, electric and mechanical effects are brought together.

Acknowledgments

JP dedicates his contribution to Prof. Peter S. Pershan (Harvard University) on his 85th birthday, remembering his mentoring in soft and condensed matter physics. We thank Keisuke Ishihara and Arghyadip Mukherjee for insightful discussions. CD thanks Szabolcs Horvát for his MaTeX package that was used to produce the figures.

Appendix A Derivation of the dynamics equations for the inner and outer radii

Using the constitutive equations introduced in the main text, we can rewrite the force balance equation (1) as a differential equation on the cell velocity only. The equation we obtain is the starting point of our study and reads:

[TABLE]

where here and in the following we drop the superscript $\rm c$ for the cell velocity. We have introduced the following effective lengths:

[TABLE]

we have also defined effective viscosities:

[TABLE]

and effective parameters:

[TABLE]

Let us now discuss the values of the effective lengths $L_{0}$, $L_{1}$ and $L_{2}$ that appear naturally in Eq. (22). We focus in this paper on the dynamics of a spherical cell aggregate whose typical radius is of the order of the hundred of micrometers: $10^{-6}\lesssim r\lesssim 10^{-3}$ m. Using experimental data and order-of-magnitude estimations (see Ref. Sarkar et al. (2019) and App. C), we are able to give estimations for the effective lengths $L_{0}$, $L_{1}$ and $L_{2}$ appearing in Eq. (22). The permeation length $L_{0}$ is of the order of the millimeter and may become relevant for larger spheroids, while $L_{1}$ is of the order of the centimeter and $L_{2}$ of the order of the decimeter. In the following, we take the limit $r\ll L_{i}$ and therefore neglect the contribution to the dynamics of the terms involving the effective lengths $L_{i}$. In this limit, Eq. (22) becomes

[TABLE]

This equation can be solved using a power-law ansatz:

[TABLE]

where $A_{1}$ and $A_{2}$ are integration constants to be determined by the boundary conditions at $r=R_{1}$ and $r=R_{2}$ (see Eqs. (7)-(12)), while the coefficients $k_{i}$ are obtained by finding a particular solution to Eq. (31) and read:

[TABLE]

and where the exponents $\beta_{1,2}$ are obtained by solving the homogeneous equation and read:

[TABLE]

In order to keep the analysis simpler, we have considered in the main text the limit where the active stress does not depend on the cell pressure, that is $\zeta_{1}\to 0$. We have moreover considered the case where the tissue shear viscosity is small compared to the tissue bulk viscosity ($\eta\ll\bar{\eta}$), which is justified by experimental values of these parameters, indicating $\eta\simeq 10^{4}$ Pa$\cdot$s while $\bar{\eta}\simeq 10^{9}$ Pa$\cdot$s. In this limit, we have in particular $\beta_{1}=-2$ and $\beta_{2}=1$. Using the boundary conditions, we finally obtain the following dimensionless equations for the radii dynamics:

[TABLE]

where we have introduced dimensionless radii: $r_{1}(\hat{t})=R_{1}(t)/R_{0}$ and $r_{2}(\hat{t})=R_{2}(t)/R_{0}$ with $R_{0}=\Lambda^{\rm f}_{1}\bar{\eta}$ and a dimensionless time $\hat{t}=t/\tau_{0}$ with $\tau_{0}=\bar{\eta}/|P^{\rm eff}_{1}|$. The dimensionless parameters have been defined in Eq. (18) (and $\hat{v}_{1}=v_{1}/\Lambda^{\rm f}_{1}|P^{\rm eff}_{1}|$).

Appendix B Bifurcations and thickness oscillations

We discuss here the bifurcations between the different regions of the state diagram displayed in Fig. 4 and indicated by solid black lines (including the $x$ and $y$ axes). In the yellow region of the lower left quadrant the inner and outer radii reach a stable steady state $(R_{1}^{*},R_{2}^{*})$. As one goes from the lower left quadrant to the upper left quadrant, the fixed point existing in the lower quadrant moves progressively to larger values and is eventually sent to infinity as one crosses the $\hat{\lambda}=0$ axis. Going from the lower left quadrant to the lower right quadrant, the value of the fixed point radius $R_{1}^{*}$ is sent progressively to 0 as $\hat{\gamma}_{1}$ increases. The fixed point radius eventually goes to negative $R_{1}^{*}$ values when the $\hat{\gamma}_{1}=0$ axis is crossed. Finally, the onset of spheroid oscillations (solid black curve between the yellow and red regions in the lower left quadrant) corresponds to the appearance of stable limit cycles through a Hopf bifurcation, as we discuss more in the following paragraphs.

We now focus our attention on the spheroid thickness oscillations (see Figs. 3(e) and 5 for illustration) and on the characterization of the bifurcation between these oscillating states and the stable steady states. Spontaneous thickness oscillations are especially interesting in our model as they appear as a fine interplay between pumping and active electric effects (see state diagram, Fig. 4). These electrohydraulic oscillations are indeed possible when flexoelectric effects spontaneously nucleate a lumen, while an outwards pumping of fluid acts for shrinking the cavity. If these two antagonistic processes are not balanced, which is the case when the system lies inside the red region in Fig. 4, then the dynamics obeys a limit cycle around an unstable fixed point and thickness oscillations are predicted. Notice that the existence of such a limit cycle is not surprising, since, as discussed in the main text for a small lumen, the radii dynamical equations have a nonvanishing curl in phase space.

The bifurcation leading to limit cycles and oscillatory solutions in Fig. 4 is found to be a Hopf bifurcation. Indeed, let us define $\mu$ and $\mu^{\prime}$ the two eigenvalues of the Jacobian associated to Eqs. (35) and (36) evaluated at the fixed point of these equations. We also define the control parameter $\Gamma$ (for instance $\Gamma=\hat{\gamma}_{1}$ for a fixed value of $\hat{\lambda}<0$ in Fig. 4), such that the system oscillates for $\Gamma_{c}^{-}\leq\Gamma\leq\Gamma_{c}^{+}$. For $\Gamma<\Gamma_{c}^{-}$ (or $\Gamma>\Gamma_{c}^{+}$), the eigenvalues of the system have a negative real part, such that the system reaches a stable steady-state at long time. Close to the boundary with the oscillating region, the eigenvalues are complex conjugated and the system spirals towards its fixed point. Precisely at the bifurcation ($\Gamma=\Gamma_{c}^{\pm}$), these eigenvalues cross the imaginary axis – a signature of a Hopf bifurcation – and the fixed point becomes unstable. The system is then driven to a limit cycle and displays periodic oscillations in time. We illustrate the crossing of the imaginary axis by the eigenvalues in Fig. 5(a), for which we have chosen the control parameter to be $\Gamma=\hat{\gamma}_{1}$.

To illustrate the behavior of the system when varying $\hat{\gamma}_{1}$, we have plotted the dynamics of the spheroid and of the lumen radii as a function of time in Fig. 5. Various kind of oscillations can be observed while $\hat{\gamma}_{1}$ is varied: small amplitude oscillations around the unstable fixed point as displayed in Fig. 5(b), large amplitude oscillations, strongly nonlinear and during which the lumen almost closes (see Fig. 5(c)), or quasi-sinusoidal oscillations with a very thin thickness of the spheroid, as shown in Fig. 5(d).

Appendix C Estimation of the parameter values

Estimation of the different phenomenological parameters used in this coarse-grained spheroid model is essential in order to study the model in a biologically relevant regime and to simplify analytic computations. Some of the parameters (such as the cell shear and bulk viscosities, the surface tension etc.) have already been estimated in experiments. However, for most of the remaining phenomenological parameters such experimental values are not yet available, and we therefore used order-of-magnitude estimations to obtain them. Experimental and estimated values are gathered in Table 1. Most of the parameters of our model were already estimated in a previous work Sarkar et al. (2019), and we detail here only the estimation of those that were not estimated, namely the coupling of the isotropic and anisotropic stresses to the velocity difference $\nu_{2}$ and $\nu_{4}$, the coupling of curvature and pumping parameter $\lambda_{4}$, and the active flexoelectricity parameter $\Lambda_{4}$.

The coefficient $\nu_{2}$ is estimated by assuming that $\nu_{2}p_{\alpha}(v^{\rm c}_{\alpha}-v^{\rm f}_{\alpha})$ is the stress due to the hydraulic friction force density $\kappa(v^{\rm c}_{\alpha}-v^{\rm f}_{\alpha})$, such that we obtain \nu_{2}\sim\ell\kappa\simeq 10^{8}~{}\text{Pa\cdots/m} where $\ell\sim 10~{}\mu$m is the typical size of a cell. The coefficient $\nu_{4}$ is estimated by assuming that its contribution to the anisotropic cell stress is of the same order of magnitude as the viscous stress in the cleft between cells, that is $\tilde{\sigma}_{xz}\sim\eta^{\rm f}\delta v_{x}/w$ where the interstitial fluid channels width is $w\sim 50$ nm and we have introduced the viscosity of the interstitial fluid $\eta^{\rm f}\sim 100$ mPa$\cdot$s. We then obtain the estimate \nu_{4}\sim\eta^{\rm f}\ell/w^{2}\simeq 10^{8}~{}\text{Pa\cdots/m}.

The coupling of curvature and pumping parameter $\lambda_{4}$ is estimated by noticing that for in a flat geometry, a cell produces a flow $v\sim\kappa^{-1}\lambda_{1}$ due to pumping. Bending this cell implies an extra flow $v^{\prime}\sim v\delta A/A$ where $A$ is the cell area and $\delta A/A\sim\ell C$ (with $C$ the curvature of the bent cell) is the extra area due to bending and which contributes to the extra flow. This extra pumping due to curvature produces a flow which is by definition of the order $v^{\prime}\sim\kappa^{-1}C\lambda_{4}$. We therefore obtain \lambda_{4}\sim\lambda_{1}\ell\sim 10^{3}~{}\text{N\cdotm}^{-2}. The flexoelectricity parameter $\Lambda_{4}$ is obtained using the a similar argument and we obtain \Lambda_{4}\sim\Lambda_{1}\ell\sim 10^{-5}~{}\text{A\cdotm}^{-1}.