A Mathematical Model of Aqueous Humor Production and Composition

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 K+ and 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.

A queous 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 proteins 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.
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 secre-tory monolayer epithelium. This view is reinforced when one considers the distribution of ion channels and transporters, and of the Na + -K + pump, as shown in Figure 2, where the NPE basolateral membrane has a channel distribution reminiscent of an epithelial apical domain. [7][8][9][10] 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, (A) Ion channels in the CE. PE basolateral membrane: NKCC: co-transporter of Na + /K + /2Cl − ; AE s : Cl − -HCO − 3 exchanger; NHE: Na + -H + exchanger; NBC s : Na + -HCO − 3 cotransporter, presumed to be of stoichiometry 1:2; K + -channels. NPE basolateral membrane: PUMP: Na + -K + ATPase; AE p : HCO − 3 -Cl − exchanger; NBC p : Na + -HCO − 3 cotransporter; K + and 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 symbolsp,s and t j denote the PC and stromal cell membranes and the tight junction. In the direction from region m to region k, J mk i denotes flux of solute i and Q mk the water flux (see main text).
p, between the stroma and the PC that induces a water flux into the PC. The second is an oncotic pressure difference, 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 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: 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 σ s and σ p are reflection coefficients for ions and proteins, respectively. Oncotic and osmotic pressures can be evaluated using van't Hoff's law = RT C, with C being difference in solute concentration (osmolar-ity) 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 Na + -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 superscripts) and the PC membrane (p). We further assume that the tight junctions act as semi-permeable membranes between stroma and the PC (t j). 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: Na + , K + , Cl − , HCO − 3 (bicarbonate), H + , CO 2 , and H 2 CO 3 (carbonic acid). We denote the concentrations of these solutes as C m i (mM), where subscript i denotes the species (i = 1, . . . , 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 nondiffusible 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 mk i (mol/s). For instance, the flux of Na + from the stroma into the cell is named J sc 1 . We consider that the stromal side of the membrane contains co-transporters of Na + , K + and 2Cl − (denoted with NKCC), anion exchangers of Cl − and HCO − 3 (AE s ), 7 Na + -H + exchangers (NHE), 8 co-transporters of Na + and 2HCO − 3 (NBC s ) and K + channels, 18 as shown in Figure 2A. The PC side of the membrane has Na + -K + ATPase (denoted with PUMP), 9 Cl − -HCO − 3 exchangers (AE p ), Na + -HCO − 3 co-transporters (NBC p ), K + and Cl − channels. 18,19 The tight junctions are assumed to be permeable to Na + , K + , Cl − and HCO − 3 . In addition, all membranes and the tight junction are permeable to CO 2 and H 2 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 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 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 mk i , 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 p i , 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 ≈ 0.2 mm/s, 4 so using typical PC length H ≈ 6 mm and a diffusion coefficient D ≈ 2 · 10 −3 mm 2 /s, 17 we can calculate the Péclet number, Pe = U H/D ≈ 600, confirming that advection has a dominant contribution at the PC exit.
The considered solutes interact through the following chemical reactions with k 1 , k −1 , k d , and k h being the reaction rate constants. The first reaction is the dissociation of H 2 CO 3 into the ions HCO − 3 and 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 m i the reaction terms (production minus consumption) of species i in region m, which we model with the mass action law 21 (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 where 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 In the stromal compartment all species concentrations C s i 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 Cell PC with i = 1, . . . , 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

Reaction rates [RR] and other parameters [Other]
Parameter Each section corresponds to a certain parameter group, each having its own marker reported in square brackets, which is used in Appendix 2 for easier referencing. For example, the reaction rates are referred to as RR. In the last section we report permeabilities and intensities of the ion channels used in the model. Left Side: Intensities associated with the fluxes through ion transporters as in (2), expressed in ·10 −6 mol/m 2 /s. Right Side: Permeability of ion channels in ·10 −8 m/s as in (8). In last row we report the number of the mathematical expression of the fluxes where these parameters appear. Justification for the choice of some of these parameters and the corresponding sources are reported in Appendix 2. details). The equations are complemented by electroneutrality in both the cell and the PC In total, we have 18 nonlinear algebraic equations for an equal number of unknowns: concentrations in the cell and the PC: C c i , C p i (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.

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

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 ).
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; Cl − flux has a magnitude compatible with measurements, whereas the Na + flux is overestimated, which is considered further in the Discussion.

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 abovementioned discrepancy with the 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 Na + -K + ATPase, thus the pump intensity obviously influences all the model results. Indeed, reduction of the parameter P 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 Na + -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  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).
vary each parameter within ±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 Na + , K + , Cl − and 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 ±25% and ±75% of the baseline values of Table 1, with no substantive qualitative changes in the results. Figure 3A shows that AH production rate is very sensitive to the permeability to K + (Pp K + ), and Cl − (Pp Cl − ) and to the intensities of the Na + -H + exchanger (P NHE ) and the cotransporter of Na + /K + /2Cl − (P 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 HCO − 3 concentration in the PC drops significantly and 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 K + and Cl − concentration in the PC and a corresponding decrease in Na + and HCO − 3 . 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 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 (Na + , K + , Cl − , HCO − 3 , H + , CO 2 , and H 2 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 B A  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 Cl − flux has a magnitude in line with measurements, whereas the Na + flux is overestimated (Table 2). We note, however, that in the measurement of Na + current, 28 the authors detected only about 25% of 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 · 10 −6 mol/m 2 /s, which is the same order of magnitude as our prediction.
Experimental evidence shows that the inhibition of the Na + -K + ATPase markedly increases the intracellular Na + content, and reduces the intracellular 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 K + and 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 H + and 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 zerodimensional lumped parameter model without the mechanism of local osmosis. of these species across the membrane are modeled using an expression derived from the integration of the diffusion equation across the membrane, 35 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 CO 2 and H 2 CO 3 across each membrane and the tight junction. Ion Channels. Ion channels, such as 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 where factor P is a membrane permeability to this ion (m/s), A is the area of the membrane, φ 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 K + channels on the cell membrane, 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 Na + -K + -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

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 · 10 −3 m. The area of the PC at the boundary with ciliary processes is estimated to be A PC = 2.6 · 10 −5 m 2 . 4 Owing to the extensive folding the area of the ciliary processes is considered to be Ap = As = 6 · 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 O(10 −3 ). P-H2O Experimental estimates of hydraulic conductivity of the CE vary from 1.1 · 10 −11 m/s/Pa 39 to 6.26-9.5·10 −11 m/s/Pa. 16 In this work, we set Ks = 2 · 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 Kp = 10Ks. RR The rates of the reactions in (4) for hydration and dehydration of CO 2 (catalyzed) are taken to be k h = 1.45 · 10 3 1/s and k d = 4.96 · 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 H 2 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 CO 2 can range from about 10 −4 m/s to 10 −2 m/s. In our model we impose Ps CO 2 = 1.5 · 10 −3 m/s, following the value used in Krahn and Weinstein, 20 for the proximal tubule brush border. The permeability to H 2 CO 3 is also taken from the same article, to be Ps H 2 CO 3 = 1.28 · 10 −5 m/s. 20 The permeabilities of the posterior membrane to both CO 2 and H 2 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 PUMP ]: 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 Na + flux was found in isolated cat ciliary body. 29 Thus, Na + flux expressed using the international system of units is in the range of approximately [6.11, 9.32] · 10 −7 mol/m 2 /s. These values might be underestimated (see the Discussion), and we have selected the magnitude of the Na + -K + ATPase to be 6 · 10 −6 mol/m 2 /s.
[Pp Cl − ]: J net Cl is measured to be approximately [3, 8.23] · 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 Cl − of the posterior membrane to be within the range Pp Cl − ∈ [2.97, 8.16] · 10 −8 m/s.
[P t j ]: Assuming that tight junctions are the main players in determining the value of transepithelial resistance, T ER = 0.6 k ·cm 2 , 18 we compute Na + flux through the tight junctions, with the values of concentration from Table 2, to obtain 3.37 · 10 −6 mol/m 2 /s (directed toward stroma). We use this value to estimate the permeability of tight junctions to ions, P t j = 9.22 · 10 −7 m/s. A higher value of 6·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 σ = 1. 42,43 The average charge and valence of fixed nondiffusible species X in the cell is taken to be z X = −1.5.