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 \(\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 proteins 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,; it flows from the PC into the AC through the pupil and exits the eye via the conventional and uveoscleral pathways. 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., 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).

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 - pump, as shown in Figure 2, where the NPE basolateral membrane has a channel distribution reminiscent of an epithelial apical domain.–

Figure 2.

(A) Ion channels in the CE. PE basolateral membrane: nkcc: co-transporter of //2; : - exchanger; nhe: - exchanger; : - cotransporter, presumed to be of stoichiometry 1:2; - channels. NPE basolateral membrane: : - ATPase; : - exchanger; : - cotransporter; and channels. (B) A sketch of the compartmental model with water and ion fluxes (not to scale). The symbols , and denote stroma, intracellular space and PC respectively. The symbols and denote the PC and stromal cell membranes and the tight junction. In the direction from region to region , denotes flux of solute and 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. 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, , between the stroma and the PC that induces a water flux into the PC. The second is an oncotic pressure difference, , 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., The main mechanism is thought to be osmotic pressure difference 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, is water flux per unit surface area (with the dimensions of velocity), is the hydraulic conductivity of the CE, and and are reflection coefficients for ions and proteins, respectively. Oncotic and osmotic pressures can be evaluated using van't Hoff's law , with being difference in solute concentration (osmolarity) across the cell layer, the universal gas constant, and the absolute temperature. Equation (1) has been used in several modeling works to evaluate AH production., 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. proposed a more detailed model that takes into account the distribution of ion channels on the cell membrane (with the exception of the - 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 ), the cellular compartment, which encompasses both NPE and PE cells () and the PC (). The compartments are separated by the stromal membrane (denoted with superscript ) and the PC membrane (). We further assume that the tight junctions act as semi-permeable membranes between stroma and the PC (). 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. 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: , , , (bicarbonate), , , and (carbonic acid). We denote the concentrations of these solutes as (mM), where subscript denotes the species (, enumerated in the order above) and superscript the region . We denote their valence and charge with . The cell contains fixed non-diffusible charged solutes, which we denote with , and their concentration and valence are denoted as and . The electric potential in region is denoted with .

Solutes are transported across the membranes of the domain via ion channels, co-transporters, and exchangers and we indicate a flux of a solute from region to region with (mol/s). For instance, the flux of from the stroma into the cell is named . We consider that the stromal side of the membrane contains co-transporters of , and 2(denoted with nkcc), anion exchangers of and (), - exchangers (nhe), co-transporters of and 2 () and channels, as shown in Figure 2A. The PC side of the membrane has - ATPase (denoted with ), - exchangers (), - co-transporters (), and channels., The tight junctions are assumed to be permeable to , , and . In addition, all membranes and the tight junction are permeable to and .

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 can be formally written as with being a non-linear, dimensionless function that involves solute concentrations and potentials on both sides of the membrane, the intensity of channel transport (mol/m/s), and (m) the surface area of the considered membrane.

The flux of solute can also be driven across the membrane between regions and by electrodiffusion through ion channels (or by diffusion for non-charged solutes), and has the generalized representation where is a membrane permeability to this solute (m/s) and is a function that has dimensions of concentration (mM), with a dependence on solute concentrations and potentials either side of the membrane (see the Appendix for details).

The flux of each solute from region to region , , is the sum of fluxes through all channels and transporters on that membrane, which involve solute . We also assume that at the outlet of the PC the solute is only driven by advection (, with [m/s] being the rate of aqueous production) and that diffusive flux can be neglected. Indeed, the velocity in thin iris-lens channel is about  mm/s, so using typical PC length  mm and a diffusion coefficient  mm/s, we can calculate the Péclet number, , confirming that advection has a dominant contribution at the PC exit.

The considered solutes interact through the following chemical reactions with , , , and being the reaction rate constants. The first reaction is the dissociation of into the ions and , 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). CA II and IV are present in the NPE cells. We denote with the reaction terms (production minus consumption) of species in region , which we model with the mass action law (see Appendix).

Water is transported across each membrane by osmosis and we denote the water flux from region to region with (m/s) and it has the following generic expression where is the difference in osmolarity between the two sides of the membrane, (m/s/Pa) is hydraulic conductivity, and 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 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 with . 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

In total, we have 18 nonlinear algebraic equations for an equal number of unknowns: concentrations in the cell and the PC: , (14 in total), , water flux , and potentials in the cell and the PC , . 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.

Model parameters are reported in Table 1 and discussed in detail along with other model parameters in the Appendix.

Table 1.

Model Parameters.

Geometrical parameters [GP]
Parameter	L (m)		H (m)		A^s˜,p˜ (m^2)		A^tj (m^2)		A^PC (m^2)
Value	10^-5		6·10^-3		6·10^-4		6·10^-7		2.6·10^-5
Hydraulic conductivities [P-H2O] and permeabilities to CO_2 and H_2CO_3 [P-C]
Parameter	K^s˜,tj (m/s/Pa)	K^p˜ (m/s/Pa)		Pco_2s˜,tj (m/s)		Pco_2p˜ (m^2)		Ph_2co_3s˜,tj (m/s)		Ph_2co_3p˜ (m/s)
Value	2·10^-11	2·10^-10		1.5·10^-3		1.5·10^-2		1.28·10^-5		1.28·10^-4
Reaction rates [RR] and other parameters [Other]
Parameter	k_d (1/s)	k_h (1/s)		K_d=k_1/k_-1 (1/mM)		T (K)		σ		z_X
Value	4.96·10^5	1.45·10^3		5.3		310		1		-1.5
Ion channels [P-ions]
Parameter	P_pump	P_nkcc	P_ae_s	P_ae_p	P_nbc_p	P_nbc_s	P_nhe	PK^+s˜	PK^+p˜	PCl^-p˜	P^tj
Value	6	1	4	0.4			3.4	5	30	6	600
Equation	(15)	(9)	(12)	(13)	(11)	(10)	(14)	(8)

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  mol/m/s. Right Side: Permeability of ion channels in  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.

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. and adapted to biological systems by Marino et al. 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.

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 ().

Table 2.

Concentrations of Ions, pH, and the Potential, in the Cell, Stroma, and the PC.

	Experimental measurements26,27,18,7
	[Na^+] (mM)	[K^+] (mM)	[Cl^-] (mM)	[HCO3-] (mM)	pH	V (mV)
Stroma	148	4	107	26		0
Cell	15 ± 3	162 ± 14	46 ± 5			−70
PC	152	3.9	131	22		−1
	Model - baseline
	[Na^+] (mM)	[K^+] (mM)	[Cl^-] (mM)	[HCO3-] (mM)	pH	V (mV)
Stroma	150	5	130	25	7.42	0
Cell	17.8	154.4	45	26.9	7.46	−75.6
PC	151.8	4.3	126.9	28.6	7.49	−1.56
	Model - baseline, no CA
	[Na^+] (mM)	[K^+] (mM)	[Cl^-] (mM)	[HCO3-] (mM)	pH	V (mV)
Cell	5.35	165.3	59.7	18.3	8.78	−73
PC	135	20.3	145	10	7.04	−22

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; flux has a magnitude compatible with measurements, whereas the 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.

Fluxes	Na^+	K^+	Cl^-	HCO3-	H_2O
Units	µmol/m^2/s				m^3/s
Experiments	0.93228		[2.8,7.86]29,30		5·10^-11 3
Model baseline	7.34	0.21	6.16	1.38	2.91·10^-11
Model baseline, no CA	3.61	0.54	3.87	0.28	1.6·10^-11

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 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 - ATPase, thus the pump intensity obviously influences all the model results. Indeed, reduction of the parameter (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 - 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 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. We study their effect on AH production rate (Fig. 3A and on the concentrations of , , and 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 and of the baseline values of Table 1, with no substantive qualitative changes in the results.

Figure 3.

Total sensitivity index for (A) production rate 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. The vertical dashed line separates the intensities of ion channels, in Equation (2), from the permeabilities in Equation (3).

Figure 3A shows that AH production rate is very sensitive to the permeability to (), and () and to the intensities of the -exchanger () and the co-transporter of //2 (). 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 and in reaction (4) by a factor of .

Results for baseline parameters except that CA is inhibited are reported in Tables 2 and 3 Model baseline - no CA. The concentration in the PC drops significantly and 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. The model also predicts an increase of and concentration in the PC and a corresponding decrease in and .

Figure 4.

Change of (A) production rate 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 goes out of the cell and flux through nhe is positive if 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 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.,

We account for the presence of seven different solutes (, , , , , , and ) 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.

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 flux has a magnitude in line with measurements, whereas the flux is overestimated (Table 2). We note, however, that in the measurement of current, the authors detected only about 25% of channels, which suggests that the measured flux may be about four times higher than that reported in Table 2, resulting in the value of  mol/m/s, which is the same order of magnitude as our prediction.

Experimental evidence shows that the inhibition of the - ATPase markedly increases the intracellular content, and reduces the intracellular content, 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 and 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. 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%. 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. explained the effect of CA inhibition on AH formation through the reduction of and 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. 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. This mechanism was found to be by far the dominant mechanism inducing water transport across the retinal pigment epithelium,, 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.