Field induced cell proliferation and death in a model epithelium

During the development of an organism from a fertilized egg, tissues are formed by the collective organization of many cells that divide or undergo apoptosis. Tissues grow by repeated rounds of cell division [1]. Cell apoptosis or programmed cell death plays a vital role in maintaining tissue homeostasis, and suppression of apoptosis can result in abnormal cell proliferation which might turn into cancer [2]. The processes of cell division and apoptosis occur during all stages of development and cells removed by apoptosis are often continuously replaced by cell division. Many studies aim to understand how gene regulatory pathways and biochemical signaling are involved in the regulation and coordination of division and apoptosis [3–5]. In recent times the mechanical properties of epithelia have been a topic of great interest, and studies have shown how tissue growth is regulated by local pressure or stiffness [6, 7], and how cell division and apoptosis depend on it [8–10]. The stress at which cell division balances cell apoptosis on average is called the homeostatic stress of the tissue and in that state the tissue can remain stationary [11]. Perturbations from the homeostatic state lead to interesting tissue dynamics that can affect tissue morphogenesis [8, 11]. The concept of stress dependent cell division and apoptosis plays an important role in the present work.

To better understand the mechanical properties of tissues, different theoretical approaches have been developed ranging from individual cell based models [8, 12–14] to hydrodynamic continuum descriptions [15–23]. Coarse grained hydrodynamic models can be used to understand the multicellular dynamics of tissues at long wavelengths. In [24], a one component continuum description of tissues has been developed which takes into account the stress distribution and the flow field generated by cell division and apoptosis and shows that the tissue effectively behaves as a viscoelastic fluid at long time scales. In [25], a two component fluid description of tissues is developed, which considers the cells together with the extracellular matrix as one component and the interstitial fluid as the other. It takes into account the material turnover as a result of cell division and apoptosis explicitly, but also the permeation of the interstitial fluid through a tissue.

The existence of transmembrane and trans-epithelium electric potentials is textbook knowledge [3, 26, 27]. In the case of plasma membranes, pumps and channels maintain an unbalance between the extracellular and intracellular ion densities. Sodium ions are kept preferentially outside cells whereas potassium ions are kept inside by specific pump activity; anions such as chloride ions which maintain electro-neutrality are transported passively by channels; at steady state, this results in a transmembrane potential, the Nernst potential, on the order of 50–100 mV, the negative pole being inside the cells. In monolayer epithelia the pumps and channels are distributed differently between the apical and basal sides. This results in a trans-epithelium potential on the order of 15–60 mV, the basal part being positive compared to the apical side [27, 28]. Thick epithelia often called stratified squamous epithelia also generate potential differences on the order of a few tens of mV [29, 30]. Such tissues are found in many places such as corneal epithelium, lining mucosa of oral cavity, esophagus, anal canal, ectocervix, vagina, foreskin, the internal portion of the lips and the ascending limb of Henle. In general they not only develop an electric potential but they also pump fluid in a physiologically relevant way like for instance in the ascending limb of Henle [29], in rabbit corneal epithelial tissues [31], in insect malpighian tubules [32] and in corneal endothelium [33]. Furthermore, fluid transport play an important role in mechanical signaling during organogenesis [34] and if the fluid flow mechanics is altered, gene expression in cells is affected which leads to disruption in organ development and congenital malfunctions [35–38]. Though a significant amount of experimental evidence is present in the literature, theoretical studies involving fluid pumping and electric field actions in thick epithelia is lacking.

Motivated by the existence of this class of tissues, we study the model case of a slab of finite thickness, homogeneous in perpendicular directions and submitted to the action of both an electric current and a fluid flow as described in figure 1. We consider both the steady state behavior and the growth dynamics. We propose a continuum description based on a two fluid framework. The cells are taken as one component, and the interstitial fluid as the other. Permeation of fluid through the tissue was modeled in [25] in the absence of ion transport. Here we generalize the theory to ion transport and fluid pumping through the epithelium. We study how the fluid flow and the electric field modify the natural homeostasis of the tissue.

This paper is organized as follows. In section 2 we introduce the different conserved quantities relevant for a thick epithelium. The constitutive equations describing tissue material properties and tissue dynamics based on symmetry considerations are introduced in section 3. A model for a thick permeating epithelium is discussed as an example in section 4, where we analyze the dynamics and steady state properties. We discuss our findings in section 5.

We first present the hydrodynamic theory of a such a tissue that is permeated by interstitial fluid flows in the presence of an electric field, starting with a discussion of conservation laws.

2.1. Volume conservation

In our coarse grained approach, we define the cell volume as Ω^c, and the fluid molecular volume as ${{\rm{\Omega }}}^{f}$. The number of cells per unit volume is denoted as n^c, and the number of interstitial fluid molecules per unit volume is denoted as n^f. If the number of cells in a volume V is N^c, and the number of interstitial fluid molecules N^f, then n^c = N^c/V, and n^f = N^f/V, and we have $V={N}^{c}{{\rm{\Omega }}}^{c}+{N}^{f}{{\rm{\Omega }}}^{f}$ [25] or equivalently

$$\begin{eqnarray}&&{n}^{c}{{\rm{\Omega }}}^{c}+{n}^{f}{{\rm{\Omega }}}^{f}=1.\end{eqnarray} \tag{ 1 }$$

The cell volume fraction is defined as ϕ = n^cΩ^c. The volume fraction of the interstitial fluid is then given by ${\phi }^{f}=1-\phi ={n}^{f}{{\rm{\Omega }}}^{f}$.

The cell number density n^c obeys a balance equation which has a flux contribution from cell flows, and also a source/sink contribution from cell division and apoptosis [24],

$$\begin{eqnarray}&&{\partial }_{t}{n}^{c}+{\partial }_{\alpha }({n}^{c}{v}_{\alpha }^{c})={n}^{c}({k}_{d}-{k}_{a}),\end{eqnarray} \tag{ 2 }$$

where ${v}_{\alpha }^{c}$ is the cell velocity, and, k_d and k_a are the rates of cell division and apoptosis respectively. Similarly, the balance equation for interstitial fluid has a flux contribution due to fluid flow and a source/sink contribution. During both cell apoptosis and cell division total mass must be conserved. Assuming the tissue is described by the cell and interstitial densities only, a dying cell must generate interstitial fluid and a dividing cell must consume interstitial fluid. More precisely, a cell volume Ω^c corresponds to Ω^c/Ω^f fluid molecules hence the first term of the right-hand side of equation (3) corresponding to the source-sink term of equation (2). The second term corresponds to the change in the number of interstitial fluid molecules as a result of cell volume change. It takes the simple expression given below, in which the convected derivative is understood to be taken with respect to cell flow, when cells and interstitial fluid have the same mass density and equation (1) is enforced [25]. The complete equation reads:

$$\begin{eqnarray}&&{\partial }_{t}{n}^{f}+{\partial }_{\alpha }({n}^{f}{v}_{\alpha }^{f})=-\displaystyle \frac{{{\rm{\Omega }}}_{c}}{{{\rm{\Omega }}}_{f}}{n}^{c}({k}_{{\rm{d}}}-{k}_{a})-{n}_{c}\displaystyle \frac{{\rm{d}}}{{\rm{d}}{t}}\displaystyle \frac{{{\rm{\Omega }}}_{c}}{{{\rm{\Omega }}}_{f}},\end{eqnarray} \tag{ 3 }$$

where ${v}_{\alpha }^{f}$ is the interstitial fluid velocity.

2.2. Charge conservation

The conservation of charges can be written as

$$\begin{eqnarray}&&\displaystyle \frac{\partial \rho }{\partial t}+{\rm{\nabla }}\cdot {\bf{I}}=0,\end{eqnarray} \tag{ 4 }$$

where ρ(x, t) is the local charge density, and I(x, t) is the electric current.

2.3. Momentum conservation

We consider the total stress to be given ${\sigma }_{\alpha \beta }={\sigma }_{\alpha \beta }^{c}+{\sigma }_{\alpha \beta }^{f}$, where ${\sigma }_{\alpha \beta }^{c}$ is the stress associated with the cells and ${\sigma }_{\alpha \beta }^{f}$ is the stress associated with the interstitial fluid. The cell stress can be decomposed into an isotropic part σ^c and a traceless anisotropic part ${\tilde{\sigma }}_{\alpha \beta }^{c}$, so that the total cell stress can be written as

$$\begin{eqnarray}&&{\sigma }_{\alpha \beta }^{c}={\tilde{\sigma }}_{\alpha \beta }^{c}+{\sigma }^{c}{\delta }_{\alpha \beta }.\end{eqnarray} \tag{ 5 }$$

For simplicity, we consider the anisotropic stress in the interstitial fluid to vanish over length scales large compared to that of the cells. Then the fluid stress is of the form

$$\begin{eqnarray}&&{\sigma }_{\alpha \beta }^{f}=-{P}^{f}{\delta }_{\alpha \beta }.\end{eqnarray} \tag{ 6 }$$

The force balance of the tissue in the absence of any external forces captures momentum conservation and can be written as

$$\begin{eqnarray}&&{\partial }_{\beta }({\sigma }_{\alpha \beta }^{c}+{\sigma }_{\alpha \beta }^{f})=0.\end{eqnarray} \tag{ 7 }$$

The force balance in equation (7), allows for a momentum transfer between the cells and the interstitial fluid which corresponds to an internal force [25]. As a consequence we write

$$\begin{eqnarray}&&{\partial }_{\beta }{\sigma }_{\alpha \beta }^{c}+{f}_{\alpha }=0,\end{eqnarray} \tag{ 8 }$$

$$\begin{eqnarray}&&{\partial }_{\beta }{\sigma }_{\alpha \beta }^{f}-{f}_{\alpha }=0.\end{eqnarray} \tag{ 9 }$$

Here f_α is the force which takes into account the momentum exchance between the two components of the system.

We now introduce, in the spirit of generalized hydrodynamics [39], the constitutive equations of a tissue permeated by a fluid, submitted to the action of an electric field and in which cells are able to generate electric currents and fluid flow. We call such a tissue a polar tissue and to conveniently indicate the directions of pumping and current generation we introduce a 'polarity' vector p with p² = 1. This further implies the existence of a nematic order parameter for the tissue given by ${q}_{\alpha \beta }={p}_{\alpha }{p}_{\beta }-(1/3){p}^{2}{\delta }_{\alpha \beta }$. Furthermore the cells can generate active stress as a result of active processes in the cytoskeleton, such as the action of molecular motors which consume a chemical fuel ATP. This generated stress can additionally have a contribution from cell division and apoptosis [24, 40].

Cell volume and the cell volume fraction are a priori function of local intensive parameters

$$\begin{eqnarray}&&{{\rm{\Omega }}}^{c}={{\rm{\Omega }}}^{c}({\sigma }^{c},{q}_{\alpha \beta }{\widetilde{\sigma }}_{\alpha \beta }^{c},{p}_{\alpha }{E}_{\alpha }),\,{\rm{a}}{\rm{n}}{\rm{d}}\,\,\phi =\phi ({\sigma }^{c},{q}_{\alpha \beta }{\widetilde{\sigma }}_{\alpha \beta }^{c},{p}_{\alpha }{E}_{\alpha }),\end{eqnarray} \tag{ 10 }$$

where Ω_c and ϕ are functions of the isotropic cell stress, σ^c, the projection of the anisotropic cell stress ${\tilde{\sigma }}_{\alpha \beta }^{c}$ on the nematic order parameter q_αβ, and the projection of the applied electric field E_α on the cell polarity p_α. In principle, intensive parameters involve fluid stress as well, in other words fluid pressure. However, to get any sizeable effect that pressure must exceed several tens of atmospheric pressure. The flows we consider here involve fractions of it. Hence, we can safely ignore the fluid pressure dependence here.

The osmotic compressibility of the cells is defined as ${\chi }^{-1}={({n}^{c})}^{-1}({{\rm{d}}{n}}^{c}/{\rm{d}}{\sigma }^{c})$. Because ${n}^{c}=\phi /{{\rm{\Omega }}}^{c}$, we also have $\chi ={n}^{c}{[{\rm{d}}(\phi /{{\rm{\Omega }}}_{c})/{\rm{d}}{\sigma }^{c}]}^{-1}$. The expansion of ${{\rm{d}}{n}}^{c}$/dt to first order using equation (10) can be written as

$$\begin{eqnarray}&&\displaystyle \frac{1}{{n}^{c}}\displaystyle \frac{{{\rm{d}}n}^{c}}{{\rm{d}}t}=-{\chi }^{-1}\displaystyle \frac{{\rm{d}}{\sigma }^{c}}{{\rm{d}}t}-{\chi }_{2}^{-1}\displaystyle \frac{{\rm{d}}({q}_{\alpha \beta }{\widetilde{\sigma }}_{\alpha \beta }^{c})}{{\rm{d}}t}-{\chi }_{3}^{-1}\displaystyle \frac{{\rm{d}}({p}_{\alpha }{E}_{\alpha })}{{\rm{d}}t},\end{eqnarray} \tag{ 11 }$$

where ${\chi }_{2}={n}^{c}{[{\rm{d}}(\phi /{{\rm{\Omega }}}_{c})/{\rm{d}}({q}_{\alpha \beta }{\tilde{\sigma }}_{\alpha \beta }^{c})]}^{-1}$, and ${\chi }_{3}={n}^{c}{[{\rm{d}}(\phi /{{\rm{\Omega }}}_{c})/{\rm{d}}({p}_{\alpha }{E}_{\alpha })]}^{-1}$. Introducing the convective derivative ${{\rm{d}}{n}}^{c}$/${\rm{d}}{t}=\partial {n}^{c}/\partial t+{v}_{\alpha }^{c}{\partial }_{\alpha }{n}^{c}$, equation (2) can then be rewritten as

$$\begin{eqnarray}&&\displaystyle \frac{1}{{n}_{c}}\displaystyle \frac{{{\rm{d}}{n}}_{c}}{{\rm{d}}{t}}=-{v}_{\gamma \gamma }^{c}+{k}_{d}-{k}_{a}.\end{eqnarray} \tag{ 12 }$$

The net growth rate k_d − k_a of the tissue is in general a function of ${\sigma }^{c},{q}_{\alpha \beta }{\tilde{\sigma }}_{\alpha \beta }^{c},{p}_{\alpha }{E}_{\alpha }$, where we focus our attention on variables that are even under time reversal [25] . It can be expressed to first order as:

$$\begin{eqnarray}&&{k}_{d}-{k}_{a}={\bar{\eta }}^{-1}({P}_{h}^{c}+{\sigma }^{c}+\nu {\tilde{\sigma }}_{\alpha \beta }^{c}{q}_{\alpha \beta }+{\nu }_{1}{p}_{\alpha }{E}_{\alpha }),\end{eqnarray} \tag{ 13 }$$

where ${\bar{\eta }}^{-1}={\rm{d}}({k}_{d}-{k}_{a})/{\rm{d}}{\sigma }^{c}{| }_{{\tilde{\sigma }}_{\alpha \beta }^{c},E},\nu$ and ν₁ are expansion coefficients. The homeostatic pressure of the tissue in the absence of anisotropic stress and electric field is denoted P^c_h. The dependence of k_d − k_a on ${q}_{\alpha \beta }{\tilde{\sigma }}_{\alpha \beta }^{c}$ has been introduced in [40] on symmetry grounds, and was shown to be experimentally relevant and measurable. This dependence captures the fact that the anisotropic part of the cell stress tends to change cell shape. It can therefore also influence cellular processes and hence modify the division rate. Similarly, we write the dependence on the p_αE_α term on symmetry grounds. Its physical origin stems from the fact that an external electric field tends to redistribute pumps and channels and also in general interacts with polar microtubules etc. As a result the cell division rate and apoptosis rate are modified. A priori, the value of the corresponding coefficient is not known and the point of this analysis is to provide means to evaluate it experimentally.

Using equations (12) and (13), we obtain a general constitutive equation for the isotropic cell stress

$$\begin{eqnarray}&&\left(1+{\tau }_{l}\displaystyle \frac{{\rm{d}}}{{\rm{d}}{t}}\right)({\sigma }^{c}+{P}_{h}^{c})+\nu \left(1+{\tau }_{\nu }\displaystyle \frac{{\rm{d}}}{{\rm{d}}{t}}\right){\tilde{\sigma }}_{\alpha \beta }^{c}{q}_{\alpha \beta }+{\nu }_{1}\left(1+{\tau }_{1}\displaystyle \frac{{\rm{d}}}{{\rm{d}}{t}}\right){p}_{\alpha }{E}_{\alpha }=\bar{\eta }{v}_{\gamma \gamma }^{c},\end{eqnarray} \tag{ 14 }$$

where ${\tau }_{l}=\bar{\eta }/\chi$ is the isotropic relaxation rate, ${\tau }_{\nu }=\bar{\eta }/{\chi }_{2}$ is the anisotropic relaxation rate, ${\tau }_{1}=\bar{\eta }/{\chi }_{3}$ is the relaxation rate arising from the electric field, and $\bar{\eta }$ is the bulk viscosity. For slowly varying states we can neglect relaxation processes and equation (14) simplifies to

$$\begin{eqnarray}&&{\sigma }^{c}+{P}_{h}^{c}=\bar{\eta }{v}_{\gamma \gamma }^{c}-\nu {\tilde{\sigma }}_{\alpha \beta }^{c}{q}_{\alpha \beta }-{\nu }_{1}{p}_{\alpha }{E}_{\alpha }.\end{eqnarray} \tag{ 15 }$$

Similarly we can write on symmetry grounds a constitutive equation for the anisotropic part of the cell stress which reads to linear order

$$\begin{eqnarray}&&{\tilde{\sigma }}_{\alpha \beta }^{c}=2\eta {\tilde{v}}_{\alpha \beta }^{c}+\zeta {q}_{\alpha \beta }-{\nu }_{2}({p}_{\alpha }{E}_{\beta }+{p}_{\beta }{E}_{\alpha }-\displaystyle \frac{2}{3}{p}_{\gamma }{E}_{\gamma }{\delta }_{\alpha \beta }),\end{eqnarray} \tag{ 16 }$$

where for simplicity we introduce an isotropic shear viscosity η and ν₂ is a coefficient, which couples the electric field to the anisotropic cell stress. The active anisotropic part of the cell stress is given by ζq_αβ. Similarly, we write a constitutive equation for the momentum exchange f_α between interstitial fluid and cells

$$\begin{eqnarray}&&{f}_{\alpha }=-\kappa ({v}_{\alpha }^{c}-{v}_{\alpha }^{f})+{\lambda }_{1}{p}_{\alpha }+{\lambda }_{2}{E}_{\alpha }+{\lambda }_{3}{q}_{\alpha \beta }{E}_{\beta }+{\lambda }_{4}{\partial }_{\beta }{q}_{\alpha \beta },\end{eqnarray} \tag{ 17 }$$

where κ describes friction between the cells and the fluid and κ⁻¹ is the effective permeability of the tissue. The second term results from the pumping activity of the tissue. This is easily seen in considering a situation such that there is no fluid pressure gradient, hence with equations (6) and (9) the momentum exchange vanishes f_α = 0. We furthermore consider a situation with no electric field and no nematic gradient. Under such conditions, there can still exist a fluid flow relative to the cells, which is determined by the balance between the first and the second terms of equation (17): this is the pumping activity observed experimentally. The third and fourth terms represent the isotropic and anisotropic parts of the force density generated by electric fields. They are characterized by the coefficients λ₂ and λ₃ respectively. The fifth term characterized by the coefficient ${\lambda }_{4}$ represents the contribution to f_α arising from a gradient of the nematic order parameter. Finally we can also write a constitutive equation for the electric current

$$\begin{eqnarray}&&{I}_{\alpha }=-\bar{\kappa }({v}_{\alpha }^{c}-{v}_{\alpha }^{f})+{{\rm{\Lambda }}}_{1}{p}_{\alpha }+{{\rm{\Lambda }}}_{2}{E}_{\alpha }+{{\rm{\Lambda }}}_{3}{q}_{\alpha \beta }{E}_{\beta }+{{\rm{\Lambda }}}_{4}{\partial }_{\beta }{q}_{\alpha \beta },\end{eqnarray} \tag{ 18 }$$

where $\bar{\kappa }$ is the coefficient describing the streaming current, Λ₁ the coefficient describing the current resulting from a polar distribution of ion pumps, and Λ₂, Λ₃ are respectively the isotropic and anisotropic part of the electric conductivity tensor. The coefficient Λ₄ is an out-of-equilibrium flexoelectric coefficient. On the coarse-graining scale we use here, the tissue is neutral. The local electric equilibrium is fast compared to tissue dynamics. Thus ∂_tρ ≃ 0, which imposes a condition of conservation of the electric current (18)

$$\begin{eqnarray}&&{\partial }_{\alpha }{I}_{\alpha }=0.\end{eqnarray} \tag{ 19 }$$

Furthermore, assuming the cell and fluid mass densities to be equal and constant [25] implies that the total volume flux ${v}_{\alpha }={n}^{c}{{\rm{\Omega }}}^{c}{v}_{\alpha }^{c}+{n}^{f}{{\rm{\Omega }}}^{f}{v}_{\alpha }^{f}$ is divergence free

$$\begin{eqnarray}&&{\partial }_{\alpha }{v}_{\alpha }=0.\end{eqnarray} \tag{ 20 }$$

The constraint of incompressibility (20) is imposed by using the fluid pressure P^f as a Lagrange multiplier. In a similar way, the constraint of current conservation (19) is imposed by using the electric potential U as the Lagrange multiplier, where ${E}_{\alpha }=-{{\rm{\partial }}}_{\alpha }U$ is the electric field.

We consider a thick planar tissue consisting of cells and interstitial fluid resting on a substrate and embedded in a fluid medium. The fluid surrounds the tissue-substrate system from all sides and can permeate the substrate whereas cells cannot. A constant fluid flow ${v}_{\mathrm{ext}}^{f}$ through the tissue is imposed for instance by an appropriate hydrostatic pressure difference. Similarly, we work at constant imposed electric current density I_ext. We study how the interplay of cell division, cell apotosis, fluid pumping, osmotic pressure and electric current in the tissue controls its dynamics and morphology. We only consider a homogeneous slab, so that dynamical quantities depend only on one variable which we call z and which describes the distance from the substrate. We postpone for future work the study of the stability of the found solutions with respect to variations of quantities parallel to the slab. We consider cells uniformly polarized along the z-direction, and we chose p_z = 1. A schematic representation of the system is given in figure 1.

Zoom In Zoom Out Reset image size

4.1. Constitutive equations in planar geometry

We first write down the constitutive equations for our model tissue in one dimension. Using (15) and (16), the constitutive equations for the cell stress ${\sigma }_{{zz}}^{c}$ can be written as

$$\begin{eqnarray}&&{\sigma }_{{zz}}^{c}=-{P}_{0}^{H}+{\eta }_{\mathrm{eff}}{\partial }_{z}{v}_{z}^{c}-{\nu }_{\mathrm{eff}}{E}_{z},\end{eqnarray} \tag{ 21 }$$

where ${P}_{0}^{H}={P}_{h}^{c}-2\zeta (1-2\nu /3)/3$ is the effective homeostatic pressure of the tissue, defined as the difference between the homeostatic pressure and the total active pressure of the tissue, and ${\eta }_{\mathrm{eff}}=\bar{\eta }+4\eta (1-2\nu /3)/3$ is the effective viscosity of the tissue. Both of them are constants for the given system. Here ${\nu }_{\mathrm{eff}}={\nu }_{1}+4{\nu }_{2}(1-2\nu /3)/3$ is the effective coefficient coupling electric field and stress. Using equations (6), (9), and (17), the force balance equation takes the form

$$\begin{eqnarray}&&{\partial }_{z}{P}_{z}^{f}=\kappa ({v}_{z}^{c}-{v}_{z}^{f})-{\lambda }_{1}-({\lambda }_{2}+\displaystyle \frac{2}{3}{\lambda }_{3}){E}_{z}.\end{eqnarray} \tag{ 22 }$$

Using the current conservation law (19) we have I_z(z) = I_ext, where ${I}_{\mathrm{ext}}$ the externally imposed electric current and the tissue electric field reads

$$\begin{eqnarray}&&{E}_{z}=\displaystyle \frac{{I}_{{\rm{e}}{\rm{x}}{\rm{t}}}-{{\rm{\Lambda }}}_{1}}{{\rm{\Lambda }}}+\displaystyle \frac{\bar{\kappa }({v}_{z}^{c}-{v}_{z}^{f})}{{\rm{\Lambda }}},\end{eqnarray} \tag{ 23 }$$

where Λ = Λ₂ + 2Λ₃/3 is the electric conductivity.

4.2. Boundary conditions

We specify boundary conditions for a planar epithelium of thickness h. We allow for narrow surface layers of thickness e, comparable to cell size, at both z = 0 and z = h to have have growth rates k_d − k_a that differ from the bulk values by a value δk. This gives rise to a difference between cell velocity and interface velocity near the interface. Thus the cell velocity at z = h is related to the interface velocity by

$$\begin{eqnarray}&&{v}_{z}^{c}(z=h)=\displaystyle \frac{{\rm{d}}{h}}{{\rm{d}}{t}}-{v}_{1},\end{eqnarray} \tag{ 24 }$$

where ${v}_{1}=\delta {k}_{1}e$. At z = 0 the tissue is attached to an immobile substrate. The boundary condition for the cell velocity reads

$$\begin{eqnarray}&&{v}_{z}^{c}(z=0)={v}_{2}\end{eqnarray} \tag{ 25 }$$

with ${v}_{2}=\delta {k}_{2}e$.

Incompressibility of the tissue (20) implies $\phi {v}_{z}^{c}+(1-\phi ){v}_{z}^{f}={v}_{\mathrm{ext}}^{f}$, where ${v}_{\mathrm{ext}}^{f}$ is the externally imposed velocity of the external fluid. Therefore the interstitial fluid velocity v_z^f reads

$$\begin{eqnarray}&&{v}_{z}^{f}=\displaystyle \frac{{v}_{\mathrm{ext}}^{f}}{1-\phi }-\displaystyle \frac{\phi }{1-\phi }{v}_{z}^{c}.\end{eqnarray} \tag{ 26 }$$

Force balance at the tissue surface implies stress continuity. Thus the total stress is balanced by the hydrostatic pressure ${P}_{\mathrm{ext}}^{1}$ of the external fluid at z = h

$$\begin{eqnarray}&&{\sigma }_{{zz}}(z=h)={\sigma }_{{zz}}^{c}-{P}^{f}=-{P}_{\mathrm{ext}}^{1}.\end{eqnarray} \tag{ 27 }$$

At z = h fluid is exchanged with the tissue at a rate that is driven by the fluid chemical potential difference between the upper outside region (denoted with the index 1, see figure 1) and the tissue, and can be written as

$$\begin{eqnarray}&&{{\rm{\Lambda }}}^{f}[({P}_{\mathrm{ext}}^{1}-{P}^{f})-{{\rm{\Pi }}}_{\mathrm{ext}}^{1}]+{J}_{p}=\displaystyle \frac{{\rm{d}}{h}}{{\rm{d}}{t}}-{v}_{\mathrm{ext}}^{f}.\end{eqnarray} \tag{ 28 }$$

Here, ${{\rm{\Pi }}}_{\mathrm{ext}}^{1}$ denotes the external osmotic pressure of osmolites that do not to enter the tissue or the interstitial fluid and Λ^f describes the permeability of the interface for water flow. This is for example the case for dextran with molecular mass exceeding 100 kDa [8]. The flux ${J}_{p}={{\rm{\Lambda }}}^{f}({{\rm{\Pi }}}_{\mathrm{ext},0}^{1}-{{\rm{\Pi }}}_{\mathrm{int},0}^{1})$, can be nonzero as a result of active pumps and transporters. It looks like a flow due to an effective water pump. Here ${{\rm{\Pi }}}_{\mathrm{ext},0}^{1}$ and ${{\rm{\Pi }}}_{\mathrm{int},0}^{1}$ denote the outside and inside osmotic pressure, respectively, of osmolites that can be exchanged between external fluid and tissue. Using equation (28), the total cell stress at (z = h) reads

$$\begin{eqnarray}&&{\sigma }_{{zz}}^{c}=-{{\rm{\Pi }}}_{\mathrm{ext}}^{1}+\displaystyle \frac{{J}_{p}}{{{\rm{\Lambda }}}^{f}}+\displaystyle \frac{1}{{{\rm{\Lambda }}}^{f}}\left({v}_{\mathrm{ext}}^{f}-\displaystyle \frac{{\rm{d}}{h}}{{\rm{d}}{t}}\right).\end{eqnarray} \tag{ 29 }$$

Using the above boundary conditions we now discuss the thickness dynamics of a thick epithelium on a substrate or basal membrane.

4.3. Internal dynamics of a thick epithelium

4.3.1. Time dependence of tissue thickness

We first write the equation for the cell velocity v^c_z for the case of constant cell volume fraction ϕ:

$$\begin{eqnarray}&&{\lambda }_{0}^{2}{\partial }_{z}^{2}{v}_{z}^{c}-a{\partial }_{z}{v}_{z}^{c}-{v}_{z}^{c}+{v}_{{\rm{\Lambda }}}=0,\end{eqnarray} \tag{ 30 }$$

where λ₀ with ${\lambda }_{0}^{2}={\eta }_{\mathrm{eff}}(1-\phi )/{\kappa }_{\mathrm{eff}}$ is a hydrodynamic screening length discussed in [25]. The length $a={\nu }_{\mathrm{eff}}\bar{\kappa }/({\rm{\Lambda }}{\kappa }_{\mathrm{eff}})$ stems from the influence of the electric field on tissue flow. The effective fluid pumping velocity is given by ${v}_{{\rm{\Lambda }}}=[\lambda ({I}_{\mathrm{ext}}-{{\rm{\Lambda }}}_{1})/{\rm{\Lambda }}+{\lambda }_{1}](1-\phi )/{\kappa }_{\mathrm{eff}}+{v}_{\mathrm{ext}}^{f}$ and ${\kappa }_{\mathrm{eff}}=\kappa -\lambda \bar{\kappa }/{\rm{\Lambda }}$ is an effective permeability where $\lambda ={\lambda }_{2}+2{\lambda }_{3}/3$.

Equation (30) can be solved for given boundary conditions (29) and (25), to determine the velocity profile and the corresponding cell stress profile

$$\begin{eqnarray}\begin{array}{rcl}{v}_{z}^{c}(z) & = & {v}_{{\rm{\Lambda }}}+({v}_{2}-{v}_{{\rm{\Lambda }}})\exp ({k}_{1}z)\\ & & +(\exp ({k}_{2}z)-\exp ({k}_{1}z))\displaystyle \frac{{\bar{P}}^{H}/{\eta }_{\mathrm{eff}}+{{av}}_{{\rm{\Lambda }}}/{\lambda }_{0}^{2}+(a/{\lambda }_{0}^{2}-{k}_{1})({v}_{2}-{v}_{{\rm{\Lambda }}})\exp ({k}_{1}h)}{({k}_{2}-a/{\lambda }_{0}^{2})\exp ({k}_{2}h)-({k}_{1}-a/{\lambda }_{0}^{2})\exp ({k}_{1}h)},\end{array}\end{eqnarray} \tag{ 31 }$$

$$\begin{eqnarray}\begin{array}{rcl}{\sigma }_{{zz}}^{c} & = & -{\sigma }_{{zz}}^{0}+{\eta }_{\mathrm{eff}}\left[\left({k}_{1}-\displaystyle \frac{{{av}}_{{\rm{\Lambda }}}}{{\lambda }_{0}^{2}}\right)({v}_{2}-{v}_{{\rm{\Lambda }}})\exp ({k}_{1}z)-\displaystyle \frac{{{av}}_{{\rm{\Lambda }}}}{{\lambda }_{0}^{2}}\right]\\ & & +\displaystyle \frac{[{\bar{P}}^{H}+{\eta }_{\mathrm{eff}}{{av}}_{{\rm{\Lambda }}}/{\lambda }_{0}^{2}-{\eta }_{\mathrm{eff}}({k}_{1}-a/{\lambda }_{0}^{2})({v}_{2}-{v}_{{\rm{\Lambda }}})\exp ({k}_{1}h)]}{[({k}_{2}-a/{\lambda }_{0}^{2})\exp ({k}_{2}h)-({k}_{1}-a/{\lambda }_{0}^{2})\exp ({k}_{1}h)]}\\ & & \times [({k}_{2}-a/{\lambda }_{0}^{2})\exp ({k}_{2}z)-({k}_{1}-a/{\lambda }_{0}^{2})\exp ({k}_{1}z)],\end{array}\end{eqnarray} \tag{ 32 }$$

where ${k}_{\mathrm{1,2}}=(a\pm \sqrt{{a}^{2}+4{\lambda }_{0}^{2}})/2{\lambda }_{0}^{2}$. Here ${\bar{P}}^{H}={P}_{h}^{c}-2\zeta (1-2\nu /3)/3+{\nu }_{\mathrm{eff}}[({I}_{\mathrm{ext}}-{{\rm{\Lambda }}}_{1})/{\rm{\Lambda }}-\bar{\kappa }{v}_{\mathrm{ext}}^{f}{\nu }_{\mathrm{eff}}/[{\rm{\Lambda }}(1-\phi )]]-{{\rm{\Pi }}}_{\mathrm{ext}}^{1}+({J}_{p}+{v}_{\mathrm{ext}}^{f})/{{\rm{\Lambda }}}^{f}$ is an effective homeostatic pressure of the tissue which takes into account the effects of electric currents and fluid flows. Equations (31) and (32) show that the instantaneous cell velocity and the cell stress profiles combine two exponentials with characteristic lengths ${k}_{1,2}^{-1}$.

Using the boundary condition (24) for the cell velocity at z = h, we obtain a dynamical equation for the tissue thickness:

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{{\rm{d}}{h}}{{\rm{d}}{t}} & = & \displaystyle \frac{{v}_{{\rm{\Lambda }}}}{(1+{{\rm{\Lambda }}}_{h})}-\displaystyle \frac{[\exp ({k}_{2}h)-\exp ({k}_{1}h)]}{(1+{{\rm{\Lambda }}}_{h})}\displaystyle \frac{[({k}_{1}-a/{\lambda }_{0}^{2})({v}_{2}-{v}_{{\rm{\Lambda }}})\exp ({k}_{1}h)-{{av}}_{{\rm{\Lambda }}}/{\lambda }_{0}^{2}-{\bar{P}}^{H}/{\eta }_{\mathrm{eff}}]}{[({k}_{2}-a/{\lambda }_{0}^{2})\exp ({k}_{2}h)-({k}_{1}-a/{\lambda }_{0}^{2})\exp ({k}_{1}h)]}\\ & & +\displaystyle \frac{({v}_{2}-{v}_{{\rm{\Lambda }}})\exp ({k}_{1}h)}{(1+{{\rm{\Lambda }}}_{h})},\end{array}\end{eqnarray} \tag{ 33 }$$

where ${{\rm{\Lambda }}}_{h}=[\exp ({k}_{2}h)-\exp ({k}_{1}h)]/[{{\rm{\Lambda }}}^{f}{\eta }_{\mathrm{eff}}(({k}_{2}-a/{\lambda }_{0}^{2})\exp ({k}_{2}h)-({k}_{1}-a/{\lambda }_{0}^{2})\exp ({k}_{1}h))]$ is non-zero and positive.

Equation (33) can give rise to three different behaviors, tissue growth, tissue collapse or finite steady state thickness. If on one hand cell division dominates over cell apoptosis dh/dt is positive and tissue growth takes place. If dh/dt remains positive at all times, the thickness increases indefinitely, leading to a complete invasion of the available space. This corresponds to what is commonly called tissue proliferation. If on the other hand apoptosis dominates over cell division, then dh/dt is negative, the tissue shrinks. If dh/dt < 0 persists at all times, the tissue finally collapses. The tissue can also reach a stable steady state with dh/dt = 0 and constant thickness. The steady state thickness then depends on the imposed electric current and fluid flow. Steady states can also be unstable. In this case the tissue thickness evolves away from the steady state value. This can give rise to either tissue collapse, indefinite growth or finite thickness, depending on initial conditions and parameter values.

For simplicity, we restrict our analysis to the case a ≪ λ₀. This is motivated by simple estimates of $\bar{\kappa }\simeq {10}^{3}\,{\rm{C}}\,{{\rm{m}}}^{-3},{\nu }_{{\rm{e}}{\rm{f}}{\rm{f}}}\simeq 1.5$ Pa m V⁻¹, and ${\rm{\Lambda }}\simeq 6\times {10}^{-3}{{\rm{\Omega }}}^{-1}\,{{\rm{m}}}^{-1}$ for which $a\simeq 5\cdot {10}^{-8}$ m. We also estimate λ₀ ≃ 1 cm, see appendix. For a ≪ λ₀, we have ${k}_{1}=-{k}_{2}=1/{\lambda }_{0}$. In this case the time dependence of tissue thickness simplifies to

$$\begin{eqnarray}&&\displaystyle \frac{(1+{{\rm{\Lambda }}}_{h})}{{v}_{{\rm{\Lambda }}}}\displaystyle \frac{{\rm{d}}{h}}{{\rm{d}}{t}}=f(h)=1-\displaystyle \frac{\alpha }{\cosh (h/{\lambda }_{0})}+\beta \tanh (h/{\lambda }_{0}),\end{eqnarray} \tag{ 34 }$$

where $\alpha =1-{v}_{2}/{v}_{{\rm{\Lambda }}}$, and $\beta ={\bar{P}}^{H}{\lambda }_{0}/({\eta }_{\mathrm{eff}}{v}_{{\rm{\Lambda }}})$.

We show in figure 2 examples of the four generic scenarios characterizing the dependence of dh/dt on tissue thickness h. Zeros of f(h) correspond to steady states which can be stable or unstable depending on the sign of the slope df/dh. Negative slopes correspond to stable and positive slopes to unstable situations. We find that either two, one or zero steady states may exist, depending on the values of α and β. The different possible scenarios are summed up in two state diagrams shown in figures 3(a) and (b), for v_Λ > 0, and v_Λ < 0 respectively. The lines delineating the states characterize continuous and discontinuous transitions. Spinodal lines are also shown. Discontinuous transitions, are reached when the conditions f(h) = 0 and df/dh = 0 are satisfied simultaneously. The corresponding line is a portion of the circle α² + β² = 1 as shown on figure 2. It will be shown as a dashed line on all figures. All other boundaries correspond either to continuous transitions shown as solid lines, or to 'spinodal' transitions shown as dotted lines. The 'spinodal' lines signal the appearance/disappearance of an unstable (resp: stable) steady state, when crossing the line from a phase in which only one stable (resp: unstable) steady state exists. In the corresponding domain, initial conditions determine the eventual fate of the system.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

As the coefficients α and β depend on the parameters defining the tissue properties and experimental conditions, we can transform the generic state diagrams corresponding to conditions in which the external flow ${v}_{\mathrm{ext}}^{f}$ or the current I_ext are imposed. We can also investigate the role of parameters such as the effective homeostatic pressure ${\bar{P}}^{H}$. We discuss the steady states and the different growth regimes as a function of ${v}_{\mathrm{ext}}^{f}$ and I_ext in the next subsection.

4.3.2. Steady states

When a steady state with dh/dt = 0 in equation (34) exists, the steady state tissue thickness is given by

$$\begin{eqnarray}&&h={\lambda }_{0}\mathrm{ln}\left[\displaystyle \frac{-B\pm \sqrt{{B}^{2}-4{AC}}}{2A}\right].\end{eqnarray} \tag{ 35 }$$

This thickness is of order the screening length λ₀ up to a logarithmic factor. Here

$$\begin{eqnarray}&&A=\alpha \left[1-\displaystyle \frac{{\lambda }_{0}{\nu }_{\mathrm{eff}}\bar{\kappa }}{{\eta }_{\mathrm{eff}}(1-\phi ){\rm{\Lambda }}}-\displaystyle \frac{{\lambda }_{0}}{{\eta }_{\mathrm{eff}}{{\rm{\Lambda }}}^{f}}\right]+{\gamma }_{1}+{\gamma }_{2}+\beta \left[1+\displaystyle \frac{{\lambda }_{0}{\nu }_{\mathrm{eff}}{\kappa }_{\mathrm{eff}}}{\lambda (1-\phi ){\eta }_{\mathrm{eff}}}\right],\end{eqnarray} \tag{ 36 }$$

$$\begin{eqnarray}&&B=2(1-\alpha -\beta -{\gamma }_{1}),\end{eqnarray} \tag{ 37 }$$

$$\begin{eqnarray}&&C=\alpha \left[1+\displaystyle \frac{{\lambda }_{0}{\nu }_{\mathrm{eff}}\bar{\kappa }}{{\eta }_{\mathrm{eff}}(1-\phi ){\rm{\Lambda }}}+\displaystyle \frac{{\lambda }_{0}}{{\eta }_{\mathrm{eff}}{{\rm{\Lambda }}}^{f}}\right]+{\gamma }_{1}-{\gamma }_{2}+\beta \left[1-\displaystyle \frac{{\lambda }_{0}{\nu }_{\mathrm{eff}}{\kappa }_{\mathrm{eff}}}{\lambda (1-\phi ){\eta }_{\mathrm{eff}}}\right],\end{eqnarray} \tag{ 38 }$$

where we have introduced the dimensionless quantities ${\tilde{v}}^{f}={v}_{\mathrm{ext}}^{f}/{v}_{2}$, ${\tilde{I}}_{\mathrm{ext}}=\lambda {I}_{\mathrm{ext}}(1-\phi )/({\rm{\Lambda }}{\kappa }_{\mathrm{eff}}{v}_{2}),{\gamma }_{1}={\hat{v}}_{{\rm{\Lambda }}}/{v}_{2}$, ${\gamma }_{2}={\lambda }_{0}({\bar{P}}^{H}-{\nu }_{\mathrm{eff}}{I}_{\mathrm{ext}}/{\rm{\Lambda }}+{\nu }_{\mathrm{eff}}\bar{\kappa }{v}_{\mathrm{ext}}^{f}/((1-\phi ){\rm{\Lambda }})+{v}_{\mathrm{ext}}^{f}/{{\rm{\Lambda }}}^{f})/({\eta }_{\mathrm{eff}}{v}_{2})$, where ${\hat{v}}_{{\rm{\Lambda }}}={v}_{{\rm{\Lambda }}}-{v}_{\mathrm{ext}}^{f}-\lambda {I}_{\mathrm{ext}}(1-\phi )/({\rm{\Lambda }}{\kappa }_{\mathrm{eff}})$.

Note that steady states exist only if A > 0, and $-B\pm \sqrt{{B}^{2}-4{AC}}\gt 2A$ or if A < 0 and $-B\,\pm \sqrt{{B}^{2}-4{AC}}\lt 2A$.

4.3.3. Tissue dynamics in the absence of electric currents

We first discuss examples of thickness dynamics in the presence of fluid flow and absence of electric currents, I_ext = 0. Figure 4 shows the growth rate dh/dt as a function of tissue thickness for v₂ > 0, i.e. increased cell division rate at the tissue-substrate interface. There are three possible behaviors: (i) In figure 4(a) the existence of a stable steady state shows that under the chosen conditions the epithelium layer will spontaneously evolve toward a finite thickness irrespective of the initial conditions; with the chosen numbers this thickness is a fraction of the screening length λ₀, (ii) figure 4(b), the existence of a pair of stable–unstable fixed points shows that for initial thicknesses smaller than the unstable one the layer goes spontaneously to the stable value whereas for larger initial values the layer proliferates; figure 4(c), dh/dt being always positive the epithelium layer always proliferates. These different scenarios, depend on the value of the external fluid flow ${v}_{\mathrm{ext}}^{f}$ and tissue homeostatic pressure P^c_h. For a negative homeostatic pressure, when cells die in the bulk, a positive external fluid flow permeating through the tissue gives rise to cell division and hence can counter the effect of cell death. When the two effects exactly balance, the tissue exhibits a steady state with a finite thickness (figure 4(a)). For a high external fluid flow, cell division surpasses cell apoptosis, and the tissue grows indefinitely (figure 4(c)). For intermediate values of ${v}_{\mathrm{ext}}^{f}$, the tissue can exhibit two different behaviors depending on the initial tissue thickness: it can either reach a stable steady state or proliferate out of bounds (figure 4(b)).

Zoom In Zoom Out Reset image size

At steady state, the tissue slab is not homogeneous. The net cell turnover, the cell velocity, the total cell stress, and the fluid pressure exhibit non-trivial profiles along the z-axis. An example of profiles is shown in figure 5 for ${v}_{\mathrm{ext}}^{f}=-3\times {10}^{-10}\,{\rm{m}}\,{{\rm{s}}}^{-1}$ and I_ext = 0 and homeostatic pressure of ${P}_{h}^{c}=-1$ KPa. In this case cells undergo apoptosis in the bulk, and only the cell division in the surface layer prevents the tissue from collapsing. The fluid flowing in the negative z direction, translates into a pressure on the tissue which increases the rate of apoptosis. This effect is cumulative and the pressure is larger at the bottom of the layer than at the top, as can be seen on the stress profile. At the bottom layer, cells divide; as they move up, they balance cell death. As required by boundary conditions the cell velocity vanishes at the free surface. Since the total apoptosis rate is larger in the presence of flow than in the absence of flow, the steady state thickness is smaller in the presence of flow than in its absence.

Zoom In Zoom Out Reset image size

The dynamics of the tissue for v₂ < 0, i.e. cells dying at the tissue-substrate interface, is shown in figure 6. There are three possible scenarios, which are again determined by the values of ${v}_{\mathrm{ext}}^{f}$, and P^c_h. For ${P}_{h}^{c}\gt 0$, and a low ${v}_{\mathrm{ext}}^{f}$, cell death still dominates, and the tissue collapses fully as shown in figure 6(a). For ${P}_{h}^{c}\gt 0$, and intermediate values of ${v}_{\mathrm{ext}}^{f}$, the tissue can either reach a stable thickness or collapse depending on the initial thickness. This possibility is shown in figure 6(b). Eventually, for large enough ${v}_{\mathrm{ext}}^{f}$, the steady state thickness regime disappears in favor of proliferation: depending on the initial thickness, the tissue either collapses or proliferates. This unstable regime is shown in figure 6(c).

Zoom In Zoom Out Reset image size

State diagrams in the ${\tilde{P}}^{H}$–${\tilde{v}}^{f}$ phase space summarize these results in figure 7(a) for v₂ > 0 and figure 7(b) for v₂ < 0. Here ${\tilde{P}}^{H}={\lambda }_{0}{P}_{h}^{c}/({\eta }_{\mathrm{eff}}{v}_{2})$. Changing ${\tilde{P}}^{H}$ can be achieved by tuning P^c_h and, changing ${\tilde{v}}^{f}$ by tuning ${v}_{\mathrm{ext}}^{f}$. As already explained for v₂ > 0, the tissue can exhibit three different behaviors, depending on the values of ${\tilde{P}}^{H}$ and ${\tilde{v}}^{f}$: (i) the green region corresponds to the existence of a stable steady state with a finite thickness irrespective of initial thickness, (ii) the orange region corresponds to tissue evolving toward a finite thickness or proliferating depending on the initial thickness, and, (iii) the pink region corresponds to uncontrolled growth. The solid line delineating the green and the pink regions indicates a continuous transition: the steady state thickness diverges upon approaching the line and reaches infinity on the line. The dotted line delineating the green and orange regions correspond to a spinodal line. It signals the appearance of a stable fixed point at finite thickness simultaneously with the appearance of an unstable fixed point for infinite thickness when crossing the line from green to orange. The dashed line which delineates the orange and pink regions is a line of discontinuous transition, where the pair of stable–unstable fixed points with finite thickness disappears when crossing the line from orange to pink.

Zoom In Zoom Out Reset image size

A negative v₂ implies a sink for cells at the tissue-substrate boundary. There are again three possible scenarios: (i) in the blue region, one has one unstable steady-state: for initial thicknesses smaller than that of the unstable fixed point the tissue collapses, while for larger initial thicknesses the tissue proliferates, (ii) in the brown region, there are two fixed points, one unstable and one stable: the tissue collapses for initial thicknesses smaller than that of the unstable fixed point, and converges to a finite stable thickness if the initial thickness is larger than that of the unstable fixed point, (iii) in the yellow region, there is no fixed point, dh/dt is always negative and the tissue collapses. The continuous line signals a continuous transition from a scenario with one unstable steady state to the total absence of steady states, the value of the unstable thickness going continuously to infinity upon approaching the line from the blue side. The dotted line between the blue and the brown region signals a spinodal transition, with the stable fixed point of the brown region going continuously to infinity upon approaching the blue region. The dashed line signals a discontinuous transition, the stable unstable pair of fixed points disappearing simultaneously for a finite common thickness.

4.3.4. Tissue dynamics in the presence of electric currents

We now discuss the state diagram as a function of the imposed electric current I_ext. Using equation (33) we plot dh/dt versus h in the regime v₂ > 0 in figure 8, for different values of the current I_ext. Like in the previous subsection, we find three possible scenarios: (a) irrespective of its initial thickness, the tissue reaches a stable steady state (b) if the initial thickness is larger than the unstable fixed point value the tissue proliferate and if the initial thickness is smaller than this value the tissue thickness goes to a finite stable value (c) dh/dt is always positive and irrespective of initial conditions the tissue proliferates. In all these cases I_ext > 0, the current flows from the substrate towards the growing interface of the tissue in the positive z direction; this provides an electric field promoting tissue growth if ν₁ > 0. Furthermore for v₂ > 0 the surface layer attached to the substrate acts also as a source of dividing cells. However if the homeostatic pressure, P^c_h, and the external fluid velocity ${v}_{\mathrm{ext}}^{f}$ are negative, both of which favor cell apoptosis, the tissue slab may reach a stable steady state whenever the opposing effects balance exactly. Such a situation is displayed on figure 8(a), where a balance between the positive I_ext, and negative ${v}_{\mathrm{ext}}^{f}$ allows for the tissue to reach a steady state. In figure 8(b), the value of the current I_ext is higher, the electric field can no longer be balanced if the initial thickness is large enough and the tissue reaches an uncontrolled growth state. However if the initial thickness is small enough the effects of the fields can still be counterbalanced, and the tissue can reach a steady state. If the value of I_ext is further increased, as in figure 8(c), the electric field is strong enough to push the tissue to the uncontrolled growth phase, where dh/dt > 0 irrespective of initial conditions.

Zoom In Zoom Out Reset image size

We have plotted steady state profiles for the cell velocity, cell turnover, cell stress, and fluid pressure, in figure 9. The fact that the cell velocity decreases continuously to zero as the ordinate goes to h, shows that cells die everywhere in the bulk. They divide only at the surface. This fact is also clear from the negative value of the turnover rate. The high negative cell stress corresponds to a high pressure which is responsible for the large apoptosis rate.

Zoom In Zoom Out Reset image size

When v₂ < 0, the surface layer attached to the substrate is a cell sink. If the homeostatic pressure is also negative, the cells tend to die in the bulk too. What can prevent the tissue from an immediate collapse is the opposite action of an electric field or of an external fluid flow provided it has the right sign. We find that in the dh/dt versus h plot of figure 10 for v₂ < 0, three different situations can arise, depending on the value of I_ext and ${v}_{\mathrm{ext}}^{f}$: (a) the tissue always collapses with dh/dt < 0, (b) with a pair of unstable–stable fixed points the tissue collapses for initial thicknesses smaller than that of the unstable fixed point, and reaches a stable thickness for initial thickness values larger than the unstable one, (c) with one unstable fixed point, the tissue collapses for an initial thickness smaller than the thickness value of the unstable fixed point, and proliferates for an initial thickness larger than that value. In figure 10(a) we see a total tissue collapse for a low value of electric current ${I}_{\mathrm{ext}}=2\,{\rm{A}}\,{{\rm{m}}}^{-{\rm{2}}}$, which is indicative of the fact that the current and external fluid velocity are not enough to counter the large apoptosis rate due to the combined effect of negative homeostatic pressure and surface apoptosis v₂ < 0. If the current is increased, keeping all other parameters constant, we find parameter values for which the electric field allows to obtain a stable tissue thickness provided the initial thickness is large enough. For a smaller initial thickness the tissue collapses as shown in figure 10(b). If the current is increased further, tissue growth dominates, if the initial thickness is larger than that of the unstable fixed point and drives the tissue to an uncontrolled growth phase. For a smaller initial thickness, we still obtain a tissue collapse. This unstable steady state is shown in figure 10(c).

Zoom In Zoom Out Reset image size

These dynamical states for both positive and negative v₂ can be visualized in two seperate diagrams in the ${\tilde{v}}^{f}$–${\tilde{I}}_{\mathrm{ext}}$ parameter space. We plot these diagrams in figures 11(a) and (b) respectively. With v₂ kept constant, tuning ${\tilde{v}}^{f}$ can be achieved by tuning the external fluid velocity ${v}_{\mathrm{ext}}^{f}$, and tuning ${\tilde{I}}_{\mathrm{ext}}$ can be achieved by tuning the current I_ext flowing through the tissue, with all other parameters kept fixed.

Zoom In Zoom Out Reset image size

Figure 11(a) shows that for positive v₂, we find three scenarios: (a) green region: the tissue reaches a stable thickness irrespective of the initial conditions (b) orange region, existence of a pair of stable–unstable fixed points: for initial thickness values smaller than that of the unstable fixed point, the tissue slab goes to a stable finite thickness whereas for larger initial values the tissue proliferates (c) uncontrolled tissue proliferation (${h}_{{ss}}\to \infty$) represented by the pink region. The solid line indicates a continuous transition between the green and the pink regions: approaching the line from the green side, the steady state thickness increases and diverges on the line. The dotted line between the green and the orange region signals the disappearance of the unstable fixed point of the orange region upon approaching the green region, the corresponding thickness diverging on the line. This is a spinodal line. The dashed line signals the simultaneous disappearance of the fixed point pair of the orange domain upon entering the pink domain. This is a line of discontinuous transition. The point where the three lines meet is again a tricritical point.

Figure 11(b) illustrates the possible scenarios for negative v₂. We again find three possibilities in the ${\tilde{v}}^{f}$–${\tilde{I}}_{\mathrm{ext}}$ space: (a) the blue region corresponds to an unstable steady state, for initial thicknesses smaller than the fixed point value, the tissue collapses whereas for larger initial thicknesses it proliferates, (b) the brown region is defined by the existence of a pair of unstable–stable fixed points, for initial thickness values smaller than that of the unstable fixed point the tissue collapses and for larger initial values it goes to a stable finite thickness corresponding to the stable fixed point (c) in the yellow region, dh/dt is always negative and the tissue collapses. The solid line separating the blue and violet domains is a line of continuous transition: the unstable thickness increases continuously to infinity as the tissue approaches the line from the blue side. The dotted line signals the disappearance of the stable fixed point as the corresponding thickness goes to infinity when the tissue enters the blue region coming from the brown side. This is a spinodal line. The dashed line signals a discontinuous transition as the pair of fixed points disappears when the tissue enters the yellow region coming from the brown side.

Table 1.  Table of parameter values used in the plots.

Parameter list
Parameter	Values
η	${10}^{4}\,\mathrm{Pa}\,\,{\rm{s}}$
$\bar{\eta }$	$2\times {10}^{9}\,\mathrm{Pa}\,{\rm{s}}$
Pch	$-5000\,\mathrm{Pa}$
ϕ	0.99
κ	${10}^{13}\,\mathrm{Pa}\,{\rm{s}}\,{{\rm{m}}}^{-2}$
$\bar{\kappa }$	${10}^{3}\,{\rm{A}}\,{\rm{s}}\,{{\rm{m}}}^{-3}$
ζ	$-1000\,\mathrm{Pa}$
v2	$20\,\mu {\rm{m}}\,{{\rm{d}}}^{-1}$
${\lambda }_{1}$	$-{10}^{8}\,{\rm{N}}\,{{\rm{m}}}^{-3}$
λ	${10}^{5}\,{\rm{N}}\ {{\rm{m}}}^{-2}\,{{\rm{V}}}^{-1}$
Λ1	$3\,{\rm{A}}\,{{\rm{m}}}^{-2}$
Λ	${10}^{-3}\,{\rm{S}}\,{{\rm{m}}}^{-1}$
Λf	${10}^{-11}\,{\mathrm{ms}}^{-1}\,{{\rm{Pa}}}^{-1}$
Jp	${10}^{-10}\,{\rm{m}}\,{{\rm{s}}}^{-1}$
ν	1
ν1	$1\,\mathrm{Pa}\,\,{\rm{m}}\,{{\rm{V}}}^{-1}$
ν2	$1.5\,\mathrm{Pa}\,{\rm{m}}\,{{\rm{V}}}^{-1}$

In this work, we have analyzed the long term growth behavior of planar thick epithelia permeated either by a constant fluid flow or by a constant electric current, or both. The predictions are striking, since one finds that the domain of stability of a finite thickness epithelium is rather small, and that a simple dc electric current or a simple fluid flow is sufficient to either lead to tissue proliferation or to tissue collapse, without any need for genetic mutation. The results should be rather robust, since the above developed arguments are based on symmetry considerations, force conservation laws and cell number balance equations. Yet, this exercice would be futile if the field values required for observing these behaviors were out of experimental reach. Even though our phenomenological theory involves a rather large number of parameters, the formulation can be cast in such a way as to involve only two control parameters, which can be estimated either from values already known experimentally, or from educated guesses. The one feature which is the most difficult to assess is the sign of the coupling parameters. These will need specific experiments to be pinned down. We expect clearly observable effects for flow fields or electric fields somewhat larger than those generated naturally in epithelia, but not orders of magnitude larger. The reason is that either flow or currents comparable to the naturally occurring ones can redistribute proteins in the cells and modify their polarity rather efficiently. The simplest result we obtain is that the steady state thickness of epithelia is proportional to the hydrodynamic screening length introduced in [25] and which one can estimate from [40], with a multiplicative logarithmic correction. The product turns out to be in the milimeter range, a very reasonable feature for stratified squamous epithelia. The main limitation of this theory stems from the assumption of homogeneity in the direction parallel to the tissue layer. It is well possible that instabilities leading to lateral structuration of the tissue exist. This possibility should be investigated in the future. In any case, well controlled experiments are clearly needed, and would give us a deeper insight on the fundamental properties of thick epithelia. The prediction of either collapse or proliferation under suitable conditions is an exciting possibility which should be tested. Last, the process of tissue collapse may take an interesting twist: our analysis does not include explicitly the cell surface layer. It is only included in a flux boundary condition. As a result, the collapse may correspond to a thickness decrease up to the last layer, which usually has very different properties [9, 40, 41], or in other word to a transition from a thick epithelium to a monolayer epithelium. This aspect could be tested experimentally as well.

The authors want to thank Charlie Duclut for reading the manuscript and critical comments. NS would also like to thank Marko Popovic, and Keisuke Ishihara for insightful comments. JP wants to thank Thuan Beng Saw for fruitful discussions.

One can obtain an order of magnitude of the coefficient λ = λ₂ + 2λ₃/3, by estimating the shear stress created by the electo-osmotic flow in the intercellular region.

Zoom In Zoom Out Reset image size

The hydrodynamic stress is as usual given by $\sigma =\eta {\partial }_{x}v$, where η is the interstitial fluid viscosity, and ${\partial }_{x}v$ the shear rate due to electro-osmosis at the cell membrane surface. Standard calculations [42] yield

$$\begin{eqnarray}&&{\partial }_{x}v\simeq \displaystyle \frac{\epsilon {\epsilon }_{0}{\zeta }_{E}}{{\lambda }_{{\rm{D}}}\eta }{E}_{z}^{\mathrm{cleft}},\end{eqnarray} \tag{ 39 }$$

where is the permittivity of the fluid (here water), ₀ is the permittivity of vacuum, ${\zeta }_{E}$ is the zeta potential of the cell membrane, ${E}_{z}^{\mathrm{cleft}}$ is the electric field in the cleft, and λ_D is the Debye screening length.

Then the stress is expressed in terms of ${E}_{z}^{\mathrm{cleft}}$ as (see figure A1)

$$\begin{eqnarray}&&\sigma \simeq \displaystyle \frac{\epsilon {\epsilon }_{0}{\zeta }_{E}}{{\lambda }_{{\rm{D}}}}{E}_{z}^{\mathrm{cleft}},\end{eqnarray} \tag{ 40 }$$

which gives the total force on the cell to be $4{{ll}}_{c}\sigma$, with l being the height of the cell, and l_c being the width of the cell. The force per unit volume can then be written as

$$\begin{eqnarray}&&{f}_{E}\simeq \displaystyle \frac{4{{ll}}_{c}\sigma }{{{ll}}_{c}^{2}}=\displaystyle \frac{4\epsilon {\epsilon }_{0}{\zeta }_{E}}{{\lambda }_{{\rm{D}}}{l}_{c}}{E}_{z}^{\mathrm{cleft}}.\end{eqnarray} \tag{ 41 }$$

The relation between the coarse grained electric field E_z, and the actual field ${E}_{z}^{\mathrm{cleft}}$ in the intercellular domain is simply given by the geometrical relation ${E}_{z}^{\mathrm{cleft}}\simeq {l}_{c}{E}_{z}/(2W)$ stemming from current conservation ignoring the current passing through the cells, where W is the width of the intercellular cleft. Using this relation, the force per unit volume can be expressed as

$$\begin{eqnarray}&&{f}_{E}\simeq \displaystyle \frac{2\epsilon {\epsilon }_{0}{\zeta }_{E}}{{\lambda }_{{\rm{D}}}W}{E}_{z},\end{eqnarray} \tag{ 42 }$$

which leads to

$$\begin{eqnarray}&&\lambda \simeq \displaystyle \frac{2}{W}\displaystyle \frac{\epsilon {\epsilon }_{0}{\zeta }_{E}}{{\lambda }_{{\rm{D}}}}.\end{eqnarray} \tag{ 43 }$$

Using ${\epsilon }_{0}=8.85\times {10}^{-12}\,{\rm{F}}\,{{\rm{m}}}^{-1},\epsilon \simeq 80$ for water, ${\zeta }_{E}\simeq 3\times {10}^{-2}\,\mathrm{mV}$, W ≃ 10⁻⁷ m, and ${\lambda }_{{\rm{D}}}\simeq {10}^{-9}\,{\rm{m}}$, we obtain $\lambda \simeq 3\times {10}^{5}\,{\rm{N}}\ {{\rm{m}}}^{-2}\,{{\rm{V}}}^{-1}$.

The term λ₁ in (17) corresponding to cell polarity has also never been measured to our knowledge. If we assume that it is comparable to the electric term, under conditions such that the average flux vanishes, we can infer an estimate based on the knowledge of typical potential differences generated by epithelial cells i.e. mV. Then λ₁ ≃ −λE_z. As the potential difference $V\simeq {10}^{-3}\,{\rm{V}}$ in MDCK blisters [43, 44] and width of a cell l_c ≃ 10⁻⁵ m, then ${E}_{z}\simeq V/{l}_{c}\simeq 0\times {10}^{2}\,{\rm{V}}\,{{\rm{m}}}^{-1}$. Thus $\lambda {E}_{z}\simeq -\lambda {E}_{z}\simeq -3\times {10}^{5}{10}^{-3}/{10}^{-5}\,{\rm{N}}\,{{\rm{m}}}^{-3}\,\simeq -{10}^{8}\,{\rm{N}}\,{{\rm{m}}}^{-3}$

We next estimate the values of $\nu ,{\nu }_{1}$, and ν₂ in equations (15) and (16). We argue that ν should be of O(1), as it is a dimensionless parameter as evident from (15). To determine the value of ν₁, we assume that ${\nu }_{1}{p}_{\alpha }{E}_{z}$ plays the role of a pressure, even though the underlying physics might be more subtle. Thus we equate ${\nu }_{1}{p}_{\alpha }{E}_{z}=\lambda {{hE}}_{z}$, so that ${\nu }_{1}=\lambda h\sim O(1)\,\mathrm{Pa}\,{\rm{m}}\,{{\rm{V}}}^{-1}$. For ν₂, we compare equation (16) with (40), from which we can express ${\nu }_{2}\simeq \epsilon {\epsilon }_{0}{\zeta }_{E}{l}_{c}/(2W{\lambda }_{{\rm{D}}})=80\times 8.85\times {10}^{-11}\times 3\times {10}^{-2}/{10}^{-9}\times 50\,\mathrm{Pa}\,\,{\rm{m}}\,{{\rm{V}}}^{-1}\simeq 1.5\,\mathrm{Pa}\,{\rm{m}}\,{{\rm{V}}}^{-1}$.

The coefficient $\bar{\kappa }$ may be estimated using standard streaming potential relations [45] adapted to the intercellular domains, similar in spirit to what was done for estimating λ. We obtain

$$\begin{eqnarray}&&\bar{\kappa }\simeq \displaystyle \frac{12\epsilon {\zeta }_{E}}{{{Wl}}_{c}}.\end{eqnarray} \tag{ 44 }$$

Using the values above and ${l}_{c}\simeq {10}^{-5}\,{\rm{m}}$, gives $\bar{\kappa }\simeq {10}^{3}\,{\rm{A}}\,{\rm{s}}\,{{\rm{m}}}^{-3}$.

The coefficient Λ^f can also be estimated using standard fluid flux calculated in the intercellular cleft and from geometric relations, just like we used in the calculation of λ, we obtain the expression

$$\begin{eqnarray}&&{{\rm{\Lambda }}}^{f}\simeq \displaystyle \frac{{W}^{3}}{6{l}_{c}^{2}\eta }.\end{eqnarray} \tag{ 45 }$$

Using the values of $W,{l}_{c}$, and η, given above, yields ${{\rm{\Lambda }}}^{f}\simeq {10}^{-11}\,{\rm{m}}\,{\mathrm{Pa}}^{-1}\,{{\rm{s}}}^{-1}$.

To estimate the value of J_p we note that at steady state the pressure difference ${P}_{\mathrm{ext}}^{1}\mbox{--}{P}^{f}$ is equal to the osmotic pressure difference ${{\rm{\Pi }}}_{\mathrm{ext},0}^{1}\mbox{--}{{\rm{\Pi }}}_{\mathrm{int},0}^{1}$, where ${J}_{p}={{\rm{\Lambda }}}^{f}({{\rm{\Pi }}}_{\mathrm{ext},0}^{1}\mbox{--}{{\rm{\Pi }}}_{\mathrm{int},0}^{1})$. Furthermore we know that the curvature radius of a cell membrane R should be much larger than the cell thickness l_c, for the cell to be stable. So the pressure difference ${P}_{\mathrm{ext}}^{1}\mbox{--}{P}^{f}\simeq \gamma /R$ should be much smaller than γ/l_c. This ensures, ${{\rm{\Pi }}}_{\mathrm{ext},0}^{1}\mbox{--}{{\rm{\Pi }}}_{\mathrm{int},0}^{1}\leqslant \gamma /{l}_{c}$. Using the values of γ ≃ 10⁻⁴ N m⁻¹, and l_c ≃ 10⁻⁵ m, we obtain ${J}_{p}\simeq {10}^{-10}\,{\rm{m}}\,{{\rm{s}}}^{-1}$.