Multiscale Multiphysics Modeling of Aqueous Humor Dynamics in the Human Eye

1. Introduction

The healthy turnover in production, flow and drainage of aqueous humor (AH) in the human eye is vital to ocular homeostasis and vision preservation. AH is a watery fluid continuously produced by the ciliary processes and slowly moving across the posterior and anterior chamber to be eventually drained out for successive recirculation. Experimental data indicate that in normal conditions the time to complete an entire recirculation cycle is about 90 minutes [1]. This allows us to mathematically assume AH dynamics to be a stationary process and safely neglect the dependence on the time variable. Homeostatic conditions (usually referred to as "baseline conditions") correspond to the situation where the volumetric flow rates (VFR) of produced, circulating and drained AH are equal to each other. In such conditions, measured data reported in [1] indicate that the baseline values of the AH VFR circulating in the eye and pressure in the eye (referred to as intraocular pressure, IOP) are $Q_{\mathrm{AH},\mathrm{b}}=2.75\mu {\mathrm{Lmin}}^{-1}$ and $p_{\mathrm{IO},\mathrm{b}}=15\mathrm{mmHg}$, respectively.

The above described healthy situation undergoes a change under pathological conditions where IOP may depart from a range associated with "normal conditions". A significant example is primary open angle glaucoma (POAG), a disease of the optic nerve representing the second leading cause of blindness worldwide [2], for which the range is usually between 10 and 21 mmHg (see [3]). Another example is thyroid eye disease (TED), an autoimmune ocular pathology characterized by swelling and bulging of the eye globe due to inflammation and remodeling of the extraocular muscles [4], for which the baseline range for IOP is between 23 and 27 mmHg (see [5]). POAG and TED are multifactorial diseases whose causes have a different nature and act at very different spatial scales thereby making them unfavorable to disentanglement and often even unaccessible to measurement. For these reasons, availing of a mechanism-driven mathematical model, possibly integrated with artificial intelligence-based (AI) data-driven approaches, may result useful to shed light on the complex etiology of the considered disease (see [6,7]).

Mechanism-driven models can be classified into two main categories: compartment models (referred to as 0D models) and continuum models (referred to as dD models, $d=1,2,3$). The construction of a 0D model usually leverages the analogy between fluid and electric variables (see [8]) in which compartments are the nodes of an equivalent electric circuit and nodes are connected between each other by a resistor and to ground by a voltage or current source (see [9,10,11,12]). 0D models are computationally efficient but have the intrinsic limitation of neglecting the spatial variation of problem unknowns. This limitation is removed by the adoption of dD models which are constituted by systems of differential equations describing the multiphysical dynamics of AH, including fluid, mechanical and electrochemical phenomena (see [13,14,15,16]). dD models are characterized by a significant information content but have the limitation of requiring large computational resources to effectively deal with the wide multiscale nature of AH dynamics.

The objective of this study was to develop a mathematical and numerical model of AH dynamics capable of describing the physiology of the problem at each spatial scale while keeping the cost of the simulation at an affordable level. To this purpose, we introduced a stationary representation by compartments consisting of three integrated modules (M): M1 for AH production, M2 for AH passive flow and M3 for AH drainage. M1 takes into account the cellular composition of the production compartment, M2 and M3 the hydrodynamic characteristics and porous structure of passive flow and drainage compartments, respectively. M1 is the result of a 0D reduction of the velocity-extended Poisson-Nernst-Planck model (see [8], Chapter 13) and simulates the transfer of neutral and charged solutes and fluid movement across the cellular structure of the ciliary epithelium (CE). M2 is the electric equivalent representation of Poiseuille flow across the series of two linear hydraulic resistors (see [13] and [8], Chapter 10). M3 is a 0D reduction of the Darcy equations for a porous medium (see [17]) and simulates AH flow across a parallel connection of a nonlinear and linear resistors.

The proposed model was numerically implemented to comprehend how the findings of each module may help to explain the onset and progression of an ocular disease and to evaluate their impact on patient’s health by contributing to develop specifically targeted therapeutical approaches. In particular, we used the integration of M1, M2 and M3 to give an answer to open questions about the role of (1) sodium channels in AH production and (2) episcleral vein pressure in AH drainage. Model predictions suggest that (1) the expression of sodium channels in the CE is essential for AH to flow from the ciliary processes into the posterior chamber; and (2) an increase of the episcleral vein pressure reduces AH drainage, potentially explaining the development of secondary open-angle glaucoma.

The structure of the paper is as follows. Section 2 illustrates the basic elements of the physiology of AH dynamics. Section 3 introduces the mathematical model of a biological semipermeable membrane. Section 4 illustrates the electric equivalent scheme of AH dynamics and its three-module composition. The multiscale equations of the model are described in Section 5 while the cellular mechanisms taking place in the CE are described in Section 6. Section 7 contains a description of the solution map that is adopted to numerically solve the model equations. Simulation results are illustrated in Section 8 and analyzed in Section 9. Concluding remarks are drawn in Section 10 while Appendix A contains the numerical values of model parameters used in the computations.

2. Physiology of AH Dynamics

Figure 1 is a multiscale view of the human eye focusing on the regions that are involved in the process of AH dynamics. The leftmost panel is a cross-section of the eye in the sagittal plane. The panel in the center of the figure is an histology of the anterior segment of the eye and shows the frontal portion of the eye occupied by the aqueous humor (anterior chamber, AC), the cornea, the iris, the ciliary body (CB) and the trabecular meshwork (TM). The rightmost panel is a zoomed view of the ciliary body and shows a ciliary process, which is the region of the CB where AH is produced. The ciliary process includes the ciliary capillaries (CC, not shown here), the ciliary stroma (CS) and the cellular double layer constituting the ciliary epithelium (CE). This latter is formed by a couplet of cells facing each other from the apical side, the pigmented epithelial (PE) and nonpigmented epithelial (NPE) cells, respectively. The region indicated with "Aqueous humor" and facing from below the CE is a subset of the posterior chamber (PC). In normal conditions, AH flows out of the CCs throughout fenestrations located on the capillary surface; then, it traverses the CS and the CE, which filters out the solutes with large molecular weight (including blood particles); eventually, the produced aqueous fluid flows along the PC until it reaches the AC. In its final course, AH exits the AC for further recirculation into the CCs by means of the TM and uveoscleral (UV) outflow pathways. In the case of the TM pathway, AH percolates through the multi-layered porous structure of the TM tissue from which it enters the Schlemm’s canal and eventually drains into the episcleral vein (EV). In the case of the UV pathway, AH diffuses through the ciliary muscle bundles into the suprachoroidal space (SCS) from which it drains into the choroidal vessels (see [19]).

3. Fluid and Solute Transport Across a Semipermeable Membrane

The description in Section 2 shows that AH travels across six biological semipermeable membranes (SM) during its flow from the CB to the TM and UV outflow pathways. The first membrane is the wall of the fenestrated CCs. The second and third membranes are the basement membranes separating the CE from the CS and PC, respectively. The fourth membrane is the pupil separating the PC from the AC. The fifth and sixth membranes are the porous structures of the TM and ciliary muscle, respectively. In the next sections we introduce a mathematical representation of AH flow across a biological membrane and its interpretation using the analogy between fluid and electric variables.

3.1. Starling’s Law of Ultrafiltration

Let a and b denote two compartments separated by a SM ${\mathcal{M}}_{ab}$ and U a physical variable with values $U_a$ and $U_b$ in a and b, respectively, with $\Delta U:=U_a-U_b$. Fluid and solute transport across ${\mathcal{M}}_{ab}$ can be mathematically described using Starling’s law of ultrafiltration (see [20] and [21])

$$\begin{array}{ccc}& Q_{\mathrm{ab}}=C_{\mathrm{ab}}\Delta \mathcal{P}.\hfill & \hfill \end{array}$$

(1a)

The quantities $Q_{\mathrm{ab}}$ and $C_{\mathrm{ab}}$ are the AH VFR across ${\mathcal{M}}_{ab}$ and the facility of ${\mathcal{M}}_{ab}$, respectively, whereas

$$\begin{array}{ccc}& \mathcal{P}:=p-\Pi \hfill & \hfill \end{array}$$

(1b)

is the total pressure, given by the difference between the hydraulic pressure p and the electro-osmo-oncotic (EOO) pressure $\Pi$, which is defined in the next section.

3.2. Electro-Osmo-Oncotic Pressure

AH is a homogeneous mixture constituted by water, charged solutes (ions) and neutral solutes (see [8], Chapters 13, 26). We denote by $S_{\alpha }$ and $S_{\beta }$ the sets containing the labels of ions and neutral solutes, respectively; for each $\alpha \in S_{\alpha }$, $c_{\alpha }$ and $z_{\alpha }$ are the molar density and charge number of ion $\alpha$, respectively; for each $\beta \in S_{\beta }$, $c_{\beta }$ is the molar density of neutral solute $\beta$; finally, $c_{\mathrm{X}}$ is the molar density of the impermeant protein charge which ensures the electroneutrality of the AH solution. For each $\alpha \in S_{\alpha }$ and $\beta \in S_{\beta }$, the quantities ${\sigma }_{\alpha }$ and ${\sigma }_{\beta }$ are the reflection coefficients of ion $\alpha$ and neutral solute $\beta$ associated with the SM, respectively. Both ${\sigma }_{\alpha }$ and ${\sigma }_{\beta }$ range in $[0,1]$, with ${\sigma }_{\alpha }$ (${\sigma }_{\beta }$) being close to 0 or 1 if the SM is highly permeable or impermeant to the considered ion (solute), respectively, implying that ${\sigma }_{\mathrm{X}}=1$. The EOO pressure difference between compartments a and b is defined as follows(see [22])

$$\begin{array}{ccc}& \Delta \Pi =F\sum _{\alpha \in S_{\alpha }}{\sigma }_{\alpha }{z}_{\alpha }\langle c_{\alpha }\rangle \Delta {\phi }_{\alpha }^{\mathrm{ec}}+RT\sum _{\beta \in S_{\beta }}{\sigma }_{\beta }\Delta c_{\beta }+RT{\sigma }_{\mathrm{X}}\Delta c_{\mathrm{X}},\hfill & \hfill \end{array}$$

(2a)

where F and R are Faraday and gas constants, respectively, T is the absolute temperature of the system, $\langle c_{\alpha }\rangle =(c_{\alpha ,a}+c_{\alpha ,b})/2$ and for each $\alpha \in S_{\alpha }$ the quantity

$$\begin{array}{ccc}& {\phi }_{\alpha }^{\mathrm{ec}}=\psi +{\displaystyle \frac{RT}{F{z}_{\alpha }}}ln\left({\displaystyle \frac{c_{\alpha }}{c_{\mathrm{ref}}}}\right)\hfill & \hfill \end{array}$$

(2b)

is the electrochemical potential of ion $\alpha$, with $\psi$ and $c_{\mathrm{ref}}$ denoting the electric potential and the ion reference molar density, respectively.

3.3. The Starling Resistor

Leveraging the analogy between electric and fluid variables (see [8], Chapter 15), we illustrate in Figure 2 the pictorial representation of equations (1). The left panel of Figure 2 is a scheme of the physical system constituted by a SM ${\mathcal{M}}_{ab}$ separating two compartments a and b containing a homogeneous solution constituted by water, ions and neutral solutes. The right panel of Figure 2 is a scheme of the "Starling resistor", constituted by the electric conductance $C_{ab}$ connecting two nodes of a circuit, denoted by a and b. Based on this electric equivalent scheme, we see that equation (1a) can be regarded as a generalized Ohm’s law in which the electric current $Q_{ab}$ flowing across a resistor is proportional to the voltage $V:=\Delta \mathcal{P}$. Interestingly, using the definition (1b) we see that the AH VFR is the sum of two terms, whose difference expresses the physical fact that AH water flow due to the EOO pressure difference has an opposite direction with respect to AH water flow due to the hydraulic pressure difference.

4. Electric Equivalent Scheme of AH Dynamics

Figure 3 is the electric equivalent representation of the anterior segment of the eye and of the three main processes that characterize AH dynamics. M1, M2 and M3 are AH production, passive flow and drainage modules, respectively. The black arrows indicate the communication between modules and the circles are the nodes of the electric equivalent circuit. The red circles are the nodes where the dependent variables of the model are given quantities while the numbered yellow circles are the nodes of the circuit where the dependent variables of the model must be computed.

Assumption A1

(Input data for the EOO). The EOO ${\overline{\Pi }}_{\mathrm{CC}}$ is given and equal to ${\Pi }_{\mathrm{bp}}$, where ${\Pi }_{\mathrm{bp}}$ is computed using (Section 3.2) in which the values of electric potential and solute molar densities are measured in the blood plasma. The EOOs ${\overline{\Pi }}_{\mathrm{EV}}$ and ${\overline{\Pi }}_{\mathrm{EV}}$ are given and equal to ${\Pi }_{\mathrm{aq}}$, where ${\Pi }_{\mathrm{aq}}$ is computed using (Section 3.2) in which the values of electric potential and solute molar densities are measured in the aqueous side. The EOO at the nodes ②, ③, ④ and ⑤ is equal to ${\Pi }_{\mathrm{aq}}$.

The blue arrow labeled "BF" represents blood flow transporting the "metabolic fuel" necessary for the production of the aqueous watery fluid indicated by the green arrow (see [10]). The block labeled "CE" represents the cellular scale mechanisms occurring in the CE. Assumption A1 implies that passive AH flow throughout the PC and AC, as well as AH flow throughout the UV outflow pathway, is described by the fluid analog of Ohm’s law which is obtained by setting $\Delta \Pi =0$ in (1). $C_{\mathrm{PC}}$, $C_{\mathrm{AC}}$, $C_{\mathrm{TM}}$ and $C_{\mathrm{UV}}$ are the hydraulic conductances (referred to as "facilities") of the PC, AC and of the TM and UV outflow pathways, respectively. The black arrow across $C_{\mathrm{TM}}$ indicates that the TM facility depends nonlinearly on the pressure difference $p_{\mathrm{AC}}-{\overline{p}}_{\mathrm{EV}}$. $Q_{\mathrm{CC}}$ is the AH VFR that filtrates across the fenestrated CCs into the CS whereas $Q_{\mathrm{CS}}$ is the AH VFR that travels between the CS and reaches the CE. $Q_{\mathrm{prod}}$ is the AH VFR that is produced by module M1 whereas $Q_{\mathrm{TM}}$ and $Q_{\mathrm{UV}}$ are the AH VFRs drained out by TM and UV outflow pathways, respectively.

5. Multiscale Mathematical Model of AH Dynamics

In this section we illustrate the equations that constitute the mathematical model of AH dynamics at the macroscale and microscale levels.

5.1. Macroscale Equations

Referring to the circuit illustrated in Figure 3, the application of Kirchhoff current law (KCL) at node ①, the continuity of current across block CE, and the application of KCL at nodes ②, ③, ④ and ⑤ give rise to the following set of algebraic equations:

$$\begin{array}{ccc}& -C_{\mathrm{wall}}\left[{\overline{\mathcal{P}}}_{\mathrm{CC}}-{\mathcal{P}}_{\mathrm{CS}}\right]+Q_{\mathrm{CS}}=0,\hfill & \hfill \end{array}$$

(3a)

$$\begin{array}{ccc}& -Q_{\mathrm{CS}}+Q_{\mathrm{prod}}=0,\hfill & \hfill \end{array}$$

(3b)

$$\begin{array}{ccc}& -Q_{\mathrm{prod}}+C_{\mathrm{PC}}\left[{\mathcal{P}}_{\mathrm{PC}}-{\mathcal{P}}_{\mathrm{IO}}\right]=0,\hfill & \hfill \end{array}$$

(3c)

$$\begin{array}{ccc}& -C_{\mathrm{PC}}\left[{\mathcal{P}}_{\mathrm{PC}}-{\mathcal{P}}_{\mathrm{IO}}\right]+C_{\mathrm{AC}}\left[{\mathcal{P}}_{\mathrm{IO}}-{\mathcal{P}}_{\mathrm{AC}}\right]=0,\hfill & \hfill \end{array}$$

(3d)

$$\begin{array}{ccc}& -C_{\mathrm{AC}}\left[{\mathcal{P}}_{\mathrm{IO}}-{\mathcal{P}}_{\mathrm{AC}}\right]+Q_{\mathrm{prod}}=0,\hfill & \hfill \end{array}$$

(3e)

$$\begin{array}{ccc}& -Q_{\mathrm{prod}}+Q_{\mathrm{TM}}+C_{\mathrm{UV}}\left[{\mathcal{P}}_{\mathrm{AC}}-{\overline{\mathcal{P}}}_{\mathrm{SCS}}\right]=0.\hfill & \hfill \end{array}$$

(3f)

System (3) is the macroscale model of AH dynamics and consists of six algebraic equations for seven dependent variables, comprising three VFRs: $Q_{\mathrm{CS}}$, $Q_{\mathrm{prod}}$ and $Q_{\mathrm{TM}}$, and four pressures: ${\mathcal{P}}_{\mathrm{CS}}$, ${\mathcal{P}}_{\mathrm{PC}}$, ${\mathcal{P}}_{\mathrm{IO}}$ and ${\mathcal{P}}_{\mathrm{AC}}$.

5.2. The Unit Cell of the CE

Figure 4 shows the compartment subdivision of the "unit cell" whose periodic repetition represents the CE. White circles indicate the number of the compartment whereas white rectangles indicate the number of the interface separating two neighbouring compartments. The compartment numbered ① is the CS. The compartments numbered ② and ④ and indicated by "1/2 PE cell" and "1/2 NPE cell", respectively, correspond to one half of each PE and NPE cell. The compartments numbered ③ and ⑤ and indicated by "1/2 PE pc" and "1/2 NPE pc", respectively, correspond to one half of the paracellular (pc) spaces separating couplets of PE and NPE cells. The compartment numbered ⑥ is the PC. The rectangle in purple color is the BM and has a thickness equal to $t_{\mathrm{BM}}$. The rectangular boundary in light green color surrounding compartments numbered ② and ④ is the BLM and has a thickness equal to $t_{\mathrm{BL}}$. The rectangle in gold color is the GJ separating the PE and NPE cells whereas the rectangle with a porous pattern is the TJ separating the paracellular spaces ③ and ⑤. Since the CE is a "leaky epithelium" (see [23,24]) the TJ allows the passage of water, ions and low molecular weight (MW) molecules and does not allow the passage of large MW molecules. This is the reason why the AH watery fluid has a low content of proteins. The quantities a, b and c are the length of the sides of the various compartments.

5.3. Micro- to Macroscale Model of the CE

The input to the unit cell is $Q_{\mathrm{CS}}$ and the output from the unit cell is $Q_{\mathrm{prod}}$. Both quantities depend on the compartment distribution of the electric potential and of each physical constituent of the metabolic fuel of AH dynamics, namely, moving charged solutes (ions), fixed charged solutes (electroneutralizing charged proteins) and moving neutral solutes (low to medium-high molecular weight molecules). To couple the microscale representation of the CE with $Q_{\mathrm{CS}}$ and the micro- to macroscale models for $Q_{\mathrm{prod}}$ described in Section 5.4, respectively, we write the mass balance equations for water and solutes in each compartment of the unit cell (see [8], Chapter 7). Denoting by u any of the problem physical variables, the stationary mass balance of u is described by the following partial differential equation (PDE)

$$\begin{array}{ccc}& \nabla ·\mathbf{F}=\mathcal{G}-\mathcal{C},\hfill & \hfill \end{array}$$

(4a)

The quantity $\mathbf{F}$ is the flux density of u whereas $\mathcal{G}$ and $\mathcal{C}$ are the production and consumption rates of u (possibly depending on u). Denoting by ${\Omega }_i$, $i=1,\dots ,5$, the i-th compartment in the scheme of Figure 4, integrating both sides of (4a) over ${\Omega }_i$ and applying the divergence theorem to the left-hand side yields the following mass balance equation in integral form

$$\begin{array}{ccc}& {\int }_{{\Sigma }_i}\mathbf{F}·\mathbf{n}d\Sigma ={\int }_{{\Omega }_i}(\mathcal{G}-\mathcal{C})d\Omega i=1,\dots ,5,\hfill & \hfill \end{array}$$

(4b)

The left-hand side of (4b) is the net flux of u across the boundary surface ${\Sigma }_i$ whereas the right-hand side of (4b) is the net production rate of u inside the volume ${\Omega }_i$. Let us denote by ${\mathcal{A}}_i$ the index set of the compartments adjacent to ${\Omega }_i$ and define ${\Sigma }_{ij}={\Sigma }_i\cap {\Sigma }_j$. Let us replace $\mathbf{F}·\mathbf{n}$ by its value $F_{ij}$ at the centroid of ${\Sigma }_{ij}$ and assume that $F_{ij}$ depends on $u_i$, $u_j$ and on other quantities, with $u_i$, $u_j$ denoting the values of u at the centroids of ${\Omega }_i$ and ${\Omega }_j$, respectively. Finally, let us replace $\mathcal{G}$ and $\mathcal{C}$ by their values ${\mathcal{G}}_i$ and ${\mathcal{C}}_i$ at the centroid of ${\Omega }_i$, with the assumption that ${\mathcal{G}}_i$ and ${\mathcal{C}}_i$ may depend on $u_i$. Then, the integral mass balance equation (4b) becomes

$$\begin{array}{ccc}& \sum _{j\in {\mathcal{A}}_i}{F}_{ij}(u_i,u_j)|{\Sigma }_{ij}|=\left({\mathcal{G}}_i\left(u_i\right)-{\mathcal{C}}_i\left(u_i\right)\right)\left|{\Omega }_i\right|i=1,\dots ,5,\hfill & \hfill \end{array}$$

(4c)

where $|{\Sigma }_{ij}|$ and $|{\Omega }_i|$ denote the area of ${\Sigma }_{ij}$ and the volume of ${\Omega }_i$, respectively. For given values of u in the CC and PC, equations (4c) represent a nonlinear system of algebraic equations for the numerical approximation $\mathbf{U}={\left[u_1,u_2,u_3,u_4,u_5\right]}^T$ of the mean values of u in each compartment ${\Omega }_i$, $i=1,\dots ,5$.

5.4. Micro- to Macroscale Model of $Q_{\mathrm{prod}}$

Assumption A2

(Total pressure in the PC). In this section, ${\mathcal{P}}_{\mathrm{PC}}$, representing the total pressure in the PC, is a given quantity.

Once the equations (4c) are solved for $i=1,\dots ,5$, the AH VFR produced by the CE can be computed as

$$\begin{array}{ccc}& Q_{\mathrm{prod}}=2N_{\mathrm{cells}}\left[Q_{\mathrm{prod}}^{\mathrm{NPE}}+Q_{\mathrm{prod}}^{\mathrm{pc}}\right].\hfill & \hfill \end{array}$$

(5a)

The quantity $N_{\mathrm{cells}}$ denotes the number of NPE cells contained in the CE while the micro- to macroscale model for the AH VFR produced by the NPE cell, $Q_{\mathrm{prod}}^{\mathrm{NPE}}$, and by the paracellular space neighbouring the NPE cell, $Q_{\mathrm{prod}}^{\mathrm{pc}}$, is represented by Starling’s law as in (1a):

$$\begin{array}{ccc}& Q_{\mathrm{prod}}^{\mathrm{NPE}}=C_8\Delta {\mathcal{P}}_8=C_8\left[{\mathcal{P}}_{④}-{\mathcal{P}}_{\mathrm{PC}}\right],\hfill & \hfill \\ & Q_{\mathrm{prod}}^{\mathrm{pc}}=C_9\Delta {\mathcal{P}}_9=C_9\left[{\mathcal{P}}_{⑤}-{\mathcal{P}}_{\mathrm{PC}}\right],\hfill & \hfill \end{array}$$

(5b)

where $C_8$ and $C_9$ are the facility of the membranes separating the NPE and the paracellular space neighbouring the NPE cell from the PC, respectively.

5.5. Micro- to Macroscale model of $C_{\mathrm{TM}}$ and $Q_{\mathrm{TM}}$

Following [25], we assume that at the microscale level the TM can be described as a porous material of thickness $t_{\mathrm{TM}}$ comprised of a solid and fluid constituents, in which the porous structure consists of a collection of $N_p$ parallel circular, cylindrical pores of diameter $d_p$. The TM hydraulic resistance $R_{\mathrm{TM},0}$ can be related to the TM permeability ${\kappa }_{\mathrm{TM}}$ by Darcy’s law as

$$\begin{array}{ccc}& R_{\mathrm{TM},0}={\displaystyle \frac{{\mu }_{\mathrm{fl}}{t}_{\mathrm{TM}}}{{\kappa }_{\mathrm{TM}}{A}_{\mathrm{TM}}}},\hfill & \hfill \end{array}$$

(6a)

where ${\mu }_{\mathrm{fl}}$ is the dynamic viscosity of the fluid, $A_{\mathrm{TM}}$ is the total surface of the TM in contact with the AC and

$$\begin{array}{ccc}& {\kappa }_{\mathrm{TM}}={\displaystyle \frac{{\Phi }_{\mathrm{TM}}{d}_p^2}{16c_K}},\hfill & \hfill \end{array}$$

(6b)

having denoted by ${\Phi }_{\mathrm{TM}}$ and $c_K$ the porosity of the TM and the Kozeny constant, respectively. Experimental data shown in [26] suggest that the hydraulic resistance exerted by the TM is not constant but increases linearly with the total pressure difference $X_{\mathrm{TM}}:={\mathcal{P}}_{\mathrm{AC}}-{\overline{\mathcal{P}}}_{\mathrm{EV}}$ according to the following relation

$$\begin{array}{ccc}& R_{\mathrm{TM}}\left(X_{\mathrm{TM}}\right)=R_{\mathrm{TM},0}\left[1+{\alpha }_Q{X}_{\mathrm{TM}}\right],\hfill & \hfill \end{array}$$

(6c)

where ${\alpha }_Q$ is a nonnegative parameter representing the obstruction coefficient. Taking the reciprocal of (6c) gives the following mathematical model for the TM outflow facility

$$\begin{array}{ccc}& C_{\mathrm{TM}}\left(X_{\mathrm{TM}}\right)={\displaystyle \frac{C_{\mathrm{TM},0}}{1+{\alpha }_Q{X}_{\mathrm{TM}}}},\hfill & \hfill \end{array}$$

(6d)

where $C_{\mathrm{TM},0}:={\left(R_{\mathrm{TM},0}\right)}^{-1}$. The element with the transverse black arrow in module M3 of Figure 3 is the equivalent electric representation of the nonlinear TM outflow facility (6d). The micro- to macroscale model of the AH VFR throughout the TM is given by the following nonlinear Starling’s law of ultrafiltration

$$\begin{array}{ccc}& Q_{\mathrm{TM}}\left({\mathcal{P}}_{\mathrm{AC}}\right)=C_{\mathrm{TM},0}{\displaystyle \frac{\left[{\mathcal{P}}_{\mathrm{AC}}-{\overline{\mathcal{P}}}_{\mathrm{EV}}\right]}{1+{\alpha }_Q\left[{\mathcal{P}}_{\mathrm{AC}}-{\overline{\mathcal{P}}}_{\mathrm{EV}}\right]}}.\hfill & \hfill \end{array}$$

(6e)

The nonlinear Ohm’s law (6e) reduces to the linear Ohm’s law by setting ${\alpha }_Q=0$.

6. Cellular Mechanisms in AH Production

Figure 5 provides a schematic description of the membrane proteins that allow the exchange of neutral solutes, charged solutes (ions) and water between neighbouring compartments and of the main metabolic activities and chemical reactions occurring in the cytoplasm of the PE and NPE cells. The scheme is an extension of that proposed in [10] and [29] and includes two types of regulatory mechanisms of ion transport across a semipermeable membrane:

• passive transport mechanisms;
• active transport mechanisms.

Passive transport is mediated by ion channels (magenta color in Figure 5) and is directed along the electrochemical gradient across the membrane. Active transport is directed against the electrochemical gradient across the membrane and is divided into primary and secondary depending on the type of energy source and the way adopted by the transporter protein to exploit it. Primary active transport is mediated by the Na+/K+ ATPase (also referred to as ionic pump, dark grey color in Figure 5) which directly exploits the ATP produced by the mitochondria (yellow color in Figure 5). Secondary active transport is mediated by the co-transporters (white color circles and ellipse in Figure 5) which indirectly exploit the potential energy source due to the transmembrane concentration gradient of another ion or by the action of specific enzymes such as the Carbonic Anhydrase (CA). Intracellular metabolism of mitochondria regulates the function of cellular respiration. Mitochondria activity is driven by the absorption of glucose and produces carbon dioxide (CO₂) and ATP. CO₂ is rapidly eliminated from the cytoplasm by the hydrolysis reaction catalyzed by the CA. ATP is used by the Na⁺/K⁺ ATPase to extrude sodium ions from the cell and intrude potassium ions into the cell. Water transport across the various membranes in the scheme of Figure 5 follows two different pathways [30]:

• passive transport of water molecules (diffusion) across the lipid bilayer of the plasma membrane and the extracellular matrix constituents of the basement membrane;
• active transport of water molecules across the pores of specialized membrane proteins (aquaporins, AQP) driven by hydraulic-osmotic pressure gradients (facilitated diffusion).

AQPs (green color in Figure 5) are solely expressed on the BLM of the NPE cells which are, as a consequence, responsible of the large majority of AH production by the ciliary processes (see [24] and [27]).

7. Solution Map

The computation of the nodal pressures and AH VFRs represented in the scheme of Figure 3 requires the solution of the fully coupled multiscale model of AH dynamics (FCMM), consisting of the union of the macroscale equations (3) with the micro- to macroscale equations (5) and (6). The complexity of solving the FCMM lies in the fact that the micro- to macroscale equations depend on the solution at the cellular level of the system of microscale equations (4c) for each constituent of the AH watery fluid (water, moving neutral and charged solutes, immobile charged solutes and electric potential). To end up with a numerical procedure that is manageable for research and clinical applications, we propose in the present work a solution map (or functional iteration, see [31], [8]), which consists of transforming the task of solving the FCMM into the successive solution of problems of lower dimension until convergence. The solution map is illustrated in Algorithm 1 where $\mathbf{X}={\left[{\mathcal{P}}_{\mathrm{CS}},{\mathcal{P}}_{\mathrm{PC}},{\mathcal{P}}_{\mathrm{IO}},{\mathcal{P}}_{\mathrm{AC}}\right]}^T$, ${\displaystyle {\parallel \mathbf{W}\parallel }_{\infty }:=\underset{i=1,\dots ,n}{max}\left|W_i\right|}$ for any $\mathbf{W}\in {\mathbb{R}}^n$, and the superscript ${\left(\right)}^{*}$ identifies the approximate solution of the FCMM at convergence. Notice that for $k>0$, the initial guess ${\mathbf{X}}^{\left(0\right)}$ for the $k+1$-th iteration of the solution map is the linear interpolation of the last computed iterate ${\mathbf{X}}^{(k+1)}$ and the previous iterate ${\mathbf{X}}^{\left(k\right)}$ by means of a relaxation parameter $\theta$ to be chosen in the interval $(0,1]$.

Algorithm 1:Solution map for the FCMM

1: Provide an initial datum ${\mathbf{X}}^{\left(0\right)}$ and a tolerance tol. 2: Set $k=0$. 3: Set ${\mathcal{P}}_{\mathrm{PC}}={\mathcal{P}}_{\mathrm{PC}}^{\left(k\right)}$ and solve the system formed by the union of equations (3a),  (3b) and (4c). 4: Compute $y^{(k+1)}=Q_{\mathrm{prod}}^{(k+1)}$ using equation (5a). 5: Set $Q_{\mathrm{prod}}=y^{(k+1)}$ and solve the system formed by equations (3c), (3d), (3e),  (3f) and (6e). 6: Compute $\delta {\mathbf{X}}^{\left(k\right)}={\mathbf{X}}^{(k+1)}-{\mathbf{X}}^{\left(k\right)}$ and compute the relative increment between two consecutive iterations $$\begin{array}{ccc}& {ϵ}^{\left(k\right)}={\displaystyle \frac{\parallel \delta {\mathbf{X}}^{\left(k\right)}{\parallel }_{\infty }}{\parallel {\mathbf{X}}^{(k+1)}{\parallel }_{\infty }}}.\hfill & \hfill \end{array}$$ (7) 7: if${ϵ}^{\left(k\right)}<\mathtt{tol}$then 8:     Set ${\mathbf{X}}^{*}={\mathbf{X}}^{(k+1)}$. 9:     Compute $Q_{\mathrm{prod}}$ using equation (5a) and set $Q_{\mathrm{prod}}^{*}=Q_{\mathrm{prod}}$. 10:     Compute $Q_{\mathrm{TM}}$ using equation (6e) and set $Q_{\mathrm{TM}}^{*}=Q_{\mathrm{TM}}$. 11:     Set $Q_{\mathrm{UV}}^{*}=Q_{\mathrm{prod}}^{*}-Q_{\mathrm{TM}}^{*}$. 12:     Exit. 13: else 14:     ${\mathbf{X}}^{\left(0\right)}=\theta {\mathbf{X}}^{(k+1)}+(1-\theta ){\mathbf{X}}^{\left(k\right)}$; 15:     $k\leftarrow k+1$; 16:     go to line 2. 17: end if

8. Results

In this section we illustrate the outcome of a series of simulations aiming at:

• to identify the optimal value of the relaxation parameter $\theta$;
• to investigate the impact of a differential expression of the Na⁺ ionic channels along the BLM of the PE and NPE cells on the production of AH by the CE;
• to investigate the impact of ${\alpha }_Q$ and $p_{\mathrm{ev}}$ on the formation of the IOP and on the drainage of AH by the TM and UV outflow pathways.

Our final goal is to leverage the results from items 2. and 3. to better understand the pathological dysfunctions behind severe ocular diseases such as primary open angle glaucoma and/or thyroid eye disease and to help design patient-specific pharmacological therapies.

8.1. Model Calibration in Baseline Conditions

In this section we perform a calibration of the set of model parameters in such a way that under given physiological values of the boundary data, Algorithm 1 returns approximations to $Q_{\mathrm{prod}}$ and $p_{\mathrm{IO}}$ that are as close as possible to the baseline values for healthy control subjects $Q_{\mathrm{prod},\mathrm{b}}=2.75\mu {\mathrm{Lmin}}^{-1}$ and $p_{\mathrm{IO},\mathrm{b}}=15\mathrm{mmHg}$ indicated in the seminal references [1] and [10].

Following [10], we set ${\overline{p}}_{\mathrm{CC}}=25\mathrm{mmHg}$ and ${\overline{p}}_{\mathrm{ec}}=8\mathrm{mmHg}$; following [32], we set ${\overline{p}}_{\mathrm{scs}}=11\mathrm{mmHg}$; we also set $p_{\mathrm{PC}}^{\left(0\right)}=0\mathrm{mmHg}$, $\theta =0.8$, the obstruction coefficient ${\alpha }_Q=0{\mathrm{mmHg}}^{-1}$ and $\mathtt{tol}={10}^{-8}$. The values of the other main model parameters are reported in Appendix A. A screenshot of the output of the computer code at convergence is illustrated in Figure 6. The percentage discrepancy between the calculated values of $Q_{\mathrm{prod}}$ and $p_{\mathrm{IO}}$ and the baseline values $Q_{\mathrm{prod},\mathrm{b}}$ and $p_{\mathrm{IO},\mathrm{b}}$ is 0.0060 % and 0.0026 %, respectively.

Figure 7 shows a pictorial view of the pathway followed by AH water flow in its motion across the CE towards the aqueous side in the PC. We see that AH flows from the stroma along parallel directions, one intercellularly from half-PE to half-NPE cell, the other from the paracellular space neighbour to the half-PE cell to the paracellular space neighbour to the half-NPE cell. We also see that no water flow is predicted to be exchanged between the half-cells of the CE and the neighbour paracellular space.

Figure 8 graphically shows the number of iterations $N_{it}$ required by the fixed-point Algorithm 1 to converge as a function of the relaxation parameter $\theta$. We see that increasing $\theta$ until 0.8 improves the convergence rate; then, furtherly increasing $\theta$ even leads to convergence failure within the maximum allowed number of iterations (equal to 5000).

8.2. The Role of the Expression of the Na⁺ Channels

In this section we use the model to theoretically support the experimental conclusions reached in [33] regarding the expression of ameloride-sensitive Na⁺ channels along the BLM of the NPE cells and their role in controlling AH production by the CE. Data shown in [33] demonstrated that the correct unidirecthetional secretion of AH from the stroma to the PC is accompanied by a net release of sodium from the NPE cells into the PC resulting from the balance between sodium extrusion by the Na+/K+ ATPase and sodium reabsorption through the ameloride-sensitive Na⁺ channels.

To check this important finding, we conducted two numerical simulations: in the first simulation (S1), we ran Algorithm 1 by mathematically setting to zero the value of the diffusional permeability to sodium on the BLM of the sole PE half cell; in the second simulation (S2), we ran Algorithm 1 by mathematically setting to zero the value of the diffusional permeability to sodium also on the BLM of the NPE half cell. No difference was detected between the convergence history of the solution map of S1 and that illustrated in Figure 6. Conversely, in the convergence history of the solution map of S2 represented in Figure 9 we see that the predicted value of Na⁺ molar density in the PE and NPE half cells is equal to zero. This depletion of intracellular sodium leads to a breakdown of water flow and pressure balance in the whole circuit of the anterior segment of the eye, as witnessed by the ’NaN’ IEEE arithmetic representation for a "Not-a-Number" obtained as a result of mathematically undefined operations like 0.0/0.0 (see [34]).

8.3. The Impact of ${\alpha }_Q$ and $p_{\mathrm{ev}}$

In this section we conduct a sensitivity analysis with respect to the obstruction coefficient ${\alpha }_Q$ in (6d) and the pressure in the episcleral vein $p_{\mathrm{ev}}$.

Figure 10 illustrates the distribution of predicted values of IOP in correspondence of 101 values of ${\alpha }_Q$ equally distributed in the interval $[0,0.2]{\mathrm{mmHg}}^{-1}$ and for 6 values of $p_{\mathrm{ev}}$ equally distributed in the interval $[8,20]\mathrm{mmHg}$. We see that IOP increases by almost 40 % with respect to baseline conditions, reaching a value of 20.9 mmHg which is very close to the upper limit of 21 mmHg that separates "normal" from "pathological" conditions.

Figure 11 illustrates the distribution of predicted values of $Q_{\mathrm{prod}}$ (top panel), $Q_{\mathrm{TM}}$ (middle panel) and $Q_{\mathrm{UV}}$ (bottom panel) in correspondence of the same values of ${\alpha }_Q$ and $p_{\mathrm{ev}}$ used in Figure 10. We see that the predicted VFRs of produced AH and drained AH through the TM outflow pathway decrease by almost 60 % and 93 % with respect to baseline conditions, whereas the predicted VFR of drained AH through the UV outflow pathway increases by more than 3 times with respect to baseline conditions.

9. Discussion

The dynamics of AH in the human eye is the focus of extensive experimental and computational research [1,10,13,15,35]. In this study, we developed a modular mathematical model of AH dynamics capable of simulating the processes of production, passive flow and drainage of AH with an equivalent electric circuit of the anterior segment of the eye. The main advantage of the model lies in the integration of multiple spatial scales ranging from ${10}^{-9}$m to ${10}^{-2}$m within a 0D representation based on the successive solution of systems of nonlinear algebraic equations of small size. This multiscale modeling approach provides an improved capability to capture how phenomena occurring at levels difficult to access experimentally affect the macroscopic function of the eye. The simulation results demonstrated that: (C1) Na⁺channels on the BLM of the NPE cells are necessary to obtain a physiological direction and intensity of AH production; and (C2) the nonlinear dependence of TM outflow facility has a substantial impact on the formation of IOP and AH production. Conclusion C1 is an original contribution of this article because it gives a theoretical confirmation of the experimental findings of [33] about the pivotal role of sodium reabsorption to compensate excessive sodium extrusion due to the action of the sodium/potassium pump on the BLM of NPE cells. Conclusion C2 is another original contribution of this article and has a twofold valence. On the one hand, it extends the range of application of the nonlinear model of TM outflow resistance proposed in [26] to the simulation of the whole process of AH dynamics. On the other hand, it allows us to parametrically investigate the pressure condition of the eye and the rate of AH production with respect to the variation of the obstruction coefficient ${\alpha }_Q$ and episcleral vein pressure ${\overline{p}}_{\mathrm{ev}}$. Results support the conjecture that deterioration of the TM tissue poroelastic properties due to aging and molecular-driven dysfunctions of TM cells may be directly related to an increase of IOP to levels compatible with primary open angle glaucoma development (see [36,37,38]). Results also support the conjecture that Thyroid Eye Disease (TED) may cause secondary glaucoma because of compression of the ophthalmic veins which results into an increase of the episcleral venous pressure and, in turn, of IOP (see [4,5,39]). Finally, results support the conjecture raised in [40] that the UV outflow pathway acts as a "relief valve" when the TM outflow pathway fails at appropriately draining AH out of the AC. We can verify this conclusion by comparing the middle and bottom panels of Figure 11 which show that $Q_{\mathrm{TM}}$ progressively decreases with ${\alpha }_Q$ and $p_{\mathrm{ev}}$ while $Q_{\mathrm{UV}}$ progressively increases.

10. Conclusions

In this article, we developed a compartment-based mathematical model of AH dynamics in the human eye, consisting of three modules (M1, M2 and M3) which describe at a macroscopic scale the processes of AH production, passive flow and drainage.

We used the model to comprehend how the findings of each module may help to explain the onset and progression of an ocular disease and to evaluate their impact on patient’s health by contributing to develop specifically targeted therapeutical approaches. In particular, we investigated the role of (1) sodium channels in AH production and (2) episcleral vein pressure in AH drainage. Model predictions suggest that (1) the expression of sodium channels in the NPE cells of the CE is essential for AH to flow from the ciliary processes into the posterior chamber; and (2) an increase of the episcleral vein pressure reduces AH drainage, potentially explaining the development of secondary open-angle glaucoma.

Compared to existing compartment models, our formulation is characterized by a genuine multiscale nature because M1 and M3 incorporate cellular-scale mechanisms which are difficult to access experimentally though significantly impacting the functions of the eye. The multiscale structure is the point of strength of the model and makes it ideal for use in a clinical setting.

The model is also affected by limitations: (L1) it is restricted to stationary simulations; (L2) it is not connected to the external structures of the anterior chamber, such as the cornea; (L3) it is not connected to the posterior structures of the eye, such as the retina. While (L1) can be easily overcome by introducing a further loop over the time variable in Algorithm 1, (L2) and (L3) require elaborations which will be the object of the next steps of our research.