**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.

^{1}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/min

^{2}

^{,}

^{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.

^{7}

^{–}

^{10}

^{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.

^{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:

^{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.

^{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.

^{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.

^{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}\).

^{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.

^{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 law

^{21}(see Appendix).

^{22}

^{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}

**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.

^{20}

^{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^-\).

^{12}

^{,}

^{13}

^{32}

^{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.

^{26}although the authors do not report quantitative results. This behavior is also reproduced by our model (not shown in the figures).

^{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.

^{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.

^{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.

**M. Dvoriashyna**, None;

**A.J.E. Foss**, None;

**E.A. Gaffney**, None;

**R. Repetto**, None

*Investig Ophthalmol Vis Sci*. 2006; 47(10): 4181–4187. [CrossRef]

*The Glaucomas*, Volume 2. St. Louis: Mosby; 1989.

*Investig Ophthalmol Vis Sci*. 1991; 32(13): 3145–3166.

*Math Med Biol*. 2017; 35: 447–467. [CrossRef]

*N Engl J Med*. 2009; 360(11): 1113–1124. [CrossRef] [PubMed]

*Invest Ophthalmol Vis Sci*. 2006; 47(10): 4254–4261. [CrossRef] [PubMed]

*Int Rev Cytol.*2002; 215: 395–431.

*Pflügers Archiv*. 2000; 440(5): 667–678. [CrossRef] [PubMed]

*Exp Eye Res*. 1986; 42(6): 559–568. [CrossRef] [PubMed]

*Invest Ophthalmol Vis Sci*. 2009; 50(4): 1791–1800. [CrossRef] [PubMed]

*Annu Rev Fluid Mech*. 2012; 44: 347–372. [CrossRef]

*Adv. Organ Biol*. 2005; 10: 127–148. [CrossRef] [PubMed]

*Prog Retin Eye Res*. 2011; 30(1): 1–17. [CrossRef] [PubMed]

*Fluid Dyn*. 2007; 42(5): 684–694. [CrossRef]

*J Model Ophthalmol*. 2016; 1(2): 29–39.

*Life*. 2020; 10(2): 8. [CrossRef]

*J Theor Biol*. 2018; 456: 233–248. [CrossRef] [PubMed]

*Am J Physiol Cell Physiol*. 1996; 271(3): C703–C720. [CrossRef]

*Am J Physiol Cell Physiol*. 1995; 268(3): C572–C579. [CrossRef]

*Am J Physiol Ren Physiol*. 1996; 270(2): F344–F355. [CrossRef]

*Mathematical Physiology*. New York:Springer; 2009: x1.

*SIAM Journal on Scientific Computing*. 1997; 18(1): 1–22. [CrossRef]

*Technometrics*. 1999; 41(1): 39–56. [CrossRef]

*J Theor Biol*. 2008; 254(1): 178–196. [CrossRef] [PubMed]

*PLOS Comput Biol*. 2020; 16(11): e1008248. [CrossRef] [PubMed]

*Exp Eye Res*. 1996; 62(2): 131–140. [CrossRef] [PubMed]

*Clin Exp Optom*. 2002; 85(6): 335–349. [PubMed]

*J Physiol*. 1989; 418(1): 83–103. [CrossRef] [PubMed]

*Invest Ophthalmol Vis Sci*. 1970; 9(1): 20–29.

*Investig Ophthalmol Vis Sci*. 2007; 48(1): 321–327. [CrossRef]

*Am J Physiol Cell Physiol*. 2001; 280(6): C1521–C1530. [CrossRef] [PubMed]

*J Gen Physiol*. 1967; 50(8): 2061–2083. [CrossRef] [PubMed]

*Invest Ophthalmol Vis Sci*. 2002; 43(6): 1897–1902. [PubMed]

*J R Soc Interface*. 2020; 17(163): 20190735. [CrossRef] [PubMed]

*Mathematical Physiology*, New York:Springer; 2009.

*PloS One*. 2007; 2(1): e144. [CrossRef] [PubMed]

*J Gen Physiol*. 1990; 96(2): 319–344. [CrossRef] [PubMed]

^{+}and Cl

^{−}fluxes across the rabbit ciliary epithelium.

*Invest Ophthalmol Vis Sci*. 1987; 28(3): 445–450. [PubMed]

^{+}-K

^{+}-2Cl

^{−}cotransporter NKCC1 in mammalian epithelial cells.

*J Physiol*. 2010; 588(21): 4089–4101. [CrossRef] [PubMed]

*Am J Physiol Renal Physiol*. 2005; 289(4): F721–F741. [CrossRef] [PubMed]

_{2}permeability of cell membranes is regulated by membrane cholesterol and protein gas channels.

*FASEB*. 2012; 26(12): 5182–5191. [CrossRef]

*Cellular Biophysics I. Transport*. Cambridge, MA: MIT Press; 1996.

*Cellular Biophysics II. Electrical Properties*. Cambridge, MA: MIT Press; 1996.

^{35}

^{35}according to which the flux of ion \(i\) from region \(m\) to region \(k\) can be written as

^{36}

^{37}according to which

^{37}

^{35}and their expressions read

^{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

**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/Pa^{39}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.

^{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.

^{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.

^{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\).