Open Access
Glaucoma  |   August 2022
A Mathematical Model of Aqueous Humor Production and Composition
Author Affiliations & Notes
  • Mariia Dvoriashyna
    Department of Applied Mathematics and Theoretical Physics, University of Cambridge, United Kingdom
    Mathematical Institute, University of Oxford, United Kingdom
  • Alexander J. E. Foss
    Department of Ophthalmology, Nottingham University Hospitals NHS Trust, United Kingdom
  • Eamonn A. Gaffney
    Mathematical Institute, University of Oxford, United Kingdom
  • Rodolfo Repetto
    Department of Civil, Chemical and Environmental Engineering, University of Genoa, Italy
  • Correspondence: Mariia Dvoriashyna, Mathematical Institute, University of Oxford, Radcliffe Observatory, Andrew Wiles Building, Woodstock Rd, Oxford OX2 6GG, UK; dvoriashyna@maths.ox.ac.uk
Investigative Ophthalmology & Visual Science August 2022, Vol.63, 1. doi:https://doi.org/10.1167/iovs.63.9.1
  • Views
  • PDF
  • Share
  • Tools
    • Alerts
      ×
      This feature is available to authenticated users only.
      Sign In or Create an Account ×
    • Get Citation

      Mariia Dvoriashyna, Alexander J. E. Foss, Eamonn A. Gaffney, Rodolfo Repetto; A Mathematical Model of Aqueous Humor Production and Composition. Invest. Ophthalmol. Vis. Sci. 2022;63(9):1. https://doi.org/10.1167/iovs.63.9.1.

      Download citation file:


      © ARVO (1962-2015); The Authors (2016-present)

      ×
  • Supplements
Abstract

Purpose: We develop a mathematical model that predicts aqueous humor (AH) production rate by the ciliary processes and aqueous composition in the posterior chamber (PC), with the aim of estimating how the aqueous production rate depends on the controlling parameters and how it can be manipulated.

Methods: We propose a compartmental mathematical model that considers the stromal region, ciliary epithelium, and PC. All domains contain an aqueous solution with different chemical species. We impose the concentration of all species on the stromal side and exploit the various ion channels present on the cell membrane to compute the water flux produced by osmosis, the solute concentrations in the AH and the transepithelial potential difference.

Results: With a feasible set of parameters, the model predictions of water flux from the stroma to the PC and of the solute concentrations in the AH are in good agreement with measurements. Key parameters which impact the aqueous production rate are identified. A relevant role is predicted to be played by cell membrane permeability to \(\text{K}^+\) and \(\text{Cl}^-\), by the level of transport due to the Na+-H+ exchanger and to the co-transporter of Na+/K+/2Cl; and by carbonic anhydrase.

Conclusions: The mathematical model predicts the formation and composition of AH, based on the structure of the ciliary epithelium. The model provides insight into the physical processes underlying the functioning of drugs that are adopted to regulate the aqueous production. It also suggests ion channels and cell membrane properties that may be targeted to manipulate the aqueous production rate.

Aqueous humor (AH) formation has been an important topic in ocular research for several decades, owing to its fundamental role in regulating IOP and in the management of glaucoma. AH is an aqueous solution containing a mixture of electrolytes, organic solutes, growth factors, and select proteins1 that fills the posterior and anterior chambers (PC and AC) of the eye (Fig. 1). It is well-established that AH is produced by the ciliary epithelium (CE) at a rate of 1 to 3 µL/min2,3; it flows from the PC into the AC through the pupil and exits the eye via the conventional and uveoscleral pathways.4 The balance between the rate of AH production and resistance to its drainage governs the IOP, which typically ranges between 12 and 22 mm Hg, in healthy subjects. An elevated IOP is correlated with the occurrence of glaucoma.5,6 Decreasing the AH production rate is one of the possible strategies to lower the IOP and treat glaucoma. 
Figure 1.
 
Sketches of a cross section of the human eye (A), of the ciliary processes (B) and of the CE (C).
Figure 1.
 
Sketches of a cross section of the human eye (A), of the ciliary processes (B) and of the CE (C).
The CE is a bilayer of cells that consists of the pigmented epithelium (PE), facing the stroma, and the non-pigmented epithelium (NPE), facing the PC (Fig. 1(B)). Selective tight junctions that connect NPE cells separate the stromal side from the PC side. The two cell layers are connected by gap junctions, effectively forming a functional syncytium (see Fig. 1(C)). Therefore, the epithelium behaves like a secretory monolayer epithelium. This view is reinforced when one considers the distribution of ion channels and transporters, and of the \(\text{Na}^+\)-\(\text{K}^+\) pump, as shown in Figure 2, where the NPE basolateral membrane has a channel distribution reminiscent of an epithelial apical domain.710 
Figure 2.
 
(A) Ion channels in the CE. PE basolateral membrane: nkcc: co-transporter of \(\text{Na}^+\)/\(\text{K}^+\)/2\(\text{Cl}^-\); \({\rm {\small AE}}_s\): \(\text{Cl}^-\)-\(\text{HCO}_3^-\) exchanger; nhe: \(\text{Na}^+\)-\(\text{H}^+\) exchanger; \({\rm {\small NBC}}_s\): \(\text{Na}^+\)- \(\text{HCO}_3^-\) cotransporter, presumed to be of stoichiometry 1:2; \(\text{K}^+\)- channels. NPE basolateral membrane: \({\rm {\small PUMP}}\): \(\text{Na}^+\)-\(\text{K}^+\) ATPase; \({\rm {\small AE}}_p\): \(\text{HCO}_3^-\)-\(\text{Cl}^-\) exchanger; \({\rm {\small NBC}}_p\): \(\text{Na}^+\)- \(\text{HCO}_3^-\) cotransporter; \(\text{K}^+\) and \(\text{Cl}^-\) channels. (B) A sketch of the compartmental model with water and ion fluxes (not to scale). The symbols \(s\), \(c\) and \(p\) denote stroma, intracellular space and PC respectively. The symbols \(\tilde{p},\ \tilde{s}\) and \(tj\) denote the PC and stromal cell membranes and the tight junction. In the direction from region \(m\) to region \(k\), \(J_i^{mk}\) denotes flux of solute \(i\) and \(Q^{mk}\) the water flux (see main text).
Figure 2.
 
(A) Ion channels in the CE. PE basolateral membrane: nkcc: co-transporter of \(\text{Na}^+\)/\(\text{K}^+\)/2\(\text{Cl}^-\); \({\rm {\small AE}}_s\): \(\text{Cl}^-\)-\(\text{HCO}_3^-\) exchanger; nhe: \(\text{Na}^+\)-\(\text{H}^+\) exchanger; \({\rm {\small NBC}}_s\): \(\text{Na}^+\)- \(\text{HCO}_3^-\) cotransporter, presumed to be of stoichiometry 1:2; \(\text{K}^+\)- channels. NPE basolateral membrane: \({\rm {\small PUMP}}\): \(\text{Na}^+\)-\(\text{K}^+\) ATPase; \({\rm {\small AE}}_p\): \(\text{HCO}_3^-\)-\(\text{Cl}^-\) exchanger; \({\rm {\small NBC}}_p\): \(\text{Na}^+\)- \(\text{HCO}_3^-\) cotransporter; \(\text{K}^+\) and \(\text{Cl}^-\) channels. (B) A sketch of the compartmental model with water and ion fluxes (not to scale). The symbols \(s\), \(c\) and \(p\) denote stroma, intracellular space and PC respectively. The symbols \(\tilde{p},\ \tilde{s}\) and \(tj\) denote the PC and stromal cell membranes and the tight junction. In the direction from region \(m\) to region \(k\), \(J_i^{mk}\) denotes flux of solute \(i\) and \(Q^{mk}\) the water flux (see main text).
The components of the AH are ultimately derived from the blood and they have to traverse the endothelium of the ciliary circulation to enter the stroma, by a process of ultrafiltration, cross from the stroma to enter the CE, passing through the basolateral surface of the PE, and cross from the CE to enter the PC, crossing the basolateral surface of the NPE. This article addresses these last two steps, to develop a model that can represent these processes. 
Mathematical modeling has proven to be useful in understanding ocular fluid flows.11 A successful model allows one to explain a particular phenomenon, separate the underlying mechanisms, determine their relative importance and make testable predictions. In the context of AH formation, a mathematical model would be useful to understand the key mechanisms governing AH production rate and suggest possible targets for future drug development. 
Water transport across the CE involves several possible mechanisms. The first is a mechanical pressure difference, \(\Delta p\), between the stroma and the PC that induces a water flux into the PC. The second is an oncotic pressure difference, \(\Delta \Pi _p\), with the oncotic pressure being higher on the stromal side. This drives fluid from the PC towards the stroma. The conventional view is that these two contributions approximately balance each other out.12,13 The main mechanism is thought to be osmotic pressure difference \(\Delta \Pi _s\) between the two sides of the CE, generated by active transport of ions, which produces a net flow towards the PC. In mathematical terms, such mechanisms of AH production are expressed by the Starling law:  
\begin{eqnarray} q=K \left( \Delta p - \sigma _p \Delta \Pi _p - \sigma _s \Delta \Pi _s \right). \quad\end{eqnarray}
(1)
In this expression, \(q\) is water flux per unit surface area (with the dimensions of velocity), \(K\) is the hydraulic conductivity of the CE, and \(\sigma _s\) and \(\sigma _p\) are reflection coefficients for ions and proteins, respectively. Oncotic and osmotic pressures can be evaluated using van’t Hoff’s law \(\Delta \Pi = RT \Delta C\), with \(\Delta C\) being difference in solute concentration (osmolarity) across the cell layer, \(R\) the universal gas constant, and \(T\) the absolute temperature. Equation (1) has been used in several modeling works to evaluate AH production.14,15 However, with this approach the distribution of ion channels and transporters present on the cell membrane is completely disregarded and thus it is not possible to obtain information about their specific role in generating the osmolarity jump across the CE. 
In a recent work, Sacco et al.16 proposed a more detailed model that takes into account the distribution of ion channels on the cell membrane (with the exception of the \(\text{Na}^+\)-\(\text{HCO}_3^-\) channel, which was omitted). The authors impose the concentration of all species both on the stromal and PC sides of the CE and use the model to compute the concentration of all species in the cell. This implies that the model cannot predict how AH production rate would be modified by inhibition of some of these channels, since the osmolarity jump across the CE is, again, prescribed. 
In the present article, we propose an alternative approach that allows us to predict not only ion and water transport, but also AH composition in the PC, given the concentrations of all species in the stroma. Because in our model fluxes and concentrations depend on the distribution of channels and transporters on the cell membranes, we can predict how they would react to the inhibition of specific channels. This is of practical and conceptual importance, because the inhibition of specific channels is the aim of certain drugs used to manipulate the AH production rate. 
Mathematical Model
In this work, we propose a compartmental 0-dimensional model of the transport of water and chemical species across the CE, a schematic of which is presented in Figure 2B. We consider three regions: the stroma (denoted with superscript \(s\)), the cellular compartment, which encompasses both NPE and PE cells (\(c\)) and the PC (\(p\)). The compartments are separated by the stromal membrane (denoted with superscript \(\tilde{s}\)) and the PC membrane (\(\tilde{p}\)). We further assume that the tight junctions act as semi-permeable membranes between stroma and the PC (\(tj\)). The spatial variation within the compartments is neglected and all variables are averaged over their respective domains. This assumption is not valid in the cleft gap between two adjacent cells, where the spatial variation of ion concentration is a key ingredient for generating standing gradient osmotic flow.17 However, in this work we only consider osmotic flow driven by the difference in ion concentration between the PC and the stroma and neglect possible local osmosis in the cleft gap, as further commented upon in the discussion. 
We consider that all domains contain an aqueous solution of seven different species: \(\text{Na}^+\), \(\text{K}^+\), \(\text{Cl}^-\), \(\text{HCO}_3^-\) (bicarbonate), \(\text{H}^+\), \(\text{CO}_2\), and \({\text{H}_2\text{CO}_3}\) (carbonic acid). We denote the concentrations of these solutes as \(C_i^m\) (mM), where subscript \(i\) denotes the species (\(i=1,\dots ,7\), enumerated in the order above) and superscript \(m\) the region \((m=s,c,p)\). We denote their valence and charge with \(z_i\). The cell contains fixed non-diffusible charged solutes, which we denote with \(X\), and their concentration and valence are denoted as \(C^c_X\) and \(z_X\). The electric potential in region \(m=s,c,p\) is denoted with \(V^m\)
Solutes are transported across the membranes of the domain via ion channels, co-transporters, and exchangers and we indicate a flux of a solute \(i\) from region \(m\) to region \(k\) with \(J_i^{mk}\) (mol/s). For instance, the flux of \(\text{Na}^+\) from the stroma into the cell is named \(J_1^{sc}\). We consider that the stromal side of the membrane contains co-transporters of \(\text{Na}^+\), \(\text{K}^+\)and 2\(\text{Cl}^-\)(denoted with nkcc), anion exchangers of \(\text{Cl}^-\) and \(\text{HCO}_3^-\) (\({\rm {\small AE}}_s\)),7 \(\text{Na}^+\)-\(\text{H}^+\) exchangers (nhe),8 co-transporters of \(\text{Na}^+\) and 2\(\text{HCO}_3^-\) (\({\rm {\small NBC}}_s\)) and \(\text{K}^+\) channels,18 as shown in Figure 2A. The PC side of the membrane has \(\text{Na}^+\)-\(\text{K}^+\) ATPase (denoted with \({\rm {\small PUMP}}\)),9 \(\text{Cl}^-\)-\(\text{HCO}_3^-\) exchangers (\({\rm {\small AE}}_p\)), \(\text{Na}^+\)-\(\text{HCO}_3^-\) co-transporters (\({\rm {\small NBC}}_p\)), \(\text{K}^+\) and \(\text{Cl}^-\) channels.18,19 The tight junctions are assumed to be permeable to \(\text{Na}^+\), \(\text{K}^+\), \(\text{Cl}^-\) and \(\text{HCO}_3^-\). In addition, all membranes and the tight junction are permeable to \(\text{CO}_2\) and \({\text{H}_2\text{CO}_3}\)
The expressions for the solute fluxes (mol/s) across each transporter are reported in the Appendix, but they all have a common representation as the flux through channel \(k\) can be formally written as  
\begin{eqnarray} \mathcal {J}_k = A P_{k} f(C,V), \quad\end{eqnarray}
(2)
with \(f\) being a non-linear, dimensionless function that involves solute concentrations \(C\) and potentials \(V\) on both sides of the membrane, \(P_k\) the intensity of channel transport (mol/m\(^2\)/s), and \(A\) (m\(^2\)) the surface area of the considered membrane. 
The flux of solute \(i\) can also be driven across the membrane between regions \(m\) and \(k\) by electrodiffusion through ion channels (or by diffusion for non-charged solutes), and has the generalized representation  
\begin{eqnarray} \mathcal {J}^{mk}_{i} = A P g(C,V), \quad\end{eqnarray}
(3)
where \(P\) is a membrane permeability to this solute (m/s) and \(g(C,V)\) is a function that has dimensions of concentration (mM), with a dependence on solute concentrations \(C\) and potentials \(V\) either side of the membrane (see the Appendix for details). 
The flux of each solute \(i\) from region \(m\) to region \(k\), \(J_i^{mk}\), is the sum of fluxes through all channels and transporters on that membrane, which involve solute \(i\). We also assume that at the outlet of the PC the solute is only driven by advection (\(QC_i^p\), with \(Q\) [m\(^3\)/s] being the rate of aqueous production) and that diffusive flux can be neglected. Indeed, the velocity in thin iris-lens channel is about \(U\approx 0.2\) mm/s,4 so using typical PC length \(H\approx 6\) mm and a diffusion coefficient \(D\approx 2\cdot 10^{-3}\) mm\(^2\)/s,17 we can calculate the Péclet number, \(Pe=UH/D \approx 600\), confirming that advection has a dominant contribution at the PC exit. 
The considered solutes interact through the following chemical reactions  
\begin{eqnarray} \text{HCO}_3^- + \text{H}^+ \mathop{\rightleftharpoons}_{k_{-1}}^{k_1} \text{H}_2\text{CO}_3 \mathop{\rightleftharpoons}_{k_{h}}^{k_d} \text{CO}_2 + \text{H}_2\text{O}, \quad\end{eqnarray}
(4)
with \(k_1\), \(k_{-1}\), \(k_d\), and \(k_h\) being the reaction rate constants. The first reaction is the dissociation of \({\text{H}_2\text{CO}_3}\) into the ions \(\text{HCO}_3^-\) and \(\text{H}^+\), which is rapid and happens almost instantaneously. The second reaction is normally slow, but can be catalyzed up to six orders of magnitude in speed by carbonic anhydrase (CA).20 CA II and IV are present in the NPE cells.10 We denote with \(R_i^m\) the reaction terms (production minus consumption) of species \(i\) in region \(m\), which we model with the mass action law21 (see Appendix). 
Water is transported across each membrane by osmosis and we denote the water flux from region \(m\) to region \(k\) with \(Q^{mk}\) (m\(^3\)/s) and it has the following generic expression  
\begin{eqnarray} Q^{mk} = - A K \sigma RT \Delta C^{mk}, \quad\end{eqnarray}
(5)
where \(\Delta C^{mk}\) is the difference in osmolarity between the two sides of the membrane, \(K\) (m/s/Pa) is hydraulic conductivity, and \(A\) is the surface area of the membrane. For example, water flux from stroma to the cell can be written as \(Q^{sc} = -A^{\tilde{s}} K^{\tilde{s}} \sigma RT \left[\sum _{i=1}^7 \left(C^s_i - C^c_i\right) - C_X \right]\)
In the stromal compartment all species concentrations \(C_i^s\) and the electric potential are imposed. 
We now may finally write the steady state conservation of mass in the cells and the PC as follows  
\begin{eqnarray} \begin{array}{@{}c@{\qquad}c@{}} \mbox{Cell} & \mbox{PC} \\ \\ J_i^{sc}- J_i^{cp} + R_i^c=0, & J_i^{cp}+J_i^{sp}-Q C_i^p + R_i^p=0, \end{array} \quad\end{eqnarray}
(6a)
 
\begin{eqnarray} Q^{sc}-Q^{cp}=0, \quad Q^{cp}+Q^{sp}-Q=0, \quad\end{eqnarray}
(6b)
with \(i=1,\dots ,7\). These equations are simplified by taking into account the fact that the first step of reaction (4) is almost instantaneous (please refer to the Appendix for details). The equations are complemented by electroneutrality in both the cell and the PC  
\begin{eqnarray} \sum _{i=1}^7 z_i C_i^c+z_X C_X=0,\quad \sum _{i=1}^7 z_i C_i^p=0. \quad\end{eqnarray}
(6c)
 
In total, we have 18 nonlinear algebraic equations for an equal number of unknowns: concentrations in the cell and the PC: \(C_i^c\), \(C_i^p\) (14 in total), \(X\), water flux \(Q\), and potentials in the cell and the PC \(V^c\), \(V^p\). These equations are solved in Matlab using the function solve for symbolic variables and the stability of the solution was verified by solving a time-dependent version of the model, using the Matlab function ODE45, which is based on an explicit Runge-Kutta (4,5) formula.22 
Model parameters are reported in Table 1 and discussed in detail along with other model parameters in the Appendix
Table 1.
 
Model Parameters.
Table 1.
 
Model Parameters.
Global Sensitivity Analysis
We used a global sensitivity analysis owing to the large number of parameters involved. This strategy allowed us to identify those that affect the model results the most. In particular, we used a variance-based method (extended Fourier amplitude sensitivity test [eFAST]), proposed in Saltelli et al.23 and adapted to biological systems by Marino et al.24 Inter alia, eFAST produces a dimensionless total sensitivity index for each parameter, which estimates the variance associated with the variation of the parameter, including its interactions with the other parameters. In simple terms, the greater the value of the index for a given parameter the higher its influence on the model result under consideration. All parameters considered in the sensitivity analysis were checked to be structurally identifiable using the scaling method proposed by Castro and de Boer.25 
Results
Baseline Values
We first show the results obtained running the model with the ‘baseline’ parameters values, estimated as detailed in the Appendix and summarized in Table 1
In the first part of Table 2 (Experimental measurements), we report available experimental measurements from the literature of the concentrations of various solutes and of the electric potentials. We show the corresponding values predicted by the model in the second section of the table (Model - baseline). Note that model variables in the stromal compartment are prescribed, and have been slightly adjusted with respect to the experimental ones, as the latter do not satisfy the electroneutrality condition, which is one of the assumptions underlying the model and a fundamental physical constraint away from Debye layers. All concentration values predicted by the model in the cell and PC are in good quantitative agreement with experiments and so is the transepithelial potential difference (\(V^p\)). 
Table 2.
 
Concentrations of Ions, pH, and the Potential, in the Cell, Stroma, and the PC.
Table 2.
 
Concentrations of Ions, pH, and the Potential, in the Cell, Stroma, and the PC.
In Table 3, we report ion and water fluxes. With a baseline set of parameters, the model is capable of predicting the AH formation rate (i.e., the flux from the stroma into the PC) with the correct order of magnitude. The model predicts that all the ion fluxes are directed towards the PC; \(\text{Cl}^-\) flux has a magnitude compatible with measurements, whereas the \(\text{Na}^+\) flux is overestimated, which is considered further in the Discussion. 
Table 3.
 
Fluxes of Ions and Water From the Literature and Predicted by the Model.
Table 3.
 
Fluxes of Ions and Water From the Literature and Predicted by the Model.
Sensitivity Analysis
In the previous section, we have shown that with the ‘baseline’ set of parameter values the model satisfactorily reproduces experimental observations, though noting the above-mentioned discrepancy with the \(\text{Na}^+\) flux. In this section, we use a global sensitivity analysis, to identify the intensities of ion fluxes and permeabilities of ion channels that affect the model results the most. 
Ion and water fluxes originate from the active pumping of the \(\text{Na}^+\)-\(\text{K}^+\) ATPase, thus the pump intensity obviously influences all the model results. Indeed, reduction of the parameter \(P_{{\rm {\small PUMP}}}\) (see Equation (2)) by 50% resulted in average decrease of AH secretion by about 25% (results not shown). 
In the sensitivity analysis reported in Figure 3, we fix the intensity of the \(\text{Na}^+\)-\(\text{K}^+\) ATPase pump and vary the intensities of the other channels and transporters in Table 1 (last section, [P-ions]), for a total number of ten parameters. We vary each parameter within \(\pm\)50% of the value in Table 1 and we perform 5973 model runs for each varied parameter, resulting in 59,730 simulations in total. The way the model parameter space is spanned is chosen by the eFAST method itself, given sampling parameters. In particular, the values of all model parameters are varied within their respective range harmonically, each one at a specific frequency.23 We study their effect on AH production rate (Fig. 3A and on the concentrations of \(\text{Na}^+\), \(\text{K}^+\), \(\text{Cl}^-\) and \(\text{HCO}_3^-\) in the PC (Fig. 3B). In the figure the bars height indicates the value of the total sensitivity index; the arrows at the top of each bar show how the given parameter influences the model output: an upwards pointing arrow indicates that the model output under consideration increases as the parameter is increased. We note that we have also varied parameters in a range of \(\pm 25\%\) and \(\pm 75\%\) of the baseline values of Table 1, with no substantive qualitative changes in the results. 
Figure 3.
 
Total sensitivity index for (A) production rate \(Q\) and (B) concentrations in the PC. The arrows indicate how the given parameter influences the model output and are based on the total sensitivity index, which is an output of the eFAST method.23 The vertical dashed line separates the intensities of ion channels, \(P_k\) in Equation (2), from the permeabilities \(P\) in Equation (3).
Figure 3.
 
Total sensitivity index for (A) production rate \(Q\) and (B) concentrations in the PC. The arrows indicate how the given parameter influences the model output and are based on the total sensitivity index, which is an output of the eFAST method.23 The vertical dashed line separates the intensities of ion channels, \(P_k\) in Equation (2), from the permeabilities \(P\) in Equation (3).
Figure 3A shows that AH production rate is very sensitive to the permeability to \(\text{K}^+\)(\(P_{K^+}^{\tilde{p}}\)), and \(\text{Cl}^-\)(\(P_{Cl^-}^{\tilde{p}}\)) and to the intensities of the \(\text{Na}^+\)-\(\text{H}^+\)exchanger (\(P_{\rm {\small NHE}}\)) and the co-transporter of \(\text{Na}^+\)/\(\text{K}^+\)/2\(\text{Cl}^-\) (\(P_{\rm {\small NKCC}}\)). The intensity of the nhe also has significant influence on ion concentrations in the AH (Fig. 3B). 
Role of CA Inhibition
In this section, we investigate how CA inhibition modifies the model predictions. We recall that CA catalyzes the second reaction in Equation (4), which in the absence of CA is slow. In our model, we impose CA inhibition by reducing reaction rates \(k_d\) and \(k_h\) in reaction (4) by a factor of \(10^6\).20 
Results for baseline parameters except that CA is inhibited are reported in Tables 2 and 3 Model baseline - no CA. The \(\text{HCO}_3^-\) concentration in the PC drops significantly and \(\text{K}^+\)concentration increases. The water flux is decreased by approximately 45%. We then perform a parameter variation similar to the one described for sensitivity analysis, but with inhibited CA. The results are reported in Figures 4 and 5, where blue bars refer to the normal case and orange bars to CA inhibition. In the figures, bar height is computed by averaging all values obtained from parameter variation and error bars represent the corresponding standard deviation. In Figure 4A, we show that the model reproduces the reduction of aqueous production as a response to CA inhibition by 40% on average, which is in agreement with experimental observations.31 The model also predicts an increase of \(\text{K}^+\) and \(\text{Cl}^-\)concentration in the PC and a corresponding decrease in \(\text{Na}^+\) and \(\text{HCO}_3^-\)
Figure 4.
 
Change of (A) production rate \(Q\) and (B) concentrations in the PC with CA inhibition.
Figure 4.
 
Change of (A) production rate \(Q\) and (B) concentrations in the PC with CA inhibition.
Figure 5.
 
Change of ion fluxes through each individual channel with CA inhibition (see Appendix for detailed expression of ion channel fluxes). The value is positive if the flux is directed outside the cell. For anion exchangers, the value is positive if \(\text{HCO}_3^-\) goes out of the cell and flux through nhe is positive if \(\text{H}^+\) is transported out of the cell.
Figure 5.
 
Change of ion fluxes through each individual channel with CA inhibition (see Appendix for detailed expression of ion channel fluxes). The value is positive if the flux is directed outside the cell. For anion exchangers, the value is positive if \(\text{HCO}_3^-\) goes out of the cell and flux through nhe is positive if \(\text{H}^+\) is transported out of the cell.
Finally, Figure 5 reports the variation of all ion fluxes through each individual channel. We note that nhe is decreased by about six times with inhibition of CA. Furthermore, nbc on both membranes reverse their transport direction, owing to a change in concentration of \(\text{HCO}_3^-\) in the cell and the PC. 
Discussion and Conclusions
We have presented a mathematical model that successfully captures many of the observed features of aqueous production. The model explains AH production on the basis of osmosis induced by selective ion pumping across the CE. This is consistent with the standard view that the contributions of oncotic and hydrostatic pressure differences, are both small and of opposite sign, and so cancel each other out.12,13 
We account for the presence of seven different solutes (\(\text{Na}^+\), \(\text{K}^+\), \(\text{Cl}^-\), \(\text{HCO}_3^-\), \(\text{H}^+\), \(\text{CO}_2\), and \({\text{H}_2\text{CO}_3}\)) that can chemically interact. We impose the concentrations of all species in the stroma (essentially a fixed boundary condition) and compute water and ion fluxes across the CE together with AH composition. This process is a significant improvement on the previous mathematical models of AH production, where the species concentrations both at the stroma and PC were prescribed. The advantage of our approach is that the model can now predict how inhibition of a specific channel would modify AH production rate and composition. 
We identify a ‘baseline’ set of parameter values, with most of the parameters being estimated on the basis of experimental measurements. With this set of parameters, the model predicts the generation of a AH flux from the stroma into the PC. This flow has the correct order of magnitude, although it slightly underestimates the experimental measurements (Tables 2 and 3). This outcome might be a consequence of the fact that we neglected the role of the cleft gaps among adjacent cells in producing a standing gradient osmotic flux, as described by Diamond and Bossert.32 
With the baseline set of parameters, the model predicts concentrations in the cell and the PC very similar to those measured experimentally (Table 3). Ion fluxes are predicted to be directed towards the PC, which confirms the generally accepted paradigm that “water follows the ions.” The \(\text{Cl}^-\) flux has a magnitude in line with measurements, whereas the \(\text{Na}^+\) flux is overestimated (Table 2). We note, however, that in the measurement of \(\text{Na}^+\) current,28 the authors detected only about 25% of \(\text{Na}^+\)channels,18 which suggests that the measured flux may be about four times higher than that reported in Table 2, resulting in the value of \(3.72\cdot 10^{-6}\) mol/m\(^2\)/s, which is the same order of magnitude as our prediction. 
Experimental evidence shows that the inhibition of the \(\text{Na}^+\)-\(\text{K}^+\) ATPase markedly increases the intracellular \(\text{Na}^+\) content, and reduces the intracellular \(\text{K}^+\) content,26 although the authors do not report quantitative results. This behavior is also reproduced by our model (not shown in the figures). 
We used a global sensitivity analysis to study the influence of each model parameter on the predicted AH production rate and AH chemical composition. We found that the AH production rate is very sensitive to the cell membrane permeabilities to \(\text{K}^+\) and \(\text{Cl}^-\) and to the intensities of the nhe exchanger and of the nkcc co-transporter. The properties of nhe also have a significant influence on ion concentrations in the AH. 
Because the IOP is the result of a balance between AH production rate and resistance to its drainage, one of the possible strategies to decrease the IOP is to use drugs that act on AH production rate. Avila et al.33 showed that topical application of different direct blockers of the nhe exchanger produce IOP reduction in mice, arguably secondary to a reduction in AH production. CA inhibitors are among the most effective drugs used to decrease the rate of AH production; for instance, CA inhibits water flux and results in reduction of water flux by 21%-41%.31 We studied with the model the effect of inhibiting CA and found that this leads to a significant decrease in AH production rate, in agreement with clinical observations. Thus, the model provides a sound physical basis to explain the role of CA on AH production. Furthermore, Counillon et al.8 explained the effect of CA inhibition on AH formation through the reduction of \(\text{H}^+\) and \(\text{HCO}_3^-\) concentrations, which leads to reduction of fluxes through NHE-1 and AE2 channels. This finding is confirmed by our model predictions and supports this explanation for the effect of CA. 
The mechanisms underlying AH production are extremely complicated and the model proposed only focuses on some aspects and we deliberately neglect effects that likely have some influence, but would complicate the picture significantly. Among the most important limitations of our model is that we neglect interactions with blood flow. Kiel et al.13 proposed a lumped parameter mathematical model that couples blood flow in the choroid and ciliary processes (including autoregulation mechanisms), oxygen delivery to and consumption by ocular tissues, and AH production. They obtain results in agreement with experimental observations that highlight the interplay between ciliary blood flow and AH production. 
Another limitation of our model is that it does not account for local osmosis, that is, water transport generated as a consequence of concentration gradients within the cleft gaps among cells.32 This mechanism was found to be by far the dominant mechanism inducing water transport across the retinal pigment epithelium,17,34 and it is likely to influence AH production. Nonetheless, although local osmosis is fundamental for water transport across the retinal pigment epithelium because there is no osmolarity jump across this cell layer, such an osmolarity jump exists across the CE and is correctly predicted by our model. Thus, the role of local osmosis is relatively diminished, justifying the use of a zero-dimensional lumped parameter model without the mechanism of local osmosis. 
Acknowledgments
Disclosure: M. Dvoriashyna, None; A.J.E. Foss, None; E.A. Gaffney, None; R. Repetto, None 
References
Fautsch MP, Johnson DH. Aqueous humor outflow: what do we know? Where will it lead us?. Investig Ophthalmol Vis Sci. 2006; 47(10): 4181–4187. [CrossRef]
Brubaker RF . Measurement of aqueous flow by fluorophotometry. In: Ritch R, Shields MB, Krupin T, eds. The Glaucomas, Volume 2. St. Louis: Mosby; 1989.
Brubaker RF . Flow of aqueous humor in humans (The Friedenwald Lecture). Investig Ophthalmol Vis Sci. 1991; 32(13): 3145–3166.
Dvoriashyna M, Repetto R, Romano MR,, Tweedy JH. Aqueous humour flow in the posterior chamber of the eye and its modifications due to pupillary block and iridotomy. Math Med Biol. 2017; 35: 447–467. [CrossRef]
Kwon YH, Fingert JH, Kuehn MH, Alward WL. Primary open-angle glaucoma. N Engl J Med. 2009; 360(11): 1113–1124. [CrossRef] [PubMed]
Rudnicka AR, Mt-Isa S, Owen CG, Cook DG, Ashby D. Variations in primary open-angle glaucoma prevalence by age, gender, and race: a Bayesian meta-analysis. Invest Ophthalmol Vis Sci. 2006; 47(10): 4254–4261. [CrossRef] [PubMed]
Hamann S . Molecular mechanisms of water transport in the eye. Int Rev Cytol. 2002; 215: 395–431.
Counillon L, Touret N, Bidet M, et al Na\(^+\)/H\(^+\) and \(^-\)/HCO\(_3^-\)–antiporters of bovine pigmented ciliary epithelial cells. Pflügers Archiv. 2000; 440(5): 667–678. [CrossRef] [PubMed]
Riley MV, Kishida K. ATPases of ciliary epithelium: cellular and subcellular distribution and probable role in secretion of aqueous humor. Exp Eye Res. 1986; 42(6): 559–568. [CrossRef] [PubMed]
Shahidullah M, To C-H, Pelis RM, Delamere NA. Studies on bicarbonate transporters and carbonic anhydrase in porcine nonpigmented ciliary epithelium. Invest Ophthalmol Vis Sci. 2009; 50(4): 1791–1800. [CrossRef] [PubMed]
Siggers JH, Ethier CR. Fluid mechanics of the eye. Annu Rev Fluid Mech. 2012; 44: 347–372. [CrossRef]
Delamere NA . Ciliary body and ciliary epithelium. Adv. Organ Biol. 2005; 10: 127–148. [CrossRef] [PubMed]
Kiel J, Hollingsworth M, Rao R, Chen M, Reitsamer H. Ciliary blood flow and aqueous humor production. Prog Retin Eye Res. 2011; 30(1): 1–17. [CrossRef] [PubMed]
Lyubimov G, Moiseeva I, Stein A. Dynamics of the intraocular fluid: Mathematical model and its main consequences. Fluid Dyn. 2007; 42(5): 684–694. [CrossRef]
Szopos M, Cassani S, Guidoboni G, et al Mathematical modeling of aqueous humor flow and intraocular pressure under uncertainty: towards individualized glaucoma management. J Model Ophthalmol. 2016; 1(2): 29–39.
Sacco R, Guidoboni G, Jerome JW, et al A theoretical approach for the electrochemical characterization of ciliary epithelium. Life. 2020; 10(2): 8. [CrossRef]
Dvoriashyna M, Foss AJ, Gaffney EA, Jensen OE, Repetto R. Osmotic and electroosmotic fluid transport across the retinal pigment epithelium: a mathematical model. J Theor Biol. 2018; 456: 233–248. [CrossRef] [PubMed]
Jacob T, Civan M. Role of ion channels in aqueous humor formation. Am J Physiol Cell Physiol. 1996; 271(3): C703–C720. [CrossRef]
Coca-Prados M, Anguita J, Chalfant M, Civan M. Pkc-sensitive cl-channels associated with ciliary epithelial homologue of picln. Am J Physiol Cell Physiol. 1995; 268(3): C572–C579. [CrossRef]
Krahn T, Weinstein AM. Acid/base transport in a model of the proximal tubule brush border: impact of carbonic anhydrase. Am J Physiol Ren Physiol. 1996; 270(2): F344–F355. [CrossRef]
Keener JP, Sneyd J. Mathematical Physiology. New York:Springer; 2009: x1.
Shampine LF, Reichelt MW. The matlab ode suite. SIAM Journal on Scientific Computing. 1997; 18(1): 1–22. [CrossRef]
Saltelli A, Tarantola S, Chan K-S. A quantitative model-independent method for global sensitivity analysis of model output. Technometrics. 1999; 41(1): 39–56. [CrossRef]
Marino S, Hogue IB, Ray CJ, Kirschner DE. A methodology for performing global uncertainty and sensitivity analysis in systems biology. J Theor Biol. 2008; 254(1): 178–196. [CrossRef] [PubMed]
Castro M, de Boer RJ. Testing structural identifiability by a simple scaling method. PLOS Comput Biol. 2020; 16(11): e1008248. [CrossRef] [PubMed]
Bowler JM, Peart D, Purves RD, Carr´e DA, Macknight AD, Civan MM. Electron probe x-ray microanalysis of rabbit ciliary epithelium. Exp Eye Res. 1996; 62(2): 131–140. [CrossRef] [PubMed]
To C, Kong C, Chan C, Shahidullah M, Do C. The mechanism of aqueous humour formation. Clin Exp Optom. 2002; 85(6): 335–349. [PubMed]
Fain GL, Farahbakhsh NA. Voltage-activated currents recorded from rabbit pigmented ciliary body epithelial cells in culture. J Physiol. 1989; 418(1): 83–103. [CrossRef] [PubMed]
Holland MG, Gipson CC. Chloride ion transport in the isolated ciliary body. Invest Ophthalmol Vis Sci. 1970; 9(1): 20–29.
Candia OA, To C-H, Law CS. Fluid transport across the isolated porcine ciliary epithelium. Investig Ophthalmol Vis Sci. 2007; 48(1): 321–327. [CrossRef]
To C-H, Do C-W, Zamudio AC, Candia OA. Model of ionic transport for bovine ciliary epithelium: effects of acetazolamide and HCO\(_3^-\). Am J Physiol Cell Physiol. 2001; 280(6): C1521–C1530. [CrossRef] [PubMed]
Diamond JM, Bossert WH. Standing-gradient osmotic flow: a mechanism for coupling of water and solute transport in epithelia. J Gen Physiol. 1967; 50(8): 2061–2083. [CrossRef] [PubMed]
Avila MY, Seidler RW, Stone RA, Civan MM. Inhibitors of NHE-1 Na\(^+\)/H\(^+\) exchange reduce mouse intraocular pressure. Invest Ophthalmol Vis Sci. 2002; 43(6): 1897–1902. [PubMed]
Dvoriashyna M, Foss AJ, Gaffney EA, Repetto R. Fluid and solute transport across the retinal pigment epithelium: a theoretical model. J R Soc Interface. 2020; 17(163): 20190735. [CrossRef] [PubMed]
Keener JP. Sneyd J, Mathematical Physiology, New York:Springer; 2009.
Beard DA, Qian H. Relationship between thermodynamic driving force and one-way fluxes in reversible processes. PloS One. 2007; 2(1): e144. [CrossRef] [PubMed]
Strieter J, Stephenson JL, Palmer LG, Weinstein AM. Volume-activated chloride permeability can mediate cell volume regulation in a mathematical model of a tight epithelium. J Gen Physiol. 1990; 96(2): 319–344. [CrossRef] [PubMed]
Chu T-C, Candia OA. Electrically silent Na+ and Cl fluxes across the rabbit ciliary epithelium. Invest Ophthalmol Vis Sci. 1987; 28(3): 445–450. [PubMed]
Hamann S, Herrera-Perez JJ, Zeuthen T, Alvarez-Leefmans FJ. Cotransport of water by the Na+-K+-2Cl cotransporter NKCC1 in mammalian epithelial cells. J Physiol. 2010; 588(21): 4089–4101. [CrossRef] [PubMed]
Weinstein AM . A mathematical model of rat distal convoluted tubule. II. Potassium secretion along the connecting segment. Am J Physiol Renal Physiol. 2005; 289(4): F721–F741. [CrossRef] [PubMed]
Itel F, Al-Samir S, Öberg F, et al. CO2 permeability of cell membranes is regulated by membrane cholesterol and protein gas channels. FASEB. 2012; 26(12): 5182–5191. [CrossRef]
Weiss T . Cellular Biophysics I. Transport. Cambridge, MA: MIT Press; 1996.
Weiss T . Cellular Biophysics II. Electrical Properties. Cambridge, MA: MIT Press; 1996.
Appendix 1
Mathematical Description of Ion Channels
In this section, we report the mathematical expressions used to model ion fluxes through the \(\text{Na}^+\)-\(\text{K}^+\) ATPase pump, ion channels and cotransporters as generalized in Equations (2) and (3)
Non charged Solute
We start with the non-charged solutes: \(\text{CO}_2\) and \({\text{H}_2\text{CO}_3}\) (\(i=6,7\)). Fluxes of these species across the membrane are modeled using an expression derived from the integration of the diffusion equation across the membrane,35 
\begin{eqnarray} \mathcal {J}_i^{mk} = A P (C^m_i - C^k_i), \quad\end{eqnarray}
(7)
where \(A\) (m\(^2\)) is the surface area of the membrane and \(P\) (m/s) the permeability of the membrane. This expression is used to model the fluxes of \(\text{CO}_2\) and \({\text{H}_2\text{CO}_3}\) across each membrane and the tight junction. 
Ion Channels
Ion channels, such as \(\text{K}^+\) channels, are selective and transport single ions down their electrochemical gradient. A frequently used expression for the membrane flux due to these channels is the Goldman-Hodgkin-Katz,35 according to which the flux of ion \(i\) from region \(m\) to region \(k\) can be written as  
\begin{eqnarray} {\cal J}^{mk}_{i} = A P z_i \phi ^{mk} \frac{C^m_i - C^k_i \exp \left(-z_i \phi ^{mk}\right)}{1-\exp \left(-z_i \phi ^{mk}\right)}, \quad\end{eqnarray}
(8)
where factor \(P\) is a membrane permeability to this ion (m/s), \(A\) is the area of the membrane, \(\phi ^{mk} = F(V^m-V^k)/RT\) is the dimensionless potential difference across the membrane, \(F\) is Faraday’s constant, \(R\) the universal gas constant, and \(T\) is temperature. This expression is used to model the flux through \(\text{K}^+\) channels on the cell membrane, \(\text{Cl}^-\) channels on the PC side of the membrane and ion fluxes across the tight junction. 
The expressions (7) and (8) are generalized in Equation (3)
Cotransporters and Antiporters
Cotransporters, such as \(\text{Na}^+\)-\(\text{K}^+\)-\(\text{Cl}^-\) (nkcc), and antiporters are passive channels and perform joint transport of solutes. We model these fluxes as a function of the difference in overall chemical potential across the membrane, an approach that is valid in a near-equilibrium system.36 
Following this approach the flux out of the cell of the nkcc channel can be written as  
\begin{eqnarray} {\cal J}_{\rm {\small NKCC}}= A^{\tilde{s}} P_{\rm {\small NKCC}} \ln \left(\frac{C_{\rm K}^c C_{\rm Na}^c (C_{\rm Cl}^c)^2}{C_{\rm K}^s C_{\rm Na}^s (C_{\rm Cl}^s)^2} \right), \quad\end{eqnarray}
(9)
where \(P_{\rm {\small NKCC}}\) is an intensity parameter. 
A similar expression can be derived for electrogenic fluxes. The flux across the channel that co-transports one \(\text{Na}^+\) and two \(\text{HCO}_3^-\) (nbc) on the stromal membrane can be written as  
\begin{eqnarray} {\cal J}_{{\rm {\small NBC}}_s} = A^{\tilde{s}} P_{{\rm {\small NBC}}_s} \left(\ln \left[ \frac{C_{{\rm Na}}^c\left(C_{\rm HCO_3}^c\right)^2}{C_{\rm Na}^s \left(C_{\rm HCO_3}^s\right)^2}\right] - \phi ^{cs} \right), \quad\end{eqnarray}
(10)
with a similar expression on the PC membrane,  
\begin{eqnarray} {\cal J}_{{\rm {\small NBC}}_p} = A^{\tilde{p}} P_{{\rm {\small NBC}}_p} \left(\ln \left[ \frac{C_{{\rm Na}}^c\left(C_{\rm HCO_3}^c\right)^2}{C_{\rm Na}^p \left(C_{\rm HCO_3}^p\right)^2}\right] - \phi ^{cp} \right). \quad\end{eqnarray}
(11)
 
Antiporters or exchangers, transport one solute species into the cell and another one out of it. For \({\rm {\small AE}}_s\) at the stromal side we may write  
\begin{eqnarray} {\cal J}_{{\rm {\small AE}}_s} = A^{\tilde{s}} P_{{\rm {\small AE}}_s} \ln \left(\frac{C_{\rm Cl}^s C_{\rm HCO_3}^c}{C_{\rm Cl}^c C_{\rm HCO_3}^s}\right), \quad\end{eqnarray}
(12)
and on the PC side  
\begin{eqnarray} {\cal J}_{{\rm {\small AE}}_p} = A^{\tilde{p}} P_{{\rm {\small AE}}_p} \ln \left(\frac{C_{\rm Cl}^p C_{\rm HCO_3}^c}{C_{\rm Cl}^c C_{\rm HCO_3}^p}\right). \quad\end{eqnarray}
(13)
Similarly, for nhe, we have the following expression  
\begin{eqnarray} {\cal J}_{\rm {\small NHE}} = A^{\tilde{s}} P_{\rm {\small NHE}} \ln \left(\frac{C_{\rm H}^s C_{\rm Na}^c}{C_{\rm H}^c C_{\rm Na}^s}\right). \quad\end{eqnarray}
(14)
 
\(\text{Na}^+\)-\(\text{K}^+\) ATPase
The \(\text{Na}^+\)-\(\text{K}^+\) ATPase pump transports three \(\text{Na}^+\) out of the cell and two \(\text{K}^+\) in, both against their concentration gradient. Because the energy levels required to fuel the pump significantly exceed the free energy of the system, near-equilibrium thermodynamics cannot be invoked to model the \(\text{Na}^+\)-\(\text{K}^+\) ATPase, in contrast with cotransporters and antiporters. Therefore, we follow the empirical approach proposed by Strieter et al.,37 according to which  
\begin{eqnarray} {\cal J}_{\rm {\small PUMP}}= A^{\tilde{p}} P_{\rm {\small PUMP}} \left(\frac{C_{\rm Na}^c}{K_{\rm Na}+C_{\rm Na}^c}\right)^3 \left( \frac{C_{\rm K}^p}{K_{\rm K}+C_{\rm K}^p}\right)^2, \quad\end{eqnarray}
(15)
where \(P_{\rm {\small PUMP}}\) (mol/m\(^2\)/s) is an intensity parameter and \(K_{\rm Na}\) and \(K_{\rm K}\) are apparent dissociation constants, which have the following expressions \(K_{\rm Na}=0.2\left(1 + C^c_{\rm K}/\left(8.33\mbox{ mM}\right)\right)\) mM and \(K_{\rm K}=0.1(1 + C^p_{\rm Na}/\left(18.5 \mbox{ mM}\right))\) mM.37 
Expressions (9)(15) are captured in general in the main text in Equation (2)
For every solute \(i\), the overall flux across the cell membrane is the sum of all fluxes through the channels, transporters and exchangers that \(i\) is involved with. For example, \(\text{Na}^+\) is transported across the stromal side of the membrane by nkcc and nbc\(_s\) (both positive if directed out of the cell) and nhe (positive for \(\text{Na}^+\) if directed into the cell); therefore, the overall \(\text{Na}^+\) flux from cell to stroma can be written as follows  
\begin{eqnarray} J_{1}^{cs} ={\cal J}_{{\rm {\small NKCC}}}+{\cal J}_{{\rm {\small NBC}}_s} -{\cal J}_{{\rm {\small NHE}}}. \quad\end{eqnarray}
(16)
\(\text{Na}^+\) is transported through the PC membrane by pump, which transports three \(\text{Na}^+\) ions out of the cell (therefore, the flux in (15) is multiplied by a factor of 3) and nbc\(_p\) (positive is directed out of the cell), hence the flux of \(\text{Na}^+\) from the cell into the PC reads  
\begin{eqnarray} J_{1}^{cp} = 3{\cal J}_{{\rm {\small PUMP}}}+{\cal J}_{{\rm {\small NBC}}_p}. \quad\end{eqnarray}
(17)
The expressions for other ions and membranes can be written analogously, based on the channels reported in Figure 2 and equations reported in this section. 
Reactions
The terms \(R_i^c\) and \(R_i^p\) in Equations (6a) represent source and sink terms owing to the chemical reactions among species; see Equation (4). They are modeled via the mass action law,35 and their expressions read  
\begin{eqnarray} \begin{array}{@{}l@{}} R_4^m = (k_{-1} C^m_7 - k_1 C^m_4 C^m_5)W^m, \\ \\ R_6^m = (- k_h C_6^m + k_d C_7^m)W^m, \\ \\ R_7^m = - R_6^m - R^m_4, \end{array} \quad\end{eqnarray}
(18)
where \(m \in \lbrace c,p\rbrace\), \(W^m\) denotes the volume of compartment \(m\) and all other \(R_k^m\) equal to zero. For the PC, \(W^p = A^{PC}H\) and for the cell, \(W^c = A^{\tilde{s}} L\)
Detailed Summary of Equations
We now rewrite the Equations (6a) using the expressions for the fluxes as reported. We also introduce some simplification based on the fact that the reaction (4) is very fast. 
For ions \(\text{Na}^+\), \(\text{K}^+\) and \(\text{Cl}^-\) (\(k=1,2\), and 3) the flux balance equations in the cell and the PC are written as:  
\begin{eqnarray} &&3{\cal J}_{{\rm {\small PUMP}}}+{\cal J}_{{\rm {\small NKCC}}}+{\cal J}_{{\rm {\small NBC}}_s}+{\cal J}_{{\rm {\small NBC}}_p} -{\cal J}_{{\rm {\small NHE}}}=0 \nonumber\\ &&\quad\text{$\text{Na}^+$ flux balance in the cell}, \end{eqnarray}
(19a)
 
\begin{eqnarray} &&3{\cal J}_{{\rm {\small PUMP}}}+{\cal J}_{{\rm {\small NBC}}_{p}}+{\cal J}^{sp}_1-QC^p_1=0 \nonumber\\ &&\quad\text{$\text{Na}^+$ flux balance in the PC}, \end{eqnarray}
(19b)
 
\begin{eqnarray} && {\cal J}^{cp}_2-2{\cal J}_{{\rm {\small PUMP}}}+{\cal J}^{cs}_2+{\cal J}_{{\rm {\small NKCC}}}=0 \nonumber\\ &&\quad \text{$\text{K}^+$ flux balance in the cell}, \end{eqnarray}
(19c)
 
\begin{eqnarray} && -2{\cal J}_{{\rm {\small PUMP}}}+{\cal J}^{cp}_2+{\cal J}^{sp}_2-QC^p_2=0 \nonumber\\ &&\quad \text{$\text{K}^+$ flux balance in the PC}, \end{eqnarray}
(19d)
 
\begin{eqnarray} && {\cal J}^{cp}_3+2{\cal J}_{{\rm {\small NKCC}}}-{\cal J}_{{\rm {\small AE}}_p}-{\cal J}_{{\rm {\small AE}}_s}=0 \nonumber\\ &&\quad \text{$\text{Cl}^-$ flux balance in the cell}, \end{eqnarray}
(19e)
 
\begin{eqnarray} && {\cal J}^{cp}_3-{\cal J}_{{\rm {\small AE}}_p}+{\cal J}^{sp}_3-Q C^p_3=0 \nonumber\\ &&\quad \text{$\text{Cl}^-$ flux balance in the PC}. \end{eqnarray}
(19f)
Note that the superscript \(sp\) in ion fluxes stands for the fluxes from stroma into the PC through the tight junctions. We also simplify the rest of the equations by using the fact that the second step of the reaction (4) happens instantaneously.20 This means that the equations for \(k=4,5\) can be replaced by a conservation of charge, which is obtained by subtracting Equation (6a) for \(k=5\) from Equation (6a) for \(k=4\) and reads  
\begin{eqnarray} && {\cal J}_{{\rm {\small AE}}_s}+{\cal J}_{{\rm {\small AE}}_p}+2{\cal J}_{{\rm {\small NBC}}_p}+2{\cal J}_{{\rm {\small NBC}}_s} - {\cal J}_{{\rm {\small NHE}}}=0 \nonumber\\ &&\quad \text{conservation of charge in the cell}, \end{eqnarray}
(19g)
 
\begin{eqnarray} && {\cal J}_{{\rm {\small AE}}_p}+2{\cal J}_{{\rm {\small NBC}}_p}+{\cal J}_4^{sp} - {\cal J}_5^{sp} -Q \left(C_4^p-C_5^p\right)=0 \nonumber\\ &&\quad \text{conservation of charge in the PC}, \end{eqnarray}
(19h)
and quasi-equilibrium of the reaction \(R_4^m=0\),  
\begin{eqnarray} C_7^{c} = K_{d} C_4^{c} C_5^{c} \quad\;\, \text{in the cell}, \end{eqnarray}
(19i)
 
\begin{eqnarray} C_7^{p} = K_{d} C_4^{p} C_5^{p} \quad \text{in the PC}. \end{eqnarray}
(19j)
In this expression, \(K_{d} = k_1/k_{-1}\) is the equilibrium constant, the choice of which is described in the next section. Furthermore, the equation for \({\text{H}_2\text{CO}_3}\) (\(k=7\)) is replaced by the conservation of total \(\text{CO}_2\), obtained by summing up Equations (6a) for \(k=5\), \(k=6\), and \(k=7\) and read  
\begin{eqnarray} \!\!\!\!\!\!\!\!\!\!\begin{array}{@{}l@{}} {\cal J}_7^{cp}+{\cal J}_6^{cp}+{\cal J}_7^{cs}+{\cal J}_6^{cs}+{\cal J}_{{\rm {\small NHE}}}=0 \;\; \text{in the cell}, \end{array} \end{eqnarray}
(19k)
 
\begin{eqnarray} \!\!\!\!\!\!\!\!\!\begin{array}{@{}l@{}} {\cal J}_7^{cp}+{\cal J}_6^{cp}+{\cal J}_7^{sp}+{\cal J}_6^{sp}+{\cal J}_5^{sp}\;-\\[.5em] \quad Q\left(C_5^p+C_6^p+C_7^p\right)=0 \;\; \text{in the PC}. \end{array}\end{eqnarray}
(19l)
 
Finally, the equations for \(\text{CO}_2\) (\(k=6\)) read  
\begin{eqnarray} &&{\cal J}_6^{cp}+{\cal J}_6^{cs} + W^c(k_h C^c_6-k_d C^c_7) =0 \nonumber\\ &&\quad \text{$\text{CO}_2$ flux balance in the cell}, \end{eqnarray}
(19m)
 
\begin{eqnarray} &&{\cal J}_6^{cp}+{\cal J}_6^{sp}-C_6^pQ -W^p(k_h C^p_6-k_d C^p_7) =0 \nonumber\\ &&\quad \text{$\text{CO}_2$ flux balance in the PC}. \end{eqnarray}
(19n)
These equations are solved along with electroneutrality Equations (6c) and water balance Equations (6b), which we rewrite here for the sake of clarity.  
\begin{eqnarray} \!\!\!\!\! \sum _{i=1}^7 z_i C_i^c+z_X C_X=0 \quad \text{electroneutrality in the cell}, \end{eqnarray}
(19o)
 
\begin{eqnarray} \sum _{i=1}^7 z_i C_i^p=0 \quad \text{electroneutrality in the PC}, \end{eqnarray}
(19p)
 
\begin{eqnarray} Q^{sc}-Q^{cp}=0 \quad \text{water flux balance in the cell}, \end{eqnarray}
(19q)
 
\begin{eqnarray} Q^{cp}+Q^{sp}-Q=0 \quad \text{water flux balance in the PC}. \end{eqnarray}
(19r)
 
We thus have 18 nonlinear algebraic equations. The unknowns are the concentrations of solutes in the cell and the PC (\(7+7\)), the potential in the cell and the PC (\(1+1\)), the concentration of fixed non-diffusible charged solutes in the cell (1) and water flux (1). 
Appendix 2
Estimation of Parameter Values
The model relies on a number of parameters, some of which are difficult to estimate, owing to the sparse availability of experimental measurements. In this section, we explain the parameters reported in the Table 1 using the markers from the table sections for easier referencing. 
  • GP We consider that the cell has length of \(L=10\) µm and the PC has length \(H = 6\cdot 10^{-3}\) m. The area of the PC at the boundary with ciliary processes is estimated to be \(A^{PC}=2.6\cdot 10^{-5}\) m\(^2\).4 Owing to the extensive folding the area of the ciliary processes is considered to be \(A^{\tilde{p}}=A^{\tilde{s}}=6\cdot 10^{-4}\) m\(^2\).38 Note that even though cell membrane at the posterior side is highly folded, this folding is accounted for in values of permeabilities to ions and water instead of the area. The area of the tight junctions is assumed to be 1000 times smaller than that of the whole epithelium owing to the small cleft to cell thickness ratio of the order \(\mathcal {O}(10^{-3})\).
  • P-H2O Experimental estimates of hydraulic conductivity of the CE vary from \(1.1\cdot 10^{-11}\) m/s/Pa39 to 6.26–9.5\(\cdot 10^{-11}\) m/s/Pa.16 In this work, we set \(K^{\tilde{s}}=2\cdot 10^{-11}\) m/s/Pa for the stromal membrane. The NPE basal membrane is highly folded; we account for this folding by increasing its permeability to water and ions by a factor 10, so that \(K^{\tilde{p}}=10 K^{\tilde{s}}\).
  • RR The rates of the reactions in (4) for hydration and dehydration of \(\text{CO}_2\) (catalyzed) are taken to be \(k_h = 1.45\cdot 10^3\) 1/s and \(k_d = 4.96 \cdot 10^{5}\) 1/s (note that a fixed water concentration is incorporated in \(k_h\)).40 In the presence of CA inhibition, these rates are reduced by a factor \(10^{-6}\). The equilibrium constant of \({\text{H}_2\text{CO}_3}\) dissociation, \(K_{d} = k_{1}/k_{-1}= 5.3\) 1/mM.34
  • P-C Itel et al.41 suggested that permeability of the cell membrane to \(\text{CO}_2\) can range from about \(10^{-4}\) m/s to \(10^{-2}\) m/s. In our model we impose \(P_{\text{$\text{CO}_2$}}^{\tilde{s}}=1.5\cdot 10^{-3}\) m/s, following the value used in Krahn and Weinstein,20 for the proximal tubule brush border. The permeability to \({\text{H}_2\text{CO}_3}\) is also taken from the same article, to be \(P_{\text{${\text{H}_2\text{CO}_3}$}}^{\tilde{s}}=1.28\cdot 10^{-5}\) m/s.20 The permeabilities of the posterior membrane to both \(\text{CO}_2\) and \({\text{H}_2\text{CO}_3}\) are assumed to be 10 times larger owing to the membrane folding.
  • P-ions The value of ion channel permeabilities are hard to estimate, owing to very sparse existing experimental measurements. The model uses as baseline the values reported in Table 1, which were carefully selected to match measured ion concentrations and fluxes and the transepithelial potential difference. We list the parameters that we took from experiments for easier referencing.
[\(P_{{\rm {\small PUMP}}}\)]: \(\text{Na}^+\) current density was measured to be approximately 9 µA/cm\(^2\),18 and a value of 0.22 µEq/h/cm\(^2\) for the net \(\text{Na}^+\) flux was found in isolated cat ciliary body.29 Thus, \(\text{Na}^+\) flux expressed using the international system of units is in the range of approximately \({[6.11,9.32]\cdot 10^{-7}}\) mol/m\(^2\)/s. These values might be underestimated (see the Discussion), and we have selected the magnitude of the \(\text{Na}^+\)-\(\text{K}^+\) ATPase to be \(6\cdot 10^{-6}\) mol/m\(^2\)/s. 
[\(P_{\text{$\text{Cl}^-$}}^{\tilde{p}}\)]: \(J_{Cl}^{net}\) is measured to be approximately \([3,8.23]\cdot 10^{-6}\) Eq/s/m\(^2\) from stroma to aqueous.30,29 Using the experimental data reported in Table 2, we estimate the range of permeability to \(\text{Cl}^-\) of the posterior membrane to be within the range \(P_{\text{$\text{Cl}^-$}}^{\tilde{p}} \in [2.97,8.16]\cdot 10^{-8}\) m/s. 
[\(P^{tj}\)]: Assuming that tight junctions are the main players in determining the value of transepithelial resistance, \(TER=0.6\) k\(\Omega \cdot\)cm\(^2\),18 we compute \(\text{Na}^+\) flux through the tight junctions, with the values of concentration from Table 2, to obtain \(3.37\cdot 10^{-6}\) mol/m\(^2\)/s (directed toward stroma). We use this value to estimate the permeability of tight junctions to ions, \(P^{tj} = 9.22\cdot 10^{-7}\) m/s. A higher value of 6\(\cdot 10^{-6}\) m/s was used in Table 1, because it provided the prediction of TEP closer to the experimentally measured value. 
  • Other The temperature considered is \(T=\) 310 K and the reflection coefficient for the ions in Equation (5) is taken to be \(\sigma =1\).42,43 The average charge and valence of fixed non-diffusible species \(X\) in the cell is taken to be \(z_X=-1.5\).
Figure 1.
 
Sketches of a cross section of the human eye (A), of the ciliary processes (B) and of the CE (C).
Figure 1.
 
Sketches of a cross section of the human eye (A), of the ciliary processes (B) and of the CE (C).
Figure 2.
 
(A) Ion channels in the CE. PE basolateral membrane: nkcc: co-transporter of \(\text{Na}^+\)/\(\text{K}^+\)/2\(\text{Cl}^-\); \({\rm {\small AE}}_s\): \(\text{Cl}^-\)-\(\text{HCO}_3^-\) exchanger; nhe: \(\text{Na}^+\)-\(\text{H}^+\) exchanger; \({\rm {\small NBC}}_s\): \(\text{Na}^+\)- \(\text{HCO}_3^-\) cotransporter, presumed to be of stoichiometry 1:2; \(\text{K}^+\)- channels. NPE basolateral membrane: \({\rm {\small PUMP}}\): \(\text{Na}^+\)-\(\text{K}^+\) ATPase; \({\rm {\small AE}}_p\): \(\text{HCO}_3^-\)-\(\text{Cl}^-\) exchanger; \({\rm {\small NBC}}_p\): \(\text{Na}^+\)- \(\text{HCO}_3^-\) cotransporter; \(\text{K}^+\) and \(\text{Cl}^-\) channels. (B) A sketch of the compartmental model with water and ion fluxes (not to scale). The symbols \(s\), \(c\) and \(p\) denote stroma, intracellular space and PC respectively. The symbols \(\tilde{p},\ \tilde{s}\) and \(tj\) denote the PC and stromal cell membranes and the tight junction. In the direction from region \(m\) to region \(k\), \(J_i^{mk}\) denotes flux of solute \(i\) and \(Q^{mk}\) the water flux (see main text).
Figure 2.
 
(A) Ion channels in the CE. PE basolateral membrane: nkcc: co-transporter of \(\text{Na}^+\)/\(\text{K}^+\)/2\(\text{Cl}^-\); \({\rm {\small AE}}_s\): \(\text{Cl}^-\)-\(\text{HCO}_3^-\) exchanger; nhe: \(\text{Na}^+\)-\(\text{H}^+\) exchanger; \({\rm {\small NBC}}_s\): \(\text{Na}^+\)- \(\text{HCO}_3^-\) cotransporter, presumed to be of stoichiometry 1:2; \(\text{K}^+\)- channels. NPE basolateral membrane: \({\rm {\small PUMP}}\): \(\text{Na}^+\)-\(\text{K}^+\) ATPase; \({\rm {\small AE}}_p\): \(\text{HCO}_3^-\)-\(\text{Cl}^-\) exchanger; \({\rm {\small NBC}}_p\): \(\text{Na}^+\)- \(\text{HCO}_3^-\) cotransporter; \(\text{K}^+\) and \(\text{Cl}^-\) channels. (B) A sketch of the compartmental model with water and ion fluxes (not to scale). The symbols \(s\), \(c\) and \(p\) denote stroma, intracellular space and PC respectively. The symbols \(\tilde{p},\ \tilde{s}\) and \(tj\) denote the PC and stromal cell membranes and the tight junction. In the direction from region \(m\) to region \(k\), \(J_i^{mk}\) denotes flux of solute \(i\) and \(Q^{mk}\) the water flux (see main text).
Figure 3.
 
Total sensitivity index for (A) production rate \(Q\) and (B) concentrations in the PC. The arrows indicate how the given parameter influences the model output and are based on the total sensitivity index, which is an output of the eFAST method.23 The vertical dashed line separates the intensities of ion channels, \(P_k\) in Equation (2), from the permeabilities \(P\) in Equation (3).
Figure 3.
 
Total sensitivity index for (A) production rate \(Q\) and (B) concentrations in the PC. The arrows indicate how the given parameter influences the model output and are based on the total sensitivity index, which is an output of the eFAST method.23 The vertical dashed line separates the intensities of ion channels, \(P_k\) in Equation (2), from the permeabilities \(P\) in Equation (3).
Figure 4.
 
Change of (A) production rate \(Q\) and (B) concentrations in the PC with CA inhibition.
Figure 4.
 
Change of (A) production rate \(Q\) and (B) concentrations in the PC with CA inhibition.
Figure 5.
 
Change of ion fluxes through each individual channel with CA inhibition (see Appendix for detailed expression of ion channel fluxes). The value is positive if the flux is directed outside the cell. For anion exchangers, the value is positive if \(\text{HCO}_3^-\) goes out of the cell and flux through nhe is positive if \(\text{H}^+\) is transported out of the cell.
Figure 5.
 
Change of ion fluxes through each individual channel with CA inhibition (see Appendix for detailed expression of ion channel fluxes). The value is positive if the flux is directed outside the cell. For anion exchangers, the value is positive if \(\text{HCO}_3^-\) goes out of the cell and flux through nhe is positive if \(\text{H}^+\) is transported out of the cell.
Table 1.
 
Model Parameters.
Table 1.
 
Model Parameters.
Table 2.
 
Concentrations of Ions, pH, and the Potential, in the Cell, Stroma, and the PC.
Table 2.
 
Concentrations of Ions, pH, and the Potential, in the Cell, Stroma, and the PC.
Table 3.
 
Fluxes of Ions and Water From the Literature and Predicted by the Model.
Table 3.
 
Fluxes of Ions and Water From the Literature and Predicted by the Model.
×
×

This PDF is available to Subscribers Only

Sign in or purchase a subscription to access this content. ×

You must be signed into an individual account to use this feature.

×