Reduced-Order Model for Cell Volume Homeostasis: Application to Aqueous Humor Production

1. Introduction

The mechanical integrity of a cell is essential to its function and the tissues and organ systems that it supports [1,2]. Cell volume stabilization is a fundamental function to preserve the mechanical integrity of a cell [3]. Cell volume stabilization along time is based on the balance among forces acting on the cell membrane which emanate from electrochemical and fluid-mechanical mechanisms [4]. The experimental characterization of force equilibrium at the nanoscale level is a highly complex task because of the difficulty of performing measurements in vivo (see [5]), and the presence of two regimes of mechanical behavior in living cells, at short and long time scales (see [6]).

Based on these considerations, we address in this article the theoretical study of cell volume dynamics and apply the formulation to the process of production of aqueous humor (AH) in the eye. AH is a slowly moving fluid that is continuously produced by the ciliary body of the eye, whose main function is to keep the anterior segment of the eye clean from metabolic wastes and to preserve the spherical shape of the eyeball by establishing the intraocular pressure (IOP) (see [7,8,9,10]).

Keeping IOP within normal levels (10 to 21 mmHg) is central to prevent the occurrence and progression of neurodegenerative diseases of the retina such as glaucoma [11,12,13]. This opens up the question to understand the biophysical mechanisms that determine the value of IOP and, consequently, the causes of the pathology, and to devise effective therapies for its cure. In this respect, the development of a mathematical model and of a computational virtual laboratory (CVL) for the simulation of AH production has been object of investigation in recent years (see [14,15,16,17]).

In this article, we propose a theoretical formulation based on homogeneous mixtures including neutral and charged solutes (see [18], Chapter 13) and utilizing a model reduction procedure from three spatial dimensions to zero spatial dimensions. The model resulting from reduction is a coupled system of highly stiff nonlinear ordinary differential equations (ODE) which are solved using the $\theta$-method and the Matlab function ode15s. Simulation tests are run to characterize the values of model parameters at baseline conditions in AH production. The model is subsequently used to investigate the relative importance of (a) impermeant charged proteins; (b) sodium-potassium pump; and (c) carbonic anhydrase in the AH production process. Results suggest that (a) and (b) play a role whereas (c) does not have a significant weight, at least for low CO₂ values.

2. Materials and Methods

In Section 2.1 we describe the geometrical representation of the cell volume and surface; in Section 2.2 we describe the heterogeneous structure of the cell membrane; in Section 2.3 we provide a conceptual scheme of transport of water and solute across the cell membrane. In Section 2.4 we introduce a continuum model of the cell based on the theory of mixtures, to describe the spatial and temporal response of the cell volume to the concomitant action of fluid and electrochemical forces. In Section 2.5 we derive a reduced-order model of the cell volume that comprises ordinary differential and algebraic equations. Finally, in Section 2.6 and Section 2.7 we describe the compact form of the reduced-order model and its numerical approximation, respectively. In Appendices Appendix A and Appendix B, we provide the expressions of the fluid velocity, molar flux densities and net production rates that are involved in cellular metabolism.

2.1. Geometrical Description of the Cell

Let $t\ge 0$ denote the time variable (units: s). We denote by ${\Omega }_t:=\Omega \left(t\right)$ the three-dimensional (3D) body representing the cell at time t and ${\mathcal{V}}_t:=\mathcal{V}\left(t\right)$ the volume of the cell at time t, such that ${\mathcal{V}}_0=4\pi R_c^3/3$ is the value of the cell volume at time $t=0$, $R_c$ being the cell radius in the initial (undeformed) configuration (units: m). From the point of view of Continuum Mechanics the body is a deformable, electrically charged, homogeneous mixture comprising a fluid constituent, $N_{\alpha }$ moving charged solute constituents (ions) and $N_{\beta }$ moving neutral solute constituents (see [18], Chapter 13). The cell volume contains also impermeant charged proteins whose chemical valency and molar density is such that electroneutrality holds in ionic homeostasis conditions (see [3] and [19]). The boundary of ${\Omega }_t$ is denoted by $\partial {\Omega }_t$ and $\mathbf{n}$ is the outward unit normal vector on $\partial {\Omega }_t$. The surface of the cell at time t is denoted by ${\mathcal{S}}_t:=\mathcal{S}\left(t\right)$, such that ${\mathcal{S}}_0=4\pi R_c^2$ is the value of the cell surface at time $t=0$. In the remainder of the article, we denote by $\mathbf{x}$ the spatial coordinate vector of any point of the cell with respect to a fixed system of reference. We also refer to the interior of ${\Omega }_t$ as the intracellular region (shortly, $in$), to the exterior of ${\Omega }_t$ as the extracellular region (shortly, $out$) and to the surface of ${\Omega }_t$ as the membrane (shortly, m).

Assumption 1.

We assume that cell deformation occurs only in the radial direction and that it does not depend on the position on the cell surface.

Assumption 2.

We assume that the electric potential, solute concentrations and fluid pressure in the extracellular region are given functions of space and time.

Figure 1 illustrates in solid line and dashed line the initial and deformed configurations of the cell, respectively. In agreement with Assumption 1, both configurations are spherical.

2.2. Porous and Lipid Structure of the Cell Membrane

The intra and extracellular regions are separated by a 3D membrane whose thickness $t_m$ is such that the ratio $\delta :=t_m/R_c$ is $\ll 1$. In the limit $\delta \to 0$, the geometrical representation of the 3D membrane degenerates into the two-dimensional (2D) manifold $\partial {\Omega }_t$. The membrane surface can be represented as an heterogeneous porous medium whose solid constituent is lipid material. The lipid surface is impermeable to water and solutes whereas the pore surface allows the exchange of fluid and solutes between intra and extracellular regions. The pore surface contains a variety of membrane proteins, denoted by $mp$. Each membrane protein $mp$ is characterized by a dimensionless function ${\Phi }^{mp}$, henceforth referred to as surface fraction, defined as

$$\begin{array}{ccc}& {\Phi }^{mp}\left(t\right):={\displaystyle \frac{{\mathcal{S}}^{mp}\left(t\right)}{\mathcal{S}\left(t\right)}}\mathcal{S}\left(t\right)>0,t\ge 0,\hfill & \end{array}$$

(1a)

where ${\mathcal{S}}^{mp}\left(t\right)$ is the total area occupied by the protein $mp$ at time t.

Assumption 3.

Let ${\mathcal{S}}_t>0$. We assume that

$$\begin{array}{ccc}& {\displaystyle \frac{d{\Phi }^{mp}\left(t\right)}{dt}}=0t>0.\hfill & \end{array}$$

(1b)

Equation (1b) implies that

$$\begin{array}{ccc}& {\Phi }^{mp}\left(t\right)={\Phi }^{mp}\left(0\right)={\displaystyle \frac{{\mathcal{S}}^{mp}\left(0\right)}{\mathcal{S}\left(0\right)}}t\ge 0.\hfill & \end{array}$$

(1c)

The first type of membrane protein that we introduce in this article is the aquaporin (AQP), which is a specialized protein for rapid transmembrane water exchange between intracellular and extracellular regions (see [20,21,22]). The quantity ${\Phi }^{AQP}$ represents the AQP surface fraction. Other membrane proteins on the cell surface comprise carrier proteins, which permit neutral solute exchange through the cell membrane by the mechanism of facilitated diffusion (see [23] and [24]), and ion channels, ion exchangers and ion pumps, which permit charged solute exchange through the cell membrane (see [25]). The total surface fractions for each considered membrane protein are defined as:

$${\Phi }^{carr}=\sum _{\beta \in S_{\beta }}{\Phi }_{\beta }^{carr},$$

(2a)

$${\Phi }^{ch}=\sum _{\alpha \in S_{\alpha }}{\Phi }_{\alpha }^{ch},$$

(2b)

$${\Phi }^{exch}=\sum _{\alpha \in S_{\alpha }}{\Phi }_{\alpha }^{exch},$$

(2c)

$${\Phi }^{pump}=\sum _{\alpha \in S_{\alpha }}{\Phi }_{\alpha }^{pump},$$

(2d)

in such a way that the total pore surface fraction is

$$\begin{array}{ccc}& {\Phi }_{tot}^p={\Phi }^{carr}+{\Phi }^{ch}+{\Phi }^{exch}+{\Phi }^{pump}+{\Phi }^{AQP}.\hfill & \end{array}$$

(3)

The remainder of the cell surface is occupied by the lipid constituent of the membrane, whose surface fraction is

$$\begin{array}{ccc}& {\Phi }^{lip}=1-{\Phi }_{tot}^p.\hfill & \end{array}$$

(4)

Assumption (3) implies that the membrane protein surface fractions in (2) do not depend on time, and that also the lipid surface fraction does not depend on time because of (4).

2.3. Water and Solute Transport Across the Cell Membrane

In this section we provide a conceptual scheme of transport of water and solute across the cell membrane and we refer to [21] and [26] for more details about the subject.

Figure 2 illustrates a schematic representation of a zoomed view of the cell membrane and of the intra and extracellular regions. The scheme shows water molecules and solutes (neutral and charged) in motion across the membrane.

Assumption 4.

Based on the scheme in Figure 2 we make the following assumptions:

A1. 

Water is transported across the membrane, in a selective manner, via aquaporins;

A2. 

Water and solutes are co-transported across the membrane via ionic channels and carrier proteins;

A3. 

Water velocity inside ionic channels and carrier proteins is the same as water velocity inside aquaporins.

2.4. Continuum-Based Model of the Cell

This section is structured in four subsections devoted to fluid motion, neutral and charged solutes, and electric potential.

2.4.1. Fluid Motion

Assumption 5.

We assume that water flow is slow and inertial forces can be neglected with respect to the viscous effects.

According to Assumption 5, the motion of water flow across the cell membrane can be described by the linear Stokes equation system (see [18], Section 10.5.1):

$$\nabla ·\mathbf{v}=0,$$

(5a)

$$\nabla ·\mathbb{T}+\mathbf{b}=\mathbf{0},$$

(5b)

$$\mathbb{T}=-p\mathbb{I}+2{\mu }_f{\nabla }_S\mathbf{v},$$

(5c)

where $\mathbf{v}$ (units: $\mathrm{m}{\mathrm{s}}^{-1}$) and p (units: $\mathrm{N}{\mathrm{m}}^{-2}$) are the fluid velocity and pressure inside the aquaporin, respectively, $\mathbb{T}$ is the fluid stress tensor (units: $\mathrm{N}{\mathrm{m}}^{-2}$), $\mathbf{b}$ is the force density acting on the fluid (units: $\mathrm{N}{\mathrm{m}}^{-3}$), ${\mu }_f$ is the fluid dynamic viscosity (units: $\mathrm{N}{\mathrm{m}}^{-2}\mathrm{s}=\mathrm{P}\mathrm{a}\mathrm{s}$) and ${\nabla }_S$ is the symmetric gradient operator.

Equation (5a) expresses mass conservation in local form; equation (5b) expresses linear momentum balance in local form; equation (5c) is the constitutive law for a Newtonian fluid.

2.4.2. Neutral solutes motion

We denote by $S_{\beta }$ the set of neutral solutes, with $\mathrm{card}\left(S_{\beta }\right)=N_{\beta }$. Each solute $\beta \in S_{\beta }$ is described by its number density $n_{\beta }$ and molar density $c_{\beta }$ (units: $\mathrm{mol}{\mathrm{m}}^{-3}=\mathrm{mM}$). The advection-diffusion model is adopted to represent the motion of neutral solutes across the cell membrane (see [18], Chapter 12):

$${\displaystyle \frac{\partial c_{\beta }}{\partial t}}+\nabla ·{\mathbf{j}}_{\beta }=P_{\beta }-C_{\beta }\beta =1,\dots ,N_{\beta },$$

(6a)

$${\mathbf{j}}_{\beta }=c_{\beta }\mathbf{v}-D_{\beta }\nabla c_{\beta }\beta =1,\dots ,N_{\beta }.$$

(6b)

Equations (6a) express mass balance in local form for each solute $\beta \in S_{\beta }$. The vector ${\mathbf{j}}_{\beta }$ is the molar flux density of solute $\beta$ (units: $\mathrm{mol}{\mathrm{m}}^{-2}{\mathrm{s}}^{-1}=\mathrm{mM}\mathrm{m}{\mathrm{s}}^{-1}$) accounting for a passive advective contribution due to water motion inside the carrier protein channel (see Assumption 4-A2), and a diffusive contribution expressed by Fick’s law in which $D_{\beta }$ is the diffusion coefficient of solute $\beta$ (units: ${\mathrm{m}}^2{\mathrm{s}}^{-1}$). The scalar functions $P_{\beta }$ and $C_{\beta }$ represent production and consumption rates of $c_{\beta }$ (units: $\mathrm{mol}{\mathrm{m}}^{-3}{\mathrm{s}}^{-1}=\mathrm{mM}{\mathrm{s}}^{-1}$).

2.4.3. Charged solutes motion

We denote by $S_{\alpha }$ the set of charged solutes (ions) with $\mathrm{card}\left(S_{\alpha }\right)=N_{\alpha }$. Each ion $\alpha \in S_{\alpha }$ is described by its number density $n_{\alpha }$, molar density $c_{\alpha }$, charge number $z_{\alpha }$, diffusion coefficient $D_{\alpha }$ and electric mobility ${\mu }_{\alpha }^{el}$ (units: ${\mathrm{m}}^2{\mathrm{V}}^{-1}{\mathrm{s}}^{-1}$). The number and molar densities are related by the equation $n_{\alpha }=N_{\mathrm{Av}}{c}_{\alpha }$ (units: ${\mathrm{m}}^{-3}$), where $N_{\mathrm{Av}}=6.02214076·{10}^{23}{\mathrm{mol}}^{-1}$ is Avogadro’s constant. The diffusion coefficient and the electric mobility are proportional through Einstein’s relation

$$\begin{array}{ccc}& D_{\alpha }={\displaystyle \frac{{\mu }_{\alpha }^{el}{V}_{th}}{|z_{\alpha }|}}\alpha =1,\dots ,N_{\alpha },\hfill & \end{array}$$

(7)

where $V_{th}=\left(K_{B}T\right)/q$ is the thermal voltage (units: V), $K_B$, T and q denoting Boltzmann’s constant, absolute temperature and electron charge, respectively. The Nernst-Planck (NP) model is adopted to represent ion transport and exchange across the cell membrane (see [18], Chapter 13):

$${\displaystyle \frac{\partial c_{\alpha }}{\partial t}}+\nabla ·{\mathbf{j}}_{\alpha }=P_{\alpha }-C_{\alpha }\alpha =1,\dots ,N_{\alpha },$$

(8a)

$${\mathbf{j}}_{\alpha }={\mathbf{j}}_{\alpha }^p+{\mathbf{j}}_{\alpha }^a\alpha =1,\dots ,N_{\alpha },$$

(8b)

$${\mathbf{j}}_{\alpha }^p={\mathbf{j}}_{\alpha }^{p,edw}+{\mathbf{j}}_{\alpha }^{p,exch}\alpha =1,\dots ,N_{\alpha },$$

(8c)

$${\mathbf{j}}_{\alpha }^{p,edw}=c_{\alpha }{\mathbf{v}}_{\alpha }-D_{\alpha }\nabla c_{\alpha }\alpha =1,\dots ,N_{\alpha },$$

(8d)

$${\mathbf{v}}_{\alpha }=\mathbf{v}+{\mu }_{\alpha }^{el}{\displaystyle \frac{|z_{\alpha }|}{z_{\alpha }}}\mathbf{E}.$$

(8e)

Equations (8a) express mass balance in local form for each ion $\alpha \in S_{\alpha }$. The vector ${\mathbf{j}}_{\alpha }$ is the molar flux density of ion $\alpha$, accounting for a passive ion molar flux density ${\mathbf{j}}_{\alpha }^p$ and an active ion molar flux density ${\mathbf{j}}_{\alpha }^a$. The vector ${\mathbf{j}}_{\alpha }^{p,edw}$ is the velocity-extended NP equation for ion electrodiffusion and accounts for a passive advective contribution and a diffusive contribution expressed by Fick’s law (see [18], Chapter 13). The vector field ${\mathbf{v}}_{\alpha }$ is the generalized drift velocity of ion $\alpha$ due to water motion inside the ion channel (see Assumption 4-A2) and Faraday’s law under quasi-static approximation (see [18], Chapter 14)

$$\begin{array}{ccc}& \mathbf{E}=-\nabla \psi ,\hfill & \end{array}$$

(8f)

where $\mathbf{E}$ and $\psi$ are the electric field (units: $\mathrm{V}{\mathrm{m}}^{-1}$) and electric potential (units: V), respectively.

The vector ${\mathbf{j}}_{\alpha }^{p,exch}$ expresses passive ion exchange mediated by transporters and co-transporters. The scalar functions $P_{\alpha }$ and $C_{\alpha }$ represent production and consumption rates of $c_{\alpha }$.

2.4.4. Electric potential

The spatial and temporal distribution of the electric potential $\psi$ across the cell membrane is described by the following differential system:

$$\nabla ·{\mathbf{J}}_{tot}=0,$$

(9a)

$${\mathbf{J}}_{tot}={\mathbf{J}}^{disp}+{\mathbf{J}}_{tot}^{cond},$$

(9b)

$${\mathbf{J}}^{disp}={\displaystyle \frac{\partial }{\partial t}}\left({\epsilon }_m\mathbf{E}\right)={\displaystyle \frac{\partial }{\partial t}}\left(-{\epsilon }_m\nabla \psi \right),$$

(9c)

$${\mathbf{J}}_{tot}^{cond}=F\sum _{\alpha \in S_{\alpha }}{z}_{\alpha }{\mathbf{j}}_{\alpha },$$

(9d)

where ${\epsilon }_m$ is the dielectric permittivity of the cell membrane lipid layer (units: $\mathrm{F}{\mathrm{m}}^{-1}$) and F is Faraday’s constant (units: $\mathrm{C}{\mathrm{mol}}^{-1}$). Equation (9a) is a consequence of Maxwell’s equations (see [27] and [18], Chapter 4) and expresses the continuity of the total electric current density ${\mathbf{J}}_{tot}$ (units: $\mathrm{C}{\mathrm{m}}^{-2}{\mathrm{s}}^{-1}=\mathrm{A}{\mathrm{m}}^{-2}$), given by the sum of the total conduction current density ${\mathbf{J}}_{tot}^{cond}$ and the displacement current density ${\mathbf{J}}^{disp}$, in which the electric field $\mathbf{E}$ is expressed as a function of the electric potential $\psi$ via equation (8f).

2.5. Cell Reduced-Order Model

This section is structured in five subsections devoted to the derivation of the reduced-order model for cell surface normal velocity, cell volume, neutral and charged solute molar densities and electric potential.

2.5.1. Time evolution of normal velocity across a single AQP

We make the following assumption on the geometry of the AQP shown in Figure 2.

Assumption 6.

The AQP is geometrically represented as a cylinder ${\omega }_p$ with radius $r_p$ and axial length $t_m$ as depicted in Figure 3.

Moreover,

Assumption 7.

We assume that:

1. 

the fluid velocity has only the axial component $V_s$;

2. 

$V_s=v_{p,n}\left(t\right)\eta \left(r\right)$, with $t\ge 0$ and $r\in [0,r_p]$;

3. 

$\eta \left(r\right)=2\left[1-{(r/r_p)}^2\right]$ (Poiseuille flow);

4. 

the force density has only the axial component $b_s$;

5. 

$b_s=b_{p,s}(s,t)$, with $s\in [0,t_m]$ and $t\ge 0$.

Inserting Assumptions A7 into equations (5) we obtain the following expression of the normal fluid velocity inside the pore channel of the AQP

$$\begin{array}{ccc}& v_{p,n}\left(t\right)={\displaystyle \frac{r_p^2}{8{\mu }_f}}\left[-{\displaystyle \frac{\partial p(s,t)}{\partial s}}+b_{p,s}(s,t)\right]s\in [0,t_m],t\ge 0.\hfill & \end{array}$$

(10a)

Integrating (10a) we get

$$\begin{array}{ccc}& p(s,t)=p(0,t)-{\displaystyle \frac{8{\mu }_f{v}_{p,n}\left(t\right)}{r_p^2}}s+{\int }_0^s{b}_{p,s}(\xi ,t)d\xi s\in [0,t_m],t\ge 0,\hfill & \end{array}$$

(10b)

For any function $U=U\left(t\right)$, we define the transmembrane difference operator

$$\begin{array}{ccc}& \Delta U\left(t\right):=U_{in}\left(t\right)-U_{out}\left(t\right)t\ge 0.\hfill & \end{array}$$

(10c)

Setting $s=t_m$ in (10b) and applying the definition (10c) to the variables p and $\Pi$, we obtain

$$\begin{array}{ccc}& v_{p,n}\left(t\right)={\mathcal{L}}_p\left[\Delta p\left(t\right)-\Delta \Pi \left(t\right)\right]t\ge 0,\hfill & \end{array}$$

(10d)

where

$$\begin{array}{ccc}& {\mathcal{L}}_p:={\displaystyle \frac{r_p^2}{8{\mu }_f{t}_m}}\hfill & \end{array}$$

(10e)

is the hydraulic conductivity of the membrane (units: $\mathrm{m}{\left(\mathrm{P}\mathrm{a}\mathrm{s}\right)}^{-1}$) and

$$\begin{array}{ccc}& \Pi \left(t\right):=\int b_{p,s}(\xi ,t)d\xi +C\hfill & \end{array}$$

(10f)

is the total osmo-oncotic pressure, C being an arbitrary integration constant. Relation (10d) is Starling’s equation [28] and represents the mathematical expression of the normal fluid velocity across the single aquaporin.

2.5.2. Time Evolution of Cell Surface Normal Velocity

In order to determine the motion of the whole cell surface, we need to define an average normal fluid velocity to represent the collective contribution of the AQPs. To this purpose, we denote by ${\sigma }_{AQP}$ the surface density of aquaporins, defined as the number $N_{AQP}$ of AQPs per square meter, and by ${\mathcal{S}}_{AQP}^{tot}\left(t\right)$ the total surface area occupied by the AQPs. Using Assumption 3, we obtain

$$\begin{array}{ccc}& {\Phi }^{AQP}={\displaystyle \frac{{\mathcal{S}}_{AQP}^{tot}\left(0\right)}{{\mathcal{S}}_0}}={\displaystyle \frac{{\sigma }_{AQP}{\mathcal{S}}_0\pi r_p^2}{{\mathcal{S}}_0}}={\sigma }_{AQP}\pi r_p^2.\hfill & \end{array}$$

(11a)

Let us introduce the total water volumetric flow rate across the cell surface (units: ${\mathrm{m}}^3{\mathrm{s}}^{-1}$)

$$\begin{array}{ccc}& Q_w\left(t\right):={\int }_{\partial {\Omega }_t}\mathbf{v}(\mathbf{x},t)·\mathbf{n}d(\partial {\Omega }_t)=v_{p,n}\left(t\right){\mathcal{S}}_{AQP}^{tot}\left(t\right)t\ge 0,\hfill & \end{array}$$

(11b)

and the average normal fluid velocity across the cell surface

$$\begin{array}{ccc}& \langle \mathbf{v}(\mathbf{x},t)·\mathbf{n}\rangle \left(t\right):={\displaystyle \frac{Q_w\left(t\right)}{{\mathcal{S}}_t}}={\displaystyle \frac{{\mathcal{S}}_{AQP}^{tot}\left(t\right)}{{\mathcal{S}}_t}}{v}_{p,n}\left(t\right)={\Phi }^{AQP}{v}_{p,n}\left(t\right)t\ge 0,\hfill & \end{array}$$

(11c)

where ${\Phi }^{AQP}$ is given by (11a) and $v_{p,n}$ is given by (10d).

The pore fluid velocity $v_{p,n}$ is typically quite large (even thousands of $\mu \mathrm{m}{\mathrm{s}}^{-1}$, see [29]) to favor a fast transmembrane exchange of water molecules. The average normal fluid velocity $\langle \mathbf{v}·\mathbf{n}\rangle$, instead, is considerably smaller (less that $1\mu \mathrm{m}{\mathrm{s}}^{-1}$, see [30]). For this reason, in the remainder of the article, we introduce the following definition.

Definition 1.

Let us denote by $v_{cell,n}$ the cell surface normal velocity. We set

$$\begin{array}{ccc}& v_{cell,n}\left(t\right):=\langle \mathbf{v}(\mathbf{x},t)·\mathbf{n}\rangle ={\Phi }^{AQP}{v}_{p,n}\left(t\right)\mathbf{x}\in \partial {\Omega }_t,t\ge 0,\hfill & \end{array}$$

(12)

where $v_{p,n}$ is defined in (10d).

2.5.3. Time evolution of cell volume

The time evolution of ${\mathcal{V}}_t$ is described by the cell mass balance equation in integral form (see [18], Chapter 6)

$$\begin{array}{ccc}& {\displaystyle \frac{dM\left(t\right)}{dt}}=-{\int }_{\partial {\Omega }_t}{\rho }_w\mathbf{v}(\mathbf{x},t)·\mathbf{n}d(\partial {\Omega }_t)+{\int }_{{\Omega }_t}{R}_w(\mathbf{x},t)d{\Omega }_{t}t\ge 0,\hfill & \end{array}$$

(13)

where ${\rho }_w$ is the mass density of water (units: $\mathrm{K}\mathrm{g}{\mathrm{m}}^{-3}$), $M\left(t\right)={\int }_{{\Omega }_t}{\rho }_{w}d{\Omega }_t$ is the mass of the cell at time t (units: $\mathrm{Kg}$), and

$$\begin{array}{ccc}& R_w(\mathbf{x},t):=P_w(\mathbf{x},t)-C_w(\mathbf{x},t)\mathbf{x}\in {\Omega }_t,t\ge 0\hfill & \end{array}$$

(14)

is the intracellular water mass density net production rate, $P_w$ and $C_w$ being water mass density production and consumption rates, respectively (units: $\mathrm{Kg}{\mathrm{m}}^{-3}{\mathrm{s}}^{-1}$).

Assumption 8.

We assume that:

$$P_w(\mathbf{x},t)={\mathcal{P}}_w\left(t\right)\Phi \left(\mathbf{x}\right)\mathbf{x}\in {\Omega }_t,t\ge 0,$$

(15a)

$$C_w(\mathbf{x},t)={\mathcal{C}}_w\left(t\right)\Phi \left(\mathbf{x}\right)\mathbf{x}\in {\Omega }_t,t\ge 0,$$

(15b)

where ${\mathcal{P}}_w$ and ${\mathcal{C}}_w$ represent the time-dependent intracellular water mass density production and consumption rates, and the shape function Φ is such that

$$\begin{array}{ccc}& {\int }_{{\Omega }_t}\Phi \left(\mathbf{x}\right)d{\Omega }_t={\mathcal{V}}_t.\hfill & \end{array}$$

(15c)

Using Tr.1 in Assumptions A4, using Definitions 1 and A8 into the right-hand side of (13), and noting that $\mathcal{S}\left(t\right)=\gamma {\left(\mathcal{V}\left(t\right)\right)}^{2/3}$, with $\gamma :={\left(36\pi \right)}^{1/3}$, we obtain the following reduced-order model of cell volume motion:

$$\begin{array}{ccc}& {\displaystyle \frac{d\mathcal{V}\left(t\right)}{dt}}=-\gamma j_{v,n}\left(t\right){\left(\mathcal{V}\left(t\right)\right)}^{2/3}+{\mathcal{R}}_w\left(t\right)\mathcal{V}\left(t\right)t>0,\hfill & \end{array}$$

(16a)

$$\begin{array}{ccc}& j_{v,n}\left(t\right):=v_{cell,n}\left(t\right)={\Phi }^{AQP}{v}_{p,n}\left(t\right)t>0,\hfill & \end{array}$$

(16b)

$$\begin{array}{ccc}& \mathcal{V}\left(0\right)={\mathcal{V}}_0,\hfill & \end{array}$$

(16c)

where $v_{p,n}$ is given by (10d) and

$$\begin{array}{ccc}& {\mathcal{R}}_w\left(t\right):=k_{w,prod}\left(t\right)-k_{w,cons}\left(t\right)t\ge 0\hfill & \end{array}$$

(16d)

is the water volume net production rate (units: ${\mathrm{s}}^{-1}$), having set $k_{w,prod}\left(t\right):={\mathcal{P}}_w\left(t\right)/{\rho }_w$ and $k_{w,cons}\left(t\right):={\mathcal{C}}_w\left(t\right)/{\rho }_w$ for every $t\ge 0$.

2.5.4. Time Evolution of Neutral Solutes

The time evolution of neutral solute molar density is described by the mass balance equation in integral form

$$\begin{array}{ccc}& {\int }_{{\Omega }_t}{\displaystyle \frac{\partial c_{\beta }}{\partial t}}d{\Omega }_t=-{\int }_{\partial {\Omega }_t}{j}_{\beta ,n}d(\partial {\Omega }_t)+{\int }_{{\Omega }_t}{R}_{\beta }d{\Omega }_t\beta \in S_{\beta },\hfill & \end{array}$$

(17)

where $j_{\beta ,n}$ is the normal molar flux density of $c_{\beta }$ over the cell surface ${\mathcal{S}}_t$, and

$$\begin{array}{ccc}& R_{\beta }(\mathbf{x},t):=P_{\beta }(\mathbf{x},t)-C_{\beta }(\mathbf{x},t)\mathbf{x}\in {\Omega }_t,t\ge 0\hfill & \end{array}$$

(18)

is the net production rate of intracellular solute $\beta$, $P_{\beta }$ and $C_{\beta }$ being production and consumption rates, respectively (units: $\mathrm{mM}{\mathrm{s}}^{-1}$).

Assumption 9.

We assume that:

$$\begin{array}{ccc}& c_{\beta }(\mathbf{x},t)=c^{\beta }\left(t\right)\Phi \left(\mathbf{x}\right)\mathbf{x}\in {\Omega }_t,t\ge 0,\hfill & \end{array}$$

(19a)

$$\begin{array}{ccc}& j_{\beta ,n}(\mathbf{x},t)=j_n^{\beta }\left(t\right){\eta }_{\beta }\left(\mathbf{x}\right)\mathbf{x}\in \partial {\Omega }_t,t\ge 0,\hfill & \end{array}$$

(19b)

where $c^{\beta }$ and $j_n^{\beta }$ are time-dependent intracellular neutral solute molar densities and normal molar flux densities over the cell surface, respectively, and the shape function ${\eta }_{\beta }$ is such that

$$\begin{array}{ccc}& {\int }_{\partial {\Omega }_t}{\eta }_{\beta }\left(\mathbf{x}\right)d(\partial {\Omega }_t)={\Phi }_{\beta }^{carr}{\mathcal{S}}_t,\hfill & \end{array}$$

(19c)

where ${\Phi }_{\beta }^{carr}$ is the surface fraction of the carrier protein that allows the transmembrane exchange of the neutral solute β. We also assume that equations similar to (15) apply to $P_{\beta }$ and $C_{\beta }$.

Using Assumptions A9 in (13c) we obtain the following reduced-order model for each neutral solute molar density:

$$\begin{array}{ccc}& {\displaystyle \frac{d{c}^{\beta }\left(t\right)}{dt}}=-\gamma j_{\beta ,n}^{tot}\left(t\right){\left(\mathcal{V}\left(t\right)\right)}^{-1/3}+{\mathcal{R}}_{\beta }\left(t\right)\beta \in S_{\beta },t>0,\hfill & \end{array}$$

(20a)

$$\begin{array}{ccc}& j_{\beta ,n}^{tot}\left(t\right):={\Phi }_{\beta }^{carr}{j}_n^{\beta }\left(t\right)\beta \in S_{\beta },t>0,\hfill & \end{array}$$

(20b)

$$\begin{array}{ccc}& c^{\beta }\left(0\right)=c_0^{\beta },\hfill & \end{array}$$

(20c)

where $c_0^{\beta }$ is the initial value of the intracellular neutral solute molar density, $\beta \in S_{\beta }$, and

$$\begin{array}{ccc}& {\mathcal{R}}_{\beta }\left(t\right):={\mathcal{P}}_{\beta }\left(t\right)-{\mathcal{C}}_{\beta }\left(t\right)\beta \in S_{\beta },t>0,\hfill & \end{array}$$

(20d)

is the net production rate of solute $\beta$ (units: $\mathrm{mM}{\mathrm{s}}^{-1}$).

2.5.5. Time evolution of charged solutes

The time evolution of charged solute molar density is described by the mass balance equation in integral form

$$\begin{array}{ccc}& {\int }_{{\Omega }_t}{\displaystyle \frac{\partial c_{\alpha }}{\partial t}}d{\Omega }_t=-{\int }_{\partial {\Omega }_t}{j}_{\alpha ,n}d(\partial {\Omega }_t)+{\int }_{{\Omega }_t}{R}_{\alpha }d{\Omega }_t\alpha \in S_{\alpha },\hfill & \end{array}$$

(21)

where $j_{\alpha ,n}$ is the normal molar flux density of $c_{\alpha }$ over the cell surface ${\mathcal{S}}_t$, and

$$\begin{array}{ccc}& R_{\alpha }(\mathbf{x},t):=P_{\alpha }(\mathbf{x},t)-C_{\alpha }(\mathbf{x},t)\mathbf{x}\in {\Omega }_t,t\ge 0\hfill & \end{array}$$

(22)

is the net production rate of intracellular ion $\alpha$, $P_{\alpha }$ and $C_{\alpha }$ being production and consumption rates, respectively (units: $\mathrm{mM}{\mathrm{s}}^{-1}$).

Assumption 10.

We assume that:

$$\begin{array}{ccc}& c_{\alpha }(\mathbf{x},t)=c^{\alpha }\left(t\right)\Phi \left(\mathbf{x}\right)\mathbf{x}\in {\Omega }_t,t\ge 0,\hfill & \end{array}$$

(23a)

$$\begin{array}{ccc}& j_{\alpha ,n}(\mathbf{x},t)=j_{\alpha ,n}^{ch}\left(t\right){\eta }_{\alpha }^{ch}\left(\mathbf{x}\right)+j_{\alpha ,n}^{exch}\left(t\right){\eta }_{\alpha }^{exch}\left(\mathbf{x}\right)+j_{\alpha ,n}^{pump}\left(t\right){\eta }_{\alpha }^{pump}\left(\mathbf{x}\right)\mathbf{x}\in \partial {\Omega }_t,t\ge 0,\hfill & \end{array}$$

(23b)

where $c^{\alpha }$ is the time-dependent intracellular molar density of ion $\alpha \in S_{\alpha }$, the functions $j_{\alpha ,n}^{ch}$ and $j_{\alpha ,n}^{exch}$ are time-dependent normal molar flux densities over the cell surface representing passive ion transport and exchange, respectively, the function $j_{\alpha ,n}^{pump}$ is a time-dependent normal molar flux density over the cell surface representing ion exchange through active pumps, and the shape functions ${\eta }_{\alpha }^{ch}$, ${\eta }_{\alpha }^{exch}$ and ${\eta }_{\alpha }^{pump}$ are such that:

$$\begin{array}{ccc}& {\int }_{\partial {\Omega }_t}{\eta }_{\alpha }^{ch}\left(\mathbf{x}\right)d(\partial {\Omega }_t)={\Phi }_{\alpha }^{ch}{\mathcal{S}}_t,\hfill & \end{array}$$

(23c)

$$\begin{array}{ccc}& {\int }_{\partial {\Omega }_t}{\eta }_{\alpha }^{exch}\left(\mathbf{x}\right)d(\partial {\Omega }_t)={\Phi }_{\alpha }^{exch}{\mathcal{S}}_t,\hfill & \end{array}$$

(23d)

$$\begin{array}{ccc}& {\int }_{\partial {\Omega }_t}{\eta }_{\alpha }^{pump}\left(\mathbf{x}\right)d(\partial {\Omega }_t)={\Phi }_{\alpha }^{pump}{\mathcal{S}}_t,\hfill & \end{array}$$

(23e)

where ${\Phi }_{\alpha }^{ch}$, ${\Phi }_{\alpha }^{exch}$ and ${\Phi }_{\alpha }^{pump}$ are the surface fractions of the membrane protein associated with passive transport, passive exchange and active pump-mediated exchange of ion α, respectively. We also assume that equations similar to (8) apply to $P_{\alpha }$ and $C_{\alpha }$.

Using Assumptions A10 in (21) we obtain the following reduced-order model for each charged solute molar density:

$$\begin{array}{ccc}& {\displaystyle \frac{d{c}^{\alpha }\left(t\right)}{dt}}=-\gamma j_{\alpha ,n}^{tot}\left(t\right){\left(\mathcal{V}\left(t\right)\right)}^{-1/3}+{\mathcal{R}}_{\alpha }\left(t\right)\alpha \in S_{\alpha },t>0,\hfill & \end{array}$$

(24a)

$$\begin{array}{ccc}& j_{\alpha ,n}^{tot}\left(t\right):=\left[{\Phi }_{\alpha }^{ch}{j}_{\alpha ,n}^{ch}\left(t\right)+{\Phi }_{\alpha }^{exch}{j}_{\alpha ,n}^{exch}\left(t\right)+{\Phi }_{\alpha }^{pump}{j}_{\alpha ,n}^{pump}\left(t\right)\right]\alpha \in S_{\alpha },t>0,\hfill & \end{array}$$

(24b)

$$\begin{array}{ccc}& c^{\alpha }\left(0\right)=c_0^{\alpha },\hfill & \end{array}$$

(24c)

where $c_0^{\alpha }$ is the initial value of the intracellular charged solute molar density, $\alpha \in S_{\alpha }$, and

$$\begin{array}{ccc}& {\mathcal{R}}_{\alpha }\left(t\right):={\mathcal{P}}_{\alpha }\left(t\right)-{\mathcal{C}}_{\alpha }\left(t\right)\alpha \in S_{\alpha },t>0,\hfill & \end{array}$$

(24d)

is the net production rate of ion $\alpha$ (units: $\mathrm{mM}{\mathrm{s}}^{-1}$).

2.5.6. Time evolution of membrane potential

For any $t\ge 0$, the membrane potential is defined as

$$\begin{array}{ccc}& {\psi }_m\left(t\right):={\psi }^{in}\left(t\right)-{\psi }^{out}\left(t\right).\hfill & \end{array}$$

(25)

The time evolution of the membrane potential is described by the charge balance equation in integral form

$$\begin{array}{ccc}& {\int }_{\partial {\Omega }_t}-{\epsilon }_m{\displaystyle \frac{\partial }{\partial t}}\left({\displaystyle \frac{\partial \psi }{\partial n}}\right)d(\partial {\Omega }_t)=-{\int }_{\partial {\Omega }_t}{J}_{cond,n}d(\partial {\Omega }_t)\hfill & \end{array}$$

(26)

where ${\displaystyle \frac{\partial \psi }{\partial n}}$ and $J_{cond,n}$ are the normal derivative of the electric potential and the normal total conduction current density over the cell surface ${\mathcal{S}}_t$, respectively.

Assumption 11.

We assume that the electric potential is a piecewise linear continuous function across the cell membrane thickness (see Figure 4).

Using Assumption 11 into the left-hand side of (26), we get

$$\begin{array}{ccc}& {\int }_{\partial {\Omega }_t}-{\epsilon }_m{\displaystyle \frac{\partial }{\partial t}}\left({\displaystyle \frac{\partial \psi }{\partial n}}\right)d(\partial {\Omega }_t)=c_m{\displaystyle \frac{d{\psi }_m\left(t\right)}{dt}}{\Phi }^{lip}\mathcal{S}\left(t\right)t\ge 0,\hfill & \end{array}$$

(27)

where

$$\begin{array}{ccc}& c_m:={\displaystyle \frac{{\epsilon }_m}{t_m}}\hfill & \end{array}$$

(28)

is the cell specific capacitance (units: $\mathrm{F}{\mathrm{m}}^{-2}$). The right-hand side of (26) is given by the following expression

$$\begin{array}{ccc}& -{\int }_{\partial {\Omega }_t}{J}_{cond,n}d(\partial {\Omega }_t)=-I_{cond}^{tot}\left(t\right)t\ge 0,\hfill & \end{array}$$

(29)

where $I_{cond}^{tot}$ is the total conduction current flowing across the cell membrane at time t (units: A), defined as

$$\begin{array}{ccc}& I_{cond}^{tot}\left(t\right)=\left[\sum _{\alpha \in S_{\alpha }}F{z}_{\alpha }{j}_{\alpha ,n}^{tot}\left(t\right)\right]\mathcal{S}\left(t\right)t\ge 0.\hfill & \end{array}$$

(30)

Replacing (27) and (29) into (26) we obtain the following reduced-order model for the membrane potential:

$$\begin{array}{ccc}& {\displaystyle \frac{d{\psi }_m\left(t\right)}{dt}}=-{\left({\Phi }^{lip}{c}_m\right)}^{-1}{j}_n^{\psi }\left(t\right)t>0,\hfill & \end{array}$$

(31a)

$$\begin{array}{ccc}& j_n^{\psi }\left(t\right):=j_{cond,n}^{tot}\left(t\right)t>0,\hfill & \end{array}$$

(31b)

$$\begin{array}{ccc}& {\psi }_m\left(0\right)={\psi }_{m,0},\hfill & \end{array}$$

(31c)

where

$$\begin{array}{ccc}& j_{cond,n}^{tot}\left(t\right):=\sum _{\alpha \in S_{\alpha }}F{z}_{\alpha }{j}_{\alpha ,n}^{tot}\left(t\right)t\ge 0\hfill & \end{array}$$

(31d)

is the normal total conduction current density over the cell surface defined in (24b), and ${\psi }_{m,0}$ is the initial value of the membrane potential.

2.6. Compact Form of the Cell Reduced-Order Model

Let us introduce the vector of dependent variables:

$$\begin{array}{ccc}& \mathbf{Y}\left(t\right):=\left[\begin{array}{c}\mathcal{V}\left(t\right)\\ {\mathbf{c}}_{\beta }\left(t\right)\\ {\mathbf{c}}_{\alpha }\left(t\right)\\ {\psi }_m\left(t\right)\end{array}\right]t\ge 0,\hfill & \end{array}$$

(32)

where:

$$\begin{array}{ccc}& {\mathbf{c}}_{\beta }\left(t\right):={\left[c_{\beta ,1}\left(t\right),\dots ,c_{\beta ,N_{\beta }}\left(t\right)\right]}^{T}t\ge 0,\hfill & \end{array}$$

(33)

$$\begin{array}{ccc}& {\mathbf{c}}_{\alpha }\left(t\right):={\left[c_{\alpha ,1}\left(t\right),\dots ,c_{\alpha ,N_{\alpha }}\left(t\right)\right]}^{T}t\ge 0,\hfill & \end{array}$$

(34)

are the column vectors containing the time values of neutral and charged solute molar densities.

Let us define the vector of source terms:

$$\begin{array}{ccc}& \mathbf{F}\left(t,\mathbf{Y}\left(t\right)\right):=\left[\begin{array}{c}{\displaystyle -\gamma j_{v,n}\left(t\right){\left(\mathcal{V}\left(t\right)\right)}^{2/3}+\left(k_{w,prod}\left(t\right)-k_{w,cons}\left(t\right)\right)\mathcal{V}\left(t\right)}\hfill \\ {\displaystyle -\gamma {\mathbf{j}}_{\beta ,n}^{tot}\left(t\right){\left(\mathcal{V}\left(t\right)\right)}^{-1/3}+\left({\mathbf{P}}_{\beta }\left(t\right)-{\mathbf{C}}_{\beta }\left(t\right)\right)}\hfill \\ {\displaystyle -\gamma {\mathbf{j}}_{\alpha ,n}^{tot}\left(t\right){\left(\mathcal{V}\left(t\right)\right)}^{-1/3}+\left({\mathbf{P}}_{\alpha }\left(t\right)-{\mathbf{C}}_{\alpha }\left(t\right)\right)}\hfill \\ {\displaystyle -{\left(c_m{\Phi }^{lip}\right)}^{-1}{j}_n^{\psi }\left(t\right)}\hfill \end{array}\right]t\ge 0,\hfill & \end{array}$$

(35)

where ${\mathbf{j}}_{\beta ,n}^{tot}$ and ${\mathbf{j}}_{\alpha ,n}^{tot}$ are column vectors of size $N_{\beta }$ and $N_{\alpha }$, respectively, containing the time values of neutral and charged solute normal molar flux densities, whereas ${\mathbf{P}}_{\beta }$ (and ${\mathbf{C}}_{\beta }$), ${\mathbf{P}}_{\alpha }$ (and ${\mathbf{C}}_{\alpha }$) are column vectors of size $N_{\beta }$ and $N_{\alpha }$, respectively, containing the time values of intracellular neutral and charged solute production (and consumption) terms.

Finally, let us define the column vector containing the initial condition of the model:

$$\begin{array}{ccc}& {\mathbf{Y}}_0:=\left[\begin{array}{c}{\mathcal{V}}_0\\ {\mathbf{c}}_0^{\beta }\\ {\mathbf{c}}_0^{\alpha }\\ {\psi }_{m,0}\end{array}\right],\hfill & \end{array}$$

(36)

where ${\mathbf{c}}_0^{\beta }$ and ${\mathbf{c}}_0^{\alpha }$ are column vectors of size $N_{\beta }$ and $N_{\alpha }$, respectively, which contain the values of the initial conditions for the intracellular neutral and charged solute molar densities.

The equation system constituting the model of cell volume motion can be written in compact form as:

$$\begin{array}{ccc}& {\displaystyle \frac{d\mathbf{Y}\left(t\right)}{dt}}=\mathbf{F}\left(t,\mathbf{Y}\left(t\right)\right)t>0,\hfill & \end{array}$$

(37a)

$$\begin{array}{ccc}& \mathbf{Y}\left(0\right)={\mathbf{Y}}_0.\hfill & \end{array}$$

(37b)

System (37) comprises differential equations and algebraic equations, represented by the constitutive laws for the normal molar flux densities of neutral and charged solutes and for the production and consumption rates of water and solutes. The expressions of the normal fluid velocity on the cell surface and normal molar flux densities are provided in Appendix A whereas the expressions of production and consumption rates are provided in Appendix B.

2.7. Numerical Approximation

The mathematical model proposed in this article is implemented within a Computational Virtual Laboratory (CVL) for the simulation of cell volume motion.

The $\theta$-method (see [18], Chapeter 3) and the Matlab tool ode suite are used for the numerical approximation of (37) (see [31]). In the case of $\theta$-method, the values $\theta =1$ and $\theta =0.5$ are utilized, $\theta =1$ corresponding to the Backward Euler (BE) method and $\theta =0.5$ corresponding to the Crank Nicolson (CN) method. In the case of the Matlab tool ode suite, the functions ode15s or ode23tb are utilized because they are specifically tailored to solve stiff problems like the one object of the present article (see [32], Chapeter 11). The function ode15s is endowed with a variable-order Backward Differentiation Formulae (BDF) method, whereas the function ode23tb uses the trapezoidal rule coupled with a BDF of order 3.

The CVL allows the user to consider physical situations of increasing complexity, starting from the solution of the sole Cauchy problem (16) in which the fluid velocity is a given function of time and the production and consumption terms are functions of time and cell volume. Simulation complexity can be increased by adding charged and neutral solutes, intracellular reactions and transmembrane ion exchange mechanisms, as well as the presence of impermeant protein charges in both intra and extracellular regions. This hierarchy of scenarios is investigated in Section 3 where the accuracy and reliability of model predictions is verified against analytical solutions and available data.

3. Results

In this section we validate the proposed CVL through the solution of case studies characterized by an increasing complexity.

3.1. The Basic Configuration

This case study is mathematically represented by the sole model for cell volume motion (16).

Assumption 12.

We make the following assumptions:

1. 

$v_{p,n}\left(t\right)=\overline{v}$ for all $t\ge 0$, $\overline{v}$ being a given constant;

2. 

$k_{w,prod}\left(t\right)={\kappa }_a$ for $t\ge 0$, ${\kappa }_a$ being a given positive constant (units: ${\mathrm{s}}^{-1}$);

3. 

$k_{w,cons}\left(t\right)={\kappa }_d{\displaystyle \frac{\mathcal{V}\left(t\right)}{{\mathcal{V}}_{ref}}}$ for $t\ge 0$, ${\kappa }_d$ (units: ${\mathrm{s}}^{-1}$) and ${\mathcal{V}}_{ref}$ (units: ${\mathrm{m}}^3$) being given positive constants.

Replacing the assumptions on the water volume production and consumption rates ino (16d) we get

$$\begin{array}{ccc}& {\mathcal{R}}_w\left(t\right)={\kappa }_a-{\kappa }_d{\displaystyle \frac{\mathcal{V}\left(t\right)}{{\mathcal{V}}_{ref}}}t\ge 0.\hfill & \end{array}$$

(38)

Replacing Assumptions 12 and (38) into (16a), we obtain the following Cauchy problem:

$$\begin{array}{ccc}& {\displaystyle \frac{d\mathcal{V}\left(t\right)}{dt}}=-\gamma {\Phi }^{AQP}\overline{v}{\left(\mathcal{V}\left(t\right)\right)}^{2/3}+\left({\kappa }_a-{\kappa }_d{\displaystyle \frac{\mathcal{V}\left(t\right)}{{\mathcal{V}}_{ref}}}\right)\mathcal{V}\left(t\right)t>0,\hfill & \end{array}$$

(39a)

$$\begin{array}{ccc}& \mathcal{V}\left(0\right)={\mathcal{V}}_{ref},\hfill & \end{array}$$

(39b)

where ${\mathcal{V}}_{ref}=(4\pi /3)R_{cell}^3$ is the cell volume in resting conditions.

The equilibrium points of system (39) are the solution of the following nonlinear algebraic equation

$$\begin{array}{ccc}& f\left(x\right)=-\gamma {\Phi }^{AQP}\overline{v}{x}^{2/3}+{\mathcal{V}}_{ref}^{1/3}x\left({\kappa }_a-{\kappa }_{d}x\right)=0,\hfill & \end{array}$$

(40)

where $x:={\mathcal{V}}_{\infty }/{\mathcal{V}}_{ref}$ is the dependent variable and ${\mathcal{V}}_{\infty }$ is the stationary value of the cell volume.

Figure 5 (left panel) shows a graph of f in correspondance of the following values of the input data: $R_{cell}=10·{10}^{-6}\mathrm{m}$, ${\Phi }^{AQP}=0.5$, $\overline{v}=-30·{10}^{-6}\mathrm{m}{\mathrm{s}}^{-1}$, ${\kappa }_a=1{\mathrm{s}}^{-1}$ and $k_d=3{\mathrm{s}}^{-1}$. System (39) admits two equilibrium points: ${\mathcal{V}}_{\infty ,1}=0{\mathrm{m}}^3$ and ${\mathcal{V}}_{\infty ,2}=1.6124{\mathcal{V}}_{ref}=6.754·{10}^{-15}{\mathrm{m}}^3$. We have

$$\begin{array}{ccc}& {\displaystyle \frac{df}{dx}}{|}_{x=0}=+\infty ,{\displaystyle \frac{df}{dx}}{|}_{x=1.6124}=-9.86·{10}^{-5}\mathrm{m}{\mathrm{s}}^{-1},\hfill & \end{array}$$

from which we can conclude that ${\mathcal{V}}_{\infty ,2}$ is the only stable equilibrium point of system (39). The theoretical conclusions of the stability analysis are confirmed by Figure 5 (left panel) which shows the normalized cell volume computed by solving system (39) in the time interval $[0,10]\mathrm{s}$ with the following values of model parameters: $R_{cell}=10·{10}^{-6}\mathrm{m}$, ${\Phi }^{AQP}=0.5$, $\overline{v}=-30·{10}^{-6}\mathrm{m}{\mathrm{s}}^{-1}$, ${\kappa }_a=1{\mathrm{s}}^{-1}$ and $k_d=3{\mathrm{s}}^{-1}$. Cell volume dynamics was studied using the BE method on a uniform time partition made of ${10}^5$ elements. Consistent with the physical configuration illustrated in Figure 1, the cell tends to increase its volume by more than 60% because of the inflow of water from the extracellular region. Cell volume reaches a stationary limit because of dominating intracellular water consumption over intracellular water production.

Figure 6 (left panel) illustrates the temporal evolution of cell volume in the case where the values of $\overline{v}$ are $[-30-18,-6,+6,+18,+30]·{10}^{-6}\mathrm{m}{\mathrm{s}}^{-1}$. The black arrow in each figure indicates the increase of $\overline{v}$, from negative to positive values. Results suggest that the cell switches between swelling and shrinking as fluid velocity changes its sign from negative to positive. Interestingly, the BE method gives rise to positive cell volumes for each $t\ge 0$ and for each considered value of $\overline{v}$. This outcome is the consequence of the positivity principle enjoyed by the BE method when applied to the linear equation $y^{\prime }\left(t\right)=-\lambda y\left(t\right)$ for $t\ge 0$, $\lambda$ being a given positive constant. The difference between using the BE method ($\theta =1$) and another $\theta$-method with $\theta \in [0,1)$ is illustrated in Figure 6 (right panel) which shows a comparison between the BE method (red line) and the Crank Nicolson method (CN, blue line) in the solution of (39) when the time interval is $[0,1]\mathrm{s}$, the fluid velocity is $\overline{v}=100·{10}^{-6}\mathrm{m}{\mathrm{s}}^{-1}$ and the number of time elements is 10, corresponding to a time step $\Delta t=0.1\mathrm{s}$. The value of the remaining model parameters is the same as in the previous example. Results indicate that the solution computed by the CN method is affected by spurious unphysical oscillations unlike that computed by the BE method. Such oscillations can be removed by increasing the number of discretization time elements, at the price of increasing the computational effort.

Figure 7 (left panel) illustrates the temporal evolution of the water volume net production rate ${\mathcal{R}}_w$ in the time interval $[0,10]\mathrm{s}$ with the fluid velocity varying in the range $[-30+30]·{10}^{-6}\mathrm{m}{\mathrm{s}}^{-1}$. We see that the stationary value of ${\mathcal{R}}_w$, for each considered normal fluid velocity in the range $[-30+30]·{10}^{-6}\mathrm{m}{\mathrm{s}}^{-1}$, changes from negative to positive, the magnitude for $v_n=-30·{10}^{-6}\mathrm{m}{\mathrm{s}}^{-1}$ being three times larger than for $v_n=+30·{10}^{-6}\mathrm{m}{\mathrm{s}}^{-1}$. Figure 7 (right panel) illustrates the temporal evolution of the cell normalized volume in the time interval $[0,1]\mathrm{s}$ (ten times smaller than before) in the case where $\mathcal{R}\left(t\right)=0$ for $t\ge 0$ and the fluid velocity varies in the range $[-30+30]·{10}^{-6}\mathrm{m}{\mathrm{s}}^{-1}$. We see that in the absence of intracellular regulatory mechanisms, model system (39) predicts an abnormal increase of cell volume for a highly negative value of normal fluid velocity.

3.2. Cell Homeostasis in the Ciliary Epithelium of the Eye

This case study is mathematically represented by the Cauchy problem (37) in which the cell represents one of the pigmented/nonpigmented couplet in the ciliary epithelium (CE) of the ciliary body of the eye (see [7,8,33]). The aim of the study is to apply the model proposed in this article to characterize the homeostatic configuration of the cell under physiological conditions of the system. According to the data reported in [7], such conditions correspond to:

(i) 

volume of the CE equal to $8\mu L$;

(ii) 

number of the cell couplets $N_{cells,CE}$ constituting the CE equal to 4 millions;

(iii) 

intraocular pressure equal to $15mmHg$;

(iv) 

AH volumetric flow rate equal to $2.75\mu \mathrm{L}\mathrm{m}\mathrm{i}{\mathrm{n}}^{-1}$.

Assumption 13.

We make the following assumptions:

1. 

the considered sets of neutral and charged solutes are:

$$\begin{array}{ccc}& S_{\beta }=\left\{{\mathrm{CO}}_2,{\mathrm{H}}_2{\mathrm{CO}}_3\right\},\hfill & \end{array}$$

(41)

$$\begin{array}{ccc}& S_{\alpha }=\left\{{\mathrm{Na}}^{+},{\mathrm{K}}^{+},{\mathrm{H}}^{+},{\mathrm{Cl}}^{-},{\mathrm{HCO}}_3^{-}\right\};\hfill & \end{array}$$

(42)

2. 

the molar densities of neutral and charged solutes are given constants denoted by ${\overline{c}}_{\beta ,ex}$, $\beta \in S_{\beta }$, and ${\overline{c}}_{\alpha ,ex}$, $\alpha \in S_{\alpha }$;

3. 

the hydraulic pressure difference is a given constant denoted by $\Delta \overline{p}$;

4. 

no transmembrane ion exchangers are considered, so that ${\Phi }^{exch}=0$;

5. 

the model of carrier membrane proteins is described in Appendix A.2;

6. 

the model of ion channels is described in Appendix A.3;

7. 

the model of net production rates for neutral and charged solutes is described in Appendix B.1;

8. 

the model of net production rate in cell volume regulation is described in Appendix B.2.

The values of all model parameters that are used in the computer simulations illustrated in this section are listed in Table 1.

The vector containing the initial condition data is

$$\begin{array}{ccc}& {\mathbf{Y}}_0={\left[2·{10}^{-15},0,0,10,140,5.0119·{10}^{-5},7,{10}^{-4},-8.0097·{10}^{-2}\right]}^T\hfill & \\ & \mathrm{units}:{\left[{\mathrm{m}}^3,\mathrm{mM},\mathrm{mM},\mathrm{mM},\mathrm{mM},\mathrm{mM},\mathrm{mM},\mathrm{mM},\mathrm{V}\right]}^T.\hfill & \end{array}$$

The vectors containing the extracellular values of the solute molar density are:

$$\begin{array}{cccc}& {\overline{\mathbf{c}}}_{\beta }={\left[0.5,1.3·{10}^{-3}\right]}^T\hfill & \hfill \mathrm{units}:{\left[\mathrm{mM},\mathrm{mM}\right]}^T,& \\ & {\overline{\mathbf{c}}}_{\alpha }={\left[130,5,3.1623·{10}^{-5},150,9.9\right]}^T\hfill & \hfill \mathrm{units}:{\left[\mathrm{mM},\mathrm{mM},\mathrm{mM},\mathrm{mM},\mathrm{mM}\right]}^T.& \end{array}$$

The values of the remaining model parameters can be found in Appendix B.1, Appendix B.2 and Appendix B.3. Computations are done using the Matlab solver ode15s setting both relative and absolute tolerances equal to ${10}^{-12}$.

3.2.1. Electroneutrality and Impermeant Charged Proteins

Electroneutrality is essential in the maintainance of cell volume as a function of time and of phenomena occurring in the intracellular and extracellular regions (see [3] and references cited therein). Computation of the total electric charge densities ${\rho }_{mob,in}$ and ${\rho }_{mob,ex}$ due to mobile ions inside and outside the cell at resting conditions yields:

$$\begin{array}{ccc}& {\rho }_{mob,in}=F\sum _{\alpha \in S_{\alpha }}{z}_{\alpha }{c}_{\alpha ,in}\left(0\right)=1.3797·{10}^7\mathrm{C}{\mathrm{m}}^{-3},\hfill & \\ & {\rho }_{mob,ex}=F\sum _{\alpha \in S_{\alpha }}{z}_{\alpha }{\overline{c}}_{\alpha }=-9.5519·{10}^5\mathrm{C}{\mathrm{m}}^{-3}.\hfill & \end{array}$$

These results indicate that the intracellular region (at $t=0$) is characterized by a high excess of positive charge whereas the extracellular region (at $t=0$) is characterized by a high excess of negative charge. This charge difference gives rise to a large osmotic pressure difference across the cell membrane which may eventually lead to disruption of cell integrity. To neutralize the excess of positive and negative charge, we need the presence of internally sequestered impermeant charges inside and outside the cell. Denoting by $c_X$ and $z_X$ the molar density and charge number of the impermeant charge, we have:

$$\begin{array}{ccc}& c_{X,in}=143\mathrm{mM}{z}_{X,in}=-1,\hfill & \end{array}$$

(43a)

$$\begin{array}{ccc}& c_{X,ex}=9.9\mathrm{mM}{z}_{X,ex}=+1.\hfill & \end{array}$$

(43b)

These results indicate that a high molar density of fixed anions is sequestered inside the cell cytoplasm whereas a much smaller molar density of cations is required to make the extracellular solution electroneutral.

Assumption 14.

We assume that the number of moles of intracellular impermeant charge $n_X$ is constant during the time evolution of the cell.

Let $c_X\left(0\right)$ denote the intracellular molar density of the impermeant charge at time $t=0$ (units: $mM$). By definition of molar density, we have

$$\begin{array}{ccc}& c_X\left(0\right)={\displaystyle \frac{n_X}{{\mathcal{V}}_0}}.\hfill & \end{array}$$

(44)

Using Assumption 14 and (44), we can express the intracellular impermeant charge molar density for any time $t\ge 0$ as

$$\begin{array}{ccc}& c_{X,in}\left(t\right)={\displaystyle \frac{n_X}{\mathcal{V}\left(t\right)}}=c_X\left(0\right){\displaystyle \frac{{\mathcal{V}}_0}{\mathcal{V}\left(t\right)}}t\ge 0.\hfill & \end{array}$$

(45)

Simulations illustrated in the next sections are conducted using the values of $z_X$ and $c_X$ in (43) and the constitutive equation (45).

3.2.2. Fast Time Scale Cell Evolution

We investigate the evolution of the cell in the time interval $[t_0,t_{end}]$, where $t_0=0\mathrm{s}$ and $t_{end}=50·{10}^{-12}\mathrm{s}$.

Figure 8 (left panel) illustrates the time evolution of the intracellular protonated hydrogen (blue curve) and the corresponding intracellular pH (red curve). Hydrogen fast diffusion from the intracellular side into the extracellular side determines a sharp decrease of $c_{{\mathrm{H}}^{+},in}$ so that the cytoplasm solution turns into a very basic condition (the maximum value of intracellular pH is 12.5). Figure 8 (right panel) illustrates the time evolution of the average cell normal surface velocity (blue curve) and the corresponding percentage variation $\Delta {\mathcal{V}}_{\%}$ of the cell volume with respect to the initial condition (red curve). Cell surface velocity is positive and very small in magnitude (less than $0.02\mu \mathrm{m}{\mathrm{s}}^{-1}$) so that the volume of the cell experiences a very little decrease (maximum magnitude equal to $2·{10}^{-11}\%$) with respect to resting conditions. Figure 8 (middle panel, bottom) illustrates the time evolution of the total AH volumetric flow rate $Q_{\mathrm{AH}}$ throughout the CE (units: $\mu \mathrm{L}{\min }^{-1}$). The quantity $Q_{\mathrm{AH}}\left(t\right)$ is computed for every $t\in [t_0,t_{end}]$ using the following relation

$$\begin{array}{ccc}& Q_{\mathrm{AH}}\left(t\right)=\left(N_{cells,CE}{v}_{cell,n}\left(t\right)\mathcal{S}\left(t\right)\right)60·{10}^9\mathrm{units}:\mu \mathrm{L}{\min }^{-1}t\in [t_0,t_{end}].\hfill & \end{array}$$

(46)

In less than 50ps, $Q_{\mathrm{AH}}$ reaches almost 72% of the volumetric flow rate that is physiologically expected in a normal tension individual (see condition [(iv]).

3.2.3. Medium Time Scale Cell Evolution

We investigate the evolution of the cell in the time interval $[t_0,t_{end}]$, where $t_0=0s$ and $t_{end}=5\mathrm{s}$.

Figure 9 (top left panel) illustrates the time evolution of the intracellular CO₂ (blue curve) and H₂CO₃ (red curve) molar densities. In less than 0.5s the hydration process gives rise to a significant production of carbon dioxide and carbonic acid, which is eventually followed up by a stationary condition corresponding to the dynamic equilibrium of reaction (A6a). The consequence of the CO₂ hydration can also be seen from Figure 9 (top right panel) that illustrates the time evolution of the intracellular H⁺ molar density and the corresponding pH (blue and red curves, respectively). Protonated hydrogen concentration increases significantly until almost $t=0.125\mathrm{s}$ because of carbonic acid dissociation (forward reaction in (A6b)). Then, the association reaction (backward reaction in (A6b)) with bicarbonate gives rise to a decrease of $c_{{\mathrm{H}}^{+},in}$ until a stationary condition is reached. In such a condition, the value of the intracellular pH (almost 6) indicates that the cytoplasm solution is acid. Figure 9 (bottom left panel) illustrates the time evolution of the average cell normal surface velocity (blue curve) and the corresponding percentage variation $\Delta {\mathcal{V}}_{\%}$ of the cell volume with respect to the initial condition (red curve). As in the fast scale cell evolution, the cell surface velocity is positive and very small in magnitude (less than $0.02\mu \mathrm{m}{\mathrm{s}}^{-1}$). In this case, however, the much larger time duration of the analysis (5 seconds instead of 50ps) allows the cell to decrease by a larger percentage amount with respect to resting conditions (less than $0.5\%$ instead of $2·{10}^{-11}\%$). Figure 9 (bottom right panel) illustrates the time evolution of the total AH volumetric flow rate $Q_{\mathrm{AH}}$ throughout the CE (units: $\mu \mathrm{L}{\min }^{-1}$). In 5s, $Q_{\mathrm{AH}}$ reaches more than 78% of the volumetric flow rate that is physiologically expected in a normal tension individual.

3.2.4. Long Time Scale Cell Evolution

We investigate the evolution of the cell in the time interval $[t_0,t_{end}]$, where $t_0=0\mathrm{s}$ and $t_{end}=5400\mathrm{s}$. The value of $t_{end}$ corresponds to the time that is needed by the eye to completely replace the AH content of the anterior chamber (see [7] and [9]).

Figure 10 (top left panel) illustrates a zoomed view of the time evolution of the membrane potential in the time interval $t\in [0,10]\mathrm{s}$. After an initial ultra-fast transient due to the mismatch between the intracellularly applied initial conditions and the conditions in the extracellular bath, the cell reaches a stationary state of $-85.9\mathrm{mV}$, corresponding to an after-hyperpolarization of $-5.9\mathrm{mV}$. Figure 10 (top right panel) illustrates a zoomed view of the time evolution of the intracellular molar densities of sodium (blue curve), potassium (red curve) and chlorine (green curve) for $t\in [0,120]\mathrm{s}$. The concentration of sodium experiences a significant depletion (from 10 to 5.4 mM) because of the NaK ATPase pump activity. Similarly, the concentration of potassium increases from 140 mM up to a stationary value of almost $144\mathrm{mM}$. The green curve in Figure 10 (top right panel) indicates that the concentration of intracellular chlorine experiences a depletion from 7mM to 5.3mM. This can be explained by Figure 10 (middle left panel) which illustrates a zoom of the chlorine molar flux density $j_{{\mathrm{Cl}}^{-},n}\left(t\right)$ (units: $\mathrm{mM}\mathrm{m}{\mathrm{s}}^{-1}$) for $t\in [0,10]\mathrm{s}$. The positive value of $j_{{\mathrm{Cl}}^{-},n}$ indicates that chlorine is swept out of the cell cytoplasm with a progressively reducing magnitude along time. Figure 10 (middle right panel illustrates a zoom of the time evolution of the average cell normal surface velocity (blue curve) and the corresponding percentage variation $\Delta {\mathcal{V}}_{\%}$ of the cell volume with respect to the initial condition (red curve) for $t\in [0,120]\mathrm{s}$. As in the previous conditions the cell surface velocity is positive and very small in magnitude (less than $0.015\mu \mathrm{m}{\mathrm{s}}^{-1}$). The stationary cell volume percentage decrease with respect to resting conditions is less than $0.57\%$. Figure 10 (bottom center panel) illustrates a zoom of the time evolution of the total AH volumetric flow rate $Q_{\mathrm{AH}}$ throughout the CE (units: $\mu \mathrm{L}{\min }^{-1}$) for $t\in [0,300]\mathrm{s}$. The stationary value of $Q_{\mathrm{AH}}$ is 2.75017 $\mu \mathrm{L}{\min }^{-1}$, with a percentage error of -0.0062 % with respect to the physiological value of 2.75 $\mu \mathrm{L}{\min }^{-1}$.

4. Discussion

Cell volume stabilization is a fundamental function to preserve the mechanical integrity of a cell [3]. Keeping the cell volume constant along time is the result of the balance among forces acting on the cell membrane from both intra and extracellular sides, which emanate from electrochemical and fluid-mechanical mechanisms [4]. The experimental characterization of force equilibrium at the scale of the cell membrane is a nontrivial task, on the one hand, because of the difficulty of performing measurements in vivo (see [5]), and, on the other hand, the presence of two distinct regimes of mechanical behavior in living cells. In fact, a combination of experimentation and theoretical analysis shows that living mammalian cytoplasm behaves as an equilibrium material at short time scales whereas it behaves as an out-of-equilibrium material at long time scales [6].

Based on these considerations, we address in this article the theoretical study of cell volume dynamics with the aim of applying the formulation to the process of production of aqueous humor (AH) in the eye. AH is a slowly moving fluid that is (1) continuously produced by the ciliary body of the eye; (2) flows from the posterior into the anterior chamber of the eye; (3) is drained out of the eye throughout two main outflow pathways (see [7,8,9,10]). The sequence: (1) AH production; (2) AH flow; and (3) AH outflow is schematically represented in Figure 11.

The main function of AH is to keep the anterior segment of the eye clean from the metabolic waste of the surrounding tissues and to preserve the spherical shape of the eyeball by establishing the intraocular pressure (IOP) as the value of AH fluid pressure in correspondance of which the volumetric flow rate of AH that is produced by the ciliary body equals the volumetric flow rate of AH that is drained out by the outflow pathways in the anterior segment of the eye (see [35]). Successful control of the IOP within normal levels (10 to 21 mmHg) is central to prevent the occurrence and progression of neurodegenerative diseases of the retina such as glaucoma [11,12,13]. Therefore, unraveling the biophysical mechanisms that concur to determine the value of IOP may be of significant help to understand the causes of glaucoma and to develop effective therapies for its cure.

The development of a mathematical model and of a CVL for the simulation of the process of production, flow and outflow of AH has been subject of investigation in recent years (see [14,15,16,17]). Our proposed formulation is characterized by the following features: (F.1) it is based on the use of homogeneous mixtures including neutral and charged solutes (see [18], Chapter 13); (F.2) it is defined at the level of the single cell; and (F.3) it utilizes a model reduction procedure from three spatial dimensions to zero spatial dimensions. Feature (F.1) confers a solid theoretical foundation to the proposed model. Features (F.2) and (F.3) make the model structure simple and the computational schemes fast and suitable for adoption in a clinical environment. The model is a consistent generalization of previous approaches [3,4] as it shares with them the same conceptual simplifying assumption of working at the level of the single cell. This assumption is applied here to evaluate the collective behaviour of the cells in the CE of the eye in the process of production of AH.

The model that we propose in these pages has also limitations: (L.1) the dependence of all the variables from the spatial coordinate is neglected; (L.2) several important transmembrane mechanisms regulating solute exchange are not considered in the simulations; (L.3) statistical variability of parameters is not accounted for. Limitation (L.1) is a consequence of the use of the 3D-to-0D reduction procedure. Limitation (L.2) is a choice to prevent proliferation of model parameters thereby rendering easier the analysis of simulation predictions. Limitation (L.3) is a consequence of the choice of using a mechanistic (continuum-based) approach. We intend to remove all of these limitations in future extensions of the formulation considered in the present article.

4.1. The Impact of the Oncotic Pressure Due to the Impermeant Charge

In this simulation we solve system (37) in the time interval $[0,600]s$ with the same set of parameters as in Section 3.2, except the reflection coefficient ${\sigma }_X$ which is set equal to $[0:0.2:1]$ (Matlab vector notation). By doing so, the weight of the contribution of the oncotic pressure difference (A3h) to the total osmo-oncotic pressure difference (A3k) increases progressively from 0% (${\sigma }_X=0$) to 100 % (${\sigma }_X=1$), thereby allowing us to investigate the impact of the oncotic pressure difference on AH simulation.

Figure 12 (left panel) illustrates the time evolution of the total volumetric flow rate $Q_{\mathrm{AH}}\left(t\right)$ as a function of time t in the interval $t\in [0,600]\mathrm{s}$. Results indicate that the smaller ${\sigma }_X$, the larger the predicted total volumetric flow rate of AH. In particular, the value of $Q_{\mathrm{AH}}$ predicted by the model that is compatible with the given intracellular and extracellular fluid pressures (20mmHg and 15 mmHg, respectively) is obtained for ${\sigma }_X=1$. To better understand the effect of including properly the impermeant charge in AH modeling, we illustrate in Figure 12 (right panel) the time evolution of the total osmo-oncotic pressure difference $\Delta \Pi \left(t\right)$ in the interval $t\in [0,600]\mathrm{s}$. We see that the magnitude of $\Delta \Pi$ largely exceeds the contribution from hydrostatic pressure difference (equal to $5\mathrm{mmHg}$) for every value of ${\sigma }_X\in [0:0.2:1]$. Moreover, for every ${\sigma }_X\in [0:0.2:1]$ the oncotic pressure difference $\Delta {\Pi }_X$ is always positive whereas the total osmotic pressure difference $\Delta {\Pi }_{osm}$ is always negative, so that, as ${\sigma }_X$ increases, the magnitude of the total osmo-oncotic pressure difference decreases, to reach the value of $-1500\mathrm{mmHg}$.

4.2. The impact of the Na⁺/K⁺ ATPase

In this simulation we solve system (37) in the time interval $[0,5400]s$ with the same set of parameters as in Section 3.2, except the amplification coefficient $M_{ATP}$ in (A10e) which is set equal to $[0:0.5:2]$ (Matlab vector notation). By doing so, the weight of the contribution of the sodium/potassium pump to the electrochemical balance of the cell increases progressively from 0% ($M_{ATP}=0$) to 100 % ($M_{ATP}=2$), with respect to the working conditions of Section 3.2, thereby allowing us to investigate the impact of the Na⁺/K⁺ ATPase on AH simulation.

Figure 13 illustrates the time evolution of the total volumetric flow rate $Q_{\mathrm{AH}}\left(t\right)$ as a function of time t in the interval $t\in [0,600]\mathrm{s}$. Results indicate that for $M_{ATP}<1$ (corresponding to values of $c_{ATP}<c_{ATP,ref}$, the predicted total volumetric flow rate of AH is negative. This means that the pump does not have enough energy to move sodium out from the cell and potassium into the cell, so that the accumulation of sodium in the cytoplasm attracts chlorine from the extracellular region, eventually leading to the inversion of fluid flow. The predicted total volumetric flow rate of AH becomes positive for $M_{ATP}=1$. This means that the pump does have enough energy to move sodium out of the cell and potassium into the cell, preventing chlorine accumulation in the cell cytoplasm. For increasing values of the ATP concentration ($M_{ATP}>1$) water flux is favored, reaching a physiological level when $M_{ATP}=2$.

4.3. The impact of carbonic anhydrase

In this simulation we solve system (37) in the time interval $[0,5400]s$ with the same set of parameters as in Section 3.2, except the amplification coefficient $A_{\mathrm{CA}}$ in Equations (A7c)- (A7d) which is set equal to $[0:1:5]$ (Matlab vector notation). By doing so, the weight of the contribution of the CA enzyme to improve the reaction rate of CO₂ hydration increases progressively from 0% ($A_{\mathrm{CA}}=0$) to 100 % ($A_{\mathrm{CA}}=5$), with respect to the working conditions of Section 3.2, thereby allowing us to investigate the impact of the CA enzyme on AH simulation.

Figure 14 illustrates a zoomed detail of the time evolution of the total volumetric flow rate $Q_{\mathrm{AH}}\left(t\right)$ as a function of time t in the interval $t\in [0,600]\mathrm{s}$. Results indicate that $Q_{\mathrm{AH}}\left(t\right)$ is practically independent of the concentration of the CA enzyme. Further tests with higher values of the amplification parameter $A_{\mathrm{CA}}$ do not show significant variation of the predicted value of $Q_{\mathrm{AH}}\left(t\right)$, this probably to be ascribed to the low value of the intracellular carbon dioxide molar density (cf. Figure 9).

4.4. The impact of IOP

In this simulation we solve system (37) in the time interval $[0,5400]\mathrm{s}$ with the same set of parameters as in Section 3.2, except the value of IOP which is taken in the range [15,   150]mmHg. By doing so, we want to investigate the response of the model in the case of highly hypertensive patients.

Figure 15 illustrates the % difference between $Q_{\mathrm{AH}}$ at 15mmHg and $Q_{\mathrm{AH}}$ evaluated in correspondence of IOP in the interval [15, 150]mmHg. The trend is linear with respect to IOP, which agrees with the following expression for fluid velocity

$$\begin{array}{ccc}& v_f=L_p\left(p_{\mathrm{in}}\left(t\right)-p_{\mathrm{ex}}\left(t\right)-\Delta \Pi \left(t\right)\right)t\ge 0,\hfill & \end{array}$$

(47)

where $p_{\mathrm{in}}$=20 mmHg is the intracellular fluid pressure (corresponding to the estimated fluid pressure in the CE upon assuming 25mmHg in the ciliary capillaries), and $\Delta \Pi$ is the total osmo-oncotic difference. Since $\Delta \Pi \left(t\right)$ turns out to be negative with respect to t, as IOP increases, we see from (47) that $v_f$ diminishes, which explains the fact that also $Q_{\mathrm{AH}}$ diminishes. We notice that for IOP=30 mmHg, which is a typical value of the intraocular pressure in hypertensive individuals, the % difference between $Q_{\mathrm{AH}}$ at 15mmHg and $Q_{\mathrm{AH}}$ is 1.09% which increases to 2.54% for IOP=50 mmHg.

5. Conclusions

In this article we have proposed, analyzed and numerically investigated a reduced-order mathematical description of cellular volume homeostasis. The model accounts for intracellular reactions and transmembrane mechanisms for neutral and charged solute exchange. Hydrostatic and osmo-oncotic pressure differences are used in the context of Starling’s Law to compute the velocity of the fluid which drives the motion and radial deformation of the cell volume.

The model is implemented within the context of a CVL that is applied to the study of the process of AH production in the human eye. The scope of the simulations is to test the potential of the CVL to assess the relative quantitative importance of the biophysical mechanisms that underlie AH production, in the perspective of its use as a supporting tool to integrate and complement in vivo experiments and artificial intelligence-based methodologies for the analysis of data on a statistically significant number of patients.

The main conclusions of the model simulations suggest the importance of impermeant charged proteins and Na⁺/K⁺ ATPase on the level and direction of AH production, whereas AH production seems to be independent on CA concentration, at least for low values of CO₂ concentration.

Appendix A.1. Normal fluid velocity

To determine the normal fluid velocity on the cell surface we consider the following two cases:

(B.C.1) 

The hydraulic pressure difference $\Delta \overline{p}\left(t\right)$ is given (units: $\mathrm{Pa}$).

(B.C.2) 

The total water volumetric flow rate across the cell surface ${\overline{Q}}_{w,cell}\left(t\right)$ is given (units: ${\mathrm{m}}^3{\mathrm{s}}^{-1}$).

In case (B.C.2) we use (11b) to obtain

$$\begin{array}{ccc}& v_{cell,n}\left(t\right)={\displaystyle \frac{{\overline{Q}}_{w,cell}\left(t\right)}{\mathcal{S}\left(t\right)}}={\displaystyle \frac{{\overline{Q}}_{w,cell}\left(t\right)}{\gamma {\left(\mathcal{V}\left(t\right)\right)}^{2/3}}}t\ge 0.\hfill & \end{array}$$

(A1)

In case (B.C.1) we use (10d) to obtain

$$\begin{array}{ccc}& v_{cell,n}\left(t\right)={\Phi }^{AQP}{v}_{p,n}\left(t\right)={\Phi }^{AQP}{\mathcal{L}}_p\left[\Delta \overline{p}\left(t\right)-\Delta \Pi \left(t\right)\right]t\ge 0.\hfill & \end{array}$$

(A2)

To determine $v_{cell,n}$ we need to provide a mathematical model for the total osmo-oncotic pressure difference $\Delta \Pi$. Let us define the electrochemical potential of ion $\alpha$

$$\begin{array}{ccc}& {\phi }_{\alpha }^{ec}=\psi +{\displaystyle \frac{V_{th}}{z_{\alpha }}}ln\left({\displaystyle \frac{c_{\alpha }}{c_{ref}}}\right)\alpha \in S_{\alpha },\hfill & \end{array}$$

(A3a)

where $c_{ref}$ is a reference molar density. Let us introduce the pressure differences:

$$\begin{array}{ccc}& \Delta {\Pi }_{\alpha }\left(t\right):={\int }_0^{t_m}-{\sigma }_{\alpha }F{z}_{\alpha }{c}_{\alpha }(s,t){\displaystyle \frac{\partial {\phi }_{\alpha }^{ec}(s,t)}{\partial s}}ds\alpha \in S_{\alpha },t\ge 0,\hfill & \end{array}$$

(A3b)

$$\begin{array}{ccc}& \Delta {\Pi }_{\beta }\left(t\right):={\int }_0^{t_m}-{\sigma }_{\beta }RT{\displaystyle \frac{\partial c_{\beta }(s,t)}{\partial s}}ds\beta \in S_{\beta }t\ge 0,\hfill & \end{array}$$

(A3c)

$$\begin{array}{ccc}& \Delta {\Pi }_X\left(t\right):={\int }_0^{t_m}-{\sigma }_{X}RT{\displaystyle \frac{\partial c_X(s,t)}{\partial s}}dst\ge 0.\hfill & \end{array}$$

(A3d)

The quantities $\Delta {\Pi }_{\alpha }$ and $\Delta {\Pi }_{\beta }$ are the osmotic pressure differences associated with the charged solute $\alpha \in S_{\alpha }$ and neutral solute $\beta \in S_{\beta }$, respectively, whereas $\Delta {\Pi }_X$ is the oncotic pressure difference associated with the impermeant charge $c_X$. The quantity R is the gas constant (units: $\mathrm{J}{\mathrm{mol}}^{-1}{\mathrm{K}}^{-1}$), and ${\sigma }_{\alpha }$, ${\sigma }_{\beta }$ and ${\sigma }_X$ are the membrane reflection coefficients associated with ion $\alpha$, solute $\beta$ and impermeant fixed charge $c_X$, respectively. We notice that each coefficient is in the range $[0,1]$, the value 0 corresponding to a completely permeable membrane and the value 1 corresponding to a completely impermeable membrane. To evaluate the integral on the right-hand side of (A3b) we replace $c_{\alpha }(s,t)$ with its spatial average

$$\begin{array}{ccc}& \langle c_{\alpha }\rangle \left(t\right):={\displaystyle \frac{{\int }_0^{t_m}{c}_{\alpha }(s,t)ds}{t_m}}t\ge 0.\hfill & \end{array}$$

(A3e)

Postponing to Appendix A.3 the computation of (A3e) and replacing it into (A3b), we get:

$$\begin{array}{ccc}& \Delta {\Pi }_{\alpha }\left(t\right)={\sigma }_{\alpha }F{z}_{\alpha }\langle c_{\alpha }\rangle \left(t\right)\Delta {\phi }_{\alpha }^{ec}\left(t\right)\alpha \in S_{\alpha }t\ge 0,\hfill & \end{array}$$

(A3f)

$$\begin{array}{ccc}& \Delta {\Pi }_{\beta }\left(t\right)={\sigma }_{\beta }RT\Delta c_{\beta }\left(t\right)\beta \in S_{\beta }t\ge 0,\hfill & \end{array}$$

(A3g)

$$\begin{array}{ccc}& \Delta {\Pi }_X\left(t\right)={\sigma }_{X}RT\Delta c_X\left(t\right)t\ge 0.\hfill & \end{array}$$

(A3h)

We define the total osmotic pressure difference as

$$\begin{array}{ccc}& \Delta {\Pi }_{osm}\left(t\right)=\sum _{\alpha \in S_{\alpha }}\Delta {\Pi }_{\alpha }\left(t\right)+\sum _{\beta \in S_{\beta }}\Delta {\Pi }_{\beta }\left(t\right)t\ge 0,\hfill & \end{array}$$

(A3i)

and the oncotic pressure difference as

$$\begin{array}{ccc}& \Delta {\Pi }_{onc}\left(t\right)=\Delta {\Pi }_X\left(t\right)t\ge 0.\hfill & \end{array}$$

(A3j)

Then, we define the total osmo-oncotic pressure difference as

$$\begin{array}{ccc}& \Delta \Pi \left(t\right)=\Delta {\Pi }_{osm}\left(t\right)+\Delta {\Pi }_{onc}\left(t\right)\hfill & \\ & =\sum _{\alpha \in S_{\alpha }}{\sigma }_{\alpha }\left[F{z}_{\alpha }\langle c_{\alpha }\rangle \left(t\right){\psi }_m\left(t\right)+RT\langle c_{\alpha }\rangle \left(t\right)ln\left({\displaystyle \frac{c_{\alpha ,in}\left(t\right)}{c_{\alpha ,ex}\left(t\right)}}\right)\right]\hfill & \\ & +RT\sum _{\beta \in S_{\beta }}{\sigma }_{\beta }\Delta c_{\beta }\left(t\right)+{\sigma }_{X}RT\Delta c_X\left(t\right)t\ge 0.\hfill & \end{array}$$

(A3k)

Relation (A3k) is the generalization of Eq. (31) in [4].

Appendix A.2. Neutral solutes

Let us consider a carrier protein whose geometrical structure is the same as the pore ${\omega }_p$ illustrated in Figure 3 in the case of an AQP. To determine the normal molar flux density of solute $\beta$ inside ${\omega }_p$, we make the following assumptions.

Assumption 15.

We assume that:

1. 

the molar flux density of the neutral solute β has only the axial component $j_{\beta ,s}$;

2. 

for any $t\ge 0$, the quantity $j_{\beta ,s}$ is spatially constant inside ${\omega }_p$, i.e., $j_{\beta ,s}=j_{\beta ,n}\left(t\right)$, with $t\ge 0$ and $s\in [0,t_m]$.

Using Assumptions A15, Assumption 4-A3 and the fact that $v_{p,n}$ is spatially constant, into the constitutive equation (6b) gives the following second-order equation for the solute molar density inside the carrier protein channel

$$\begin{array}{ccc}& v_{p,n}\left(t\right){\displaystyle \frac{\partial c_{\beta }(s,t)}{\partial s}}-D_{\beta }{\displaystyle \frac{{\partial }^2{c}_{\beta }(s,t)}{\partial s^2}}=0t\ge 0,s\in [0,t_m],\hfill & \end{array}$$

(A4a)

whose solution is

$$\begin{array}{ccc}& c_{\beta }(s,t)=A\left(t\right)+B\left(t\right)exp\left({\displaystyle \frac{v_{p,n}\left(t\right)}{D_{\beta }}}s\right)t\ge 0,s\in [0,t_m],\hfill & \end{array}$$

(A4b)

$A=A\left(t\right)$ and $B=B\left(t\right)$ being time-dependent arbitrary functions. Inserting (A4b) into the constitutive equation (6b) yields

$$\begin{array}{ccc}& j_{\beta ,n}\left(t\right)=v_{p,n}\left(t\right)A\left(t\right)t\ge 0.\hfill & \end{array}$$

(A4c)

To determine $A\left(t\right)$ we need to enforce the boundary conditions:

$$\begin{array}{ccc}& c_{\beta }(0,t)=c_{\beta ,in}\left(t\right),c_{\beta }(t_m,t)=c_{\beta ,ex}\left(t\right)t\ge 0,\hfill & \end{array}$$

(A4d)

where $c_{\beta ,in}\left(t\right)$ and $c_{\beta ,ex}\left(t\right)$ are the intra and extracellular values of the neutral solute molar density for any $t\ge 0$. Using (A4d) into (A4b), we obtain the following expression of the normal molar flux density of $\beta$ inside the carrier protein channel

$$\begin{array}{ccc}& j_{\beta ,n}\left(t\right)=-P_{\beta }\left[c_{\beta ,ex}\left(t\right)\mathrm{Be}\left(X_{\beta }\left(t\right)\right)-c_{\beta ,in}\left(t\right)\mathrm{Be}(-X_{\beta }\left(t\right))\right]t\ge 0,\hfill & \end{array}$$

(A4e)

where $P_{\beta }:=D_{\beta }/t_m$ is the membrane permeability to the neutral solute $\beta$ (units: $\mathrm{m}{\mathrm{s}}^{-1}$),

$$\begin{array}{ccc}& \mathrm{Be}\left(W\right):={\displaystyle \frac{W}{e^W-1}}W\in \mathbb{R}\hfill & \end{array}$$

is the inverse of the Bernoulli function and

$$\begin{array}{ccc}& X_{\beta }\left(t\right):={\displaystyle \frac{v_{p,n}\left(t\right)}{P_{\beta }}}t\ge 0.\hfill & \end{array}$$

(A4f)

Appendix A.3. Charged solutes

We can proceed in the same way as in Appendix A.2 to obtain the following expression for the normal molar flux density of ion $\alpha \in S_{\alpha }$ inside an ionic channel

$$\begin{array}{ccc}& j_{\alpha ,n}^{edw}\left(t\right)=-P_{\alpha }\left[c_{\alpha ,ex}\left(t\right)\mathrm{Be}\left(X_{\alpha }\left(t\right)\right)-c_{\alpha ,in}\left(t\right)\mathrm{Be}(-X_{\alpha }\left(t\right))\right]t\ge 0,\hfill & \end{array}$$

(A5a)

where $P_{\alpha }:=D_{\alpha }/t_m$ is the membrane permeability to the charged solute $\alpha$ and

$$\begin{array}{ccc}& X_{\alpha }\left(t\right):={\displaystyle \frac{v_{p,n}\left(t\right)}{P_{\alpha }}}+z_{\alpha }{\displaystyle \frac{{\psi }_m\left(t\right)}{V_{th}}}t\ge 0.\hfill & \end{array}$$

(A5b)

Using the definition (8e) of the generalized drift velocity of ion $\alpha$ into (A5b) we find the spatial distribution of $c_{\alpha }$ inside the ionic channel

$$\begin{array}{ccc}& c_{\alpha }(s,t)=\tilde{A}\left(t\right)+\tilde{B}\left(t\right)exp\left({\displaystyle \frac{v_{\alpha ,n}\left(t\right)}{D_{\alpha }}}s\right)t\ge 0,s\in [0,t_m],\hfill & \end{array}$$

(A5c)

where $\tilde{A}=\tilde{A}\left(t\right)$ and $\tilde{B}=\tilde{B}\left(t\right)$ are time-dependent arbitrary functions, and $v_{\alpha ,n}\left(t\right)$ is the average normal component of ${\mathbf{v}}_{\alpha }={\mathbf{v}}_{\alpha }(\mathbf{x},t)$ over the cell surface, $\mathbf{x}\in \partial {\Omega }_t$ and $t\ge 0$. Replacing (A5c) into (A3e), we obtain

$$\begin{array}{ccc}& \langle c_{\alpha }\rangle \left(t\right)=c_{\alpha ,in}\left(t\right){\xi }_{in}\left(t\right)+c_{\alpha ,ex}\left(t\right){\xi }_{ex}\left(t\right)t\ge 0,\hfill & \end{array}$$

(A5d)

where:

$$\begin{array}{ccc}& {\xi }_{in}\left(t\right)={\displaystyle \frac{\mathrm{Be}\left(-X_{\alpha }\left(t\right)\right)-1}{X_{\alpha }\left(t\right)}}t\ge 0,\hfill & \end{array}$$

(A5e)

$$\begin{array}{ccc}& {\xi }_{ex}\left(t\right)={\displaystyle \frac{1-\mathrm{Be}\left(X_{\alpha }\left(t\right)\right)}{X_{\alpha }\left(t\right)}}t\ge 0.\hfill & \end{array}$$

(A5f)

Figure A1. Red solid line: plot of $\langle c_{\alpha }\rangle$ for $\tau =[0:20]$ and $X_{\alpha }\left(\tau \right)=\tau -10$, $\tau$ being a dimensionless time. The endpoint values of $c_{\alpha }$ are $c_{\alpha ,in}\left(\tau \right)=1\mathrm{mM}$ and $c_{\alpha ,ex}\left(t\right)=5\mathrm{mM}$ for every $\tau \in [0:20]$. Black dashed line: arithmetic average of $c_{\alpha }$.

Figure A1. Red solid line: plot of $\langle c_{\alpha }\rangle$ for $\tau =[0:20]$ and $X_{\alpha }\left(\tau \right)=\tau -10$, $\tau$ being a dimensionless time. The endpoint values of $c_{\alpha }$ are $c_{\alpha ,in}\left(\tau \right)=1\mathrm{mM}$ and $c_{\alpha ,ex}\left(t\right)=5\mathrm{mM}$ for every $\tau \in [0:20]$. Black dashed line: arithmetic average of $c_{\alpha }$.

We illustrate in Figure A1 the graphical representation of $\langle c_{\alpha }\rangle \left(\tau \right)$ (red solid line) compared to the arithmetic value $c_{\alpha }^{arith}\left(\tau \right)=(c_{\alpha ,in}\left(\tau \right)+c_{\alpha ,ex}\left(\tau \right))/2$ for $\tau =[0:20]$ (Matlab vector notation), $X_{\alpha }\left(\tau \right)=\tau -10$, $c_{\alpha ,in}\left(\tau \right)=1\mathrm{mM}$ and $c_{\alpha ,ex}\left(\tau \right)=5\mathrm{mM}$. We see that $\langle c_{\alpha }\rangle$ and $c_{\alpha }^{arith}$ are comparably close only if $|X_{\alpha }|\simeq 0$, which corresponds to a diffusive regime of transport inside the channel. Conversely, their distance increases for larger values of $|X_{\alpha }|$, which corresponds to an advection-dominated regime of transport inside the channel. On the basis of these considerations, the adoption of (A5d) instead of the customary choice $c_{\alpha }^{arith}$ (as done, for example, in [4]) is expected to warrant an accurate and robust model simulation.

Appendix B.1. Mathematical model of CA-mediated CO₂ hydrolysis

In this section we illustrate the mathematical modeling of the CA-enzyme-mediated ${\mathrm{CO}}_2$ hydrolysis. This process can be conveniently represented as the following two-step chemical reaction (see [36,37]):

$$\begin{array}{ccc}& \mathrm{STEP}1:{\mathrm{CO}}_2+{\mathrm{H}}_2\mathrm{O}\rightleftarrows {\mathrm{H}}_2{\mathrm{CO}}_3.\hfill & \end{array}$$

(A6a)

$$\begin{array}{ccc}& \mathrm{STEP}2:{\mathrm{H}}_2{\mathrm{CO}}_3\rightleftarrows {\mathrm{H}}^{+}+{\mathrm{HCO}}_3^{-}.\hfill & \end{array}$$

(A6b)

STEP 1 is the conversion of intracellular carbon dioxide into carbonic acid under mediation of the CA enzyme, which is a very fast catalyzer of the CO₂ hydration process (see [38,39]). STEP 1 is characterized by the quantities $k_{hydr}$ and $k_{dehydr}$ (units: ${\mathrm{s}}^{-1}$) depending on the molar density of the carbonic anhydrase enzyme (CA) and representing the rate constants of the forward and backward reactions in STEP 1, respectively.

Assumption 16.

Let $c_{{\mathrm{CO}}_{2\left(aq\right),in}}$ denote the molar density of intracellular CO₂ that is hydrated by water molecules. We assume that each ${\mathrm{CO}}_2$ molecule that is produced by mitochondria respiration is hydrated by a corresponding water molecule. This allows us to denote by $c_{{\mathrm{CO}}_{2,in}}$ the molar density of the intracellular $c_{{\mathrm{CO}}_{2\left(aq\right)}}$.

Using the Law of Conservation of Mass and Assumption 16, the net production rates in the mass balance equation for $\beta ={\mathrm{CO}}_2$ and $\beta ={\mathrm{H}}_2{\mathrm{CO}}_3$ have the following expressions:

$$\begin{array}{ccc}& {\mathcal{R}}_{{\mathrm{CO}}_2}\left(t\right)=k_{dehydr}{c}_{{\mathrm{H}}_2{\mathrm{CO}}_3,in}\left(t\right)-k_{hydr}{c}_{{\mathrm{CO}}_{2,in}}\left(t\right)t\ge 0,\hfill & \end{array}$$

(A7a)

$$\begin{array}{ccc}& {\mathcal{R}}_{{\mathrm{H}}_2{\mathrm{CO}}_3}\left(t\right)=-{\mathcal{R}}_{{\mathrm{CO}}_2}\left(t\right)t\ge 0.\hfill & \end{array}$$

(A7b)

The hydration and dehydration rate constants experimentally depend on the amount of CA that is present in the compartment where the hydration reaction takes place (see [40]). Following [36,37], in this article we use the following model for the hydration and dehydration rate constants:

$$\begin{array}{ccc}& k_{hydr}=A_{\mathrm{CA}}{k}_{hydr,ref},\hfill & \end{array}$$

(A7c)

$$\begin{array}{ccc}& k_{dehydr}=A_{\mathrm{CA}}{k}_{dehydr,ref},\hfill & \end{array}$$

(A7d)

where $A_{\mathrm{CA}}$ is a nonnegative given constant, whereas $k_{hydr,ref}=0.037{\mathrm{s}}^{-1}$ and $k_{dehydr,ref}=13.7{\mathrm{s}}^{-1}$ are experimental reference values measured at $T={25}^{\circ }\mathrm{C}$ reported in [41]. Taking $A_{\mathrm{CA}}>1$ is a way to represent the catalyzing effect of CA compared to the uncatalyzed reaction corresponding to $A_{\mathrm{CA}}=1$. In the numerical simulations illustrated in Section 3.2 we set $A_{\mathrm{CA}}=5$.

STEP 2 is the dissociation of intracellular carbonic acid into bicarbonate and protonated hydrogen, and is characterized by the quantities $k_{diss}$ (units: ${\mathrm{mM}}^{-1}{\mathrm{s}}^{-1}$) and $k_{assoc}$ (units: ${\mathrm{mM}}^{-2}{\mathrm{s}}^{-1}$), the rate constants of the forward and backward reactions in STEP 2, respectively. The dissociation reaction of H₂CO₃ is extremely rapid, so that the values of the forward and backward rate constants in STEP 2 are very large. In the numerical simulations illustrated in Section 3.2 we use the data of [36,37], and set $k_{diss}={10}^{16}{\mathrm{s}}^{-1}$ and $k_{assoc}=k_{diss}/K_{eq}$, where $K_{eq}=0.2804mM$ is the equilibrium constant of STEP 2. Using the Law of Conservation of Mass, we obtain the following expressions for the net production rates in the mass balance equations for $\alpha ={\mathrm{H}}^{+}$ and $\alpha ={\mathrm{HCO}}_3^{-}$:

$$\begin{array}{ccc}& {\mathcal{R}}_{\alpha }\left(t\right)=k_{diss}{c}_{{\mathrm{H}}_2{\mathrm{CO}}_3}\left(t\right)-k_{assoc}{c}_{{\mathrm{H}}^{+}}\left(t\right)c_{{\mathrm{HCO}}_3^{-}}\left(t\right)\alpha =\left\{{\mathrm{H}}^{+},{\mathrm{HCO}}_3^{-}\right\},t\ge 0.\hfill & \end{array}$$

(A8)

Assumption 17.

We assume that ions $\alpha \in \left\{{\mathrm{Na}}^{+},{\mathrm{K}}^{+},{\mathrm{Cl}}^{-}\right\}$ are non-reacting (see [42], Section 8.3.3). Therefore, we have ${\mathcal{R}}_{\alpha }\left(t\right)=0$, $t\ge 0$.

Appendix B.2. Net production rate in cell volume regulation

The mechanisms which govern intracellular water production/consumption are object of considerable debate and investigation because of their importance in cell life and survival. One example is provided by the process of normotonic cell shrinkage which is the major hallmark of cellular apoptosis [43]. Another example is provided by the TCA cycle (Krebs cycle) in cellular respiration, whose products of metabolism of fuels are ATP, CO₂, and the so-called “metabolic” water [44].

The difficulty of accessing and sampling the contents of intact cells makes the study of the intracellular fluid environment, and, more generally, of body fluid content, a challenging problem. Specific approaches to accurately detect metabolic water content as a result of intracellular metabolic activity have been proposed in [45] and [46].

In this article, we focus on the role of net production and consumption of intracellular water ${\mathcal{R}}_w$ in driving the motion of the cell surface. Our proposed model is

$$\begin{array}{ccc}& {\mathcal{R}}_w\left(t\right)=k_{w,p}-k_{w,c}{\displaystyle \frac{\mathcal{V}\left(t\right)}{{\mathcal{V}}_0}}t\ge 0,\hfill & \end{array}$$

(A9)

where $k_{w,p}$ and $k_{w,c}$ are given constants (units: ${\mathrm{s}}^{-1}$) and ${\mathcal{V}}_0$ is the value of cell volume in resting conditions. In the simulations illustrated in Section 3.2, we set $k_{w,p}=k_{w,c}=1{\mathrm{s}}^{-1}$.

Appendix B.3. Na + /K + ATPase

The mathematical model of the sodium-potassium pump (Na⁺/K⁺ ATPase) adopted in this article follows the idea proposed in [3]. The molar flux densities for sodium and potassium are:

$$\begin{array}{ccc}& j_{\mathrm{Na},n}^a\left(t\right)=j_{pump}\left(t\right)t\ge 0,\hfill & \end{array}$$

(A10a)

$$\begin{array}{ccc}& j_{\mathrm{K},n}^a\left(t\right)=-{\displaystyle \frac{2}{3}}{j}_{pump}\left(t\right)t\ge 0,\hfill & \end{array}$$

(A10b)

where:

$$\begin{array}{cc}& j_{pump}\left(t\right)=j_{pump}^{MAX}{\left({\displaystyle \frac{c_{{\mathrm{Na}}_{in}}\left(t\right)}{c_{{\mathrm{Na}}_{in}}\left(t\right)+c_{\mathrm{Na},1/2}}}\right)}^3{\left({\displaystyle \frac{c_{{\mathrm{K}}_{ex}}\left(t\right)}{c_{{\mathrm{K}}_{ex}}\left(t\right)+c_{\mathrm{K},1/2}}}\right)}^{2}t\ge 0,\hfill \end{array}$$

(A10c)

$$\begin{array}{ccc}& j_{pump}^{MAX}:=c_{ATP}{v}_{pump}.\hfill & \end{array}$$

(A10d)

The quantity $j_{pump}^{MAX}$ is the maximum pump molar flux density (units: $\mathrm{mM}\mathrm{m}{\mathrm{s}}^{-1}$, $c_{ATP}$ is the intracellular molar density of ATP (units: $\mathrm{mM}$) and $v_{pump}$ is the ion transfer velocity of the pump (units: $\mathrm{m}{\mathrm{s}}^{-1}$). The quantities $c_{ATP}$ and $v_{pump}$ are defined as:

$$\begin{array}{ccc}& c_{ATP}=M_{ATP}{c}_{ATP,ref},\hfill & \end{array}$$

(A10e)

$$\begin{array}{ccc}& v_{pump}=r_{turn}{t}_M,\hfill & \end{array}$$

(A10f)

where $c_{ATP,ref}$ is the reference value of ATP molar density (units: $\mathrm{mM}$), $r_{turn}$ is pump turnover rate (units: ${\mathrm{s}}^{-1}$) and $M_{ATP}$ is a nonnegative given constant. The quantities $c_{\mathrm{Na},1/2}$ and $c_{\mathrm{K},1/2}$ are the Michaelis constants of the pump model (units: $\mathrm{mM}$).

In the simulations illustrated in Section 3.2 we set $c_{ATP,ref}={10}^{-2}mM$, $M_{ATP}=2$, $r_{turn}={10}^7{\mathrm{s}}^{-1}$, $c_{\mathrm{Na},1/2}=1.3\mathrm{mM}$ and $c_{\mathrm{K},1/2}=0.14\mathrm{mM}$.