## Abstract

This paper presents a numerical study on the transport of ions and ionic solution in human corneas and the corresponding influences on corneal hydration. The transport equations for each ionic species and ionic solution within the corneal stroma are derived based on the transport processes developed for electrolytic solutions, whereas the transport across epithelial and endothelial membranes is modelled by using phenomenological equations derived from the thermodynamics of irreversible processes. Numerical examples are provided for both human and rabbit corneas, from which some important features are highlighted.

## 1. Introduction

The cornea is a transparent tissue that is highly specialized to refract and transmit light. Corneal transparency is heavily dependent upon the regulation of normal stromal hydration. Isolated stromal tissues swell freely in isotonic saline and become opaque, while the *in vivo* corneas remain thin and clear. How the limiting membranes control the stromal hydration and thus maintain the hydrophilic stroma in a state of relative deturgescence is one of the most interesting questions in corneal physiology (Ruberti & Klyce 2003).

Modelling of corneal swelling and its interaction with its environment is critical to our understanding of corneal function, particularly when important physiological parameters are refractory to experimental investigation. For the cornea, endothelial membrane transport parameters are difficult to access directly, because the membrane cannot be separated atraumatically from the underlying stromal tissue. Thus, development of an effective stromal transport model could help isolate the transport properties of the endothelial membrane from those of the stroma. Considerable efforts have been made in past decades for developing dynamic models to simulate the whole corneal transport system and to describe how the stroma and limiting layers interact to control stromal hydration (Ruberti & Klyce 2002). These models can be generally divided into two groups. The first group is the transport model (e.g. Stanley *et al*. 1966; Friedman 1973; Klyce & Russell 1979; Ruberti 1999; Li *et al*. 2004) in which the increase or decrease of stromal thickness at a particular point was determined based on the increase or decrease of the water volume at that point. The latter was calculated based on assumed stromal transport models for the flow of ions and ionic solution together with assumed boundary conditions for the flow of ions and ionic solution across epithelial and endothelial membranes. The second group is the triphasic model based on the poroelastic theory (Lai *et al*. 1991) in which the solid matrix is regarded as a linear isotropic elastic material and the corneal swelling is considered to be governed by the transport of ions in the electrolyte within the poroelastic medium controlled by the coupled electrochemical and mechanical processes (Bryant & McDonnell 1998). The volume change of the cornea is thus related to both the fluid flow and the elastic deformation of the corneal matrix. However, though this triphasic model is well grounded with regard to the fundamental processes responsible for transport of fluid and ions in the cornea, the model itself is actually static and is, therefore, incapable of performing dynamic simulations of the corneal response to external osmotic perturbations (Ruberti & Klyce 2002).

The work presented in this paper is the further development of our previous model (Li *et al*. 2004) by considering the interactions between ionic species and ions and ionic solution when they are in transport within the corneal stroma and across the epithelial and endothelial membranes. The transport equations for each ionic species and ionic solution within the corneal stroma are derived based on the transport processes developed for electrolytic solutions, whereas the transport across epithelial and endothelial membranes is modelled by using phenomenological equations derived based on the thermodynamics of irreversible processes.

## 2. Brief review of corneal structure

Geometrically, the human cornea can be approximately regarded as a spherical cap with a variable thickness. The base diameter of the cornea is about 12 mm and the average radius of curvature in the corneal central region is about 7.8 mm. The thickness of the cornea is roughly 0.52 mm at its centre and 0.65 mm in the periphery. The diffusion distance from the periphery to the corneal centre is over 6 mm, which is more than 10 times that from the posterior to the anterior at the centre (0.52 mm). Thus, it is clear that as far as the transport is concerned, the cornea can be modelled as a one-dimensional tissue in the anterior-to-posterior direction.

In the aspect of corneal structure and material, the human cornea is comprised of an outer stratified squamous non-keratinized epithelium, an inner connective tissue stroma with resident keratocytes and a low cuboid endothelium. All of the three layers have a uniform and consistent arrangement throughout the tissue (Gipson 1994).

The epithelium has five to seven cell layers, approximately 40 μm thick and is located at the outside surface of the cornea. The epithelium is extraordinarily regular in thickness over the entire cornea and it has an absolutely smooth, wet, apical surface that serves as the major refractive surface of the eye. The cell layers of the epithelium include three to four outer flattened squamous cells; one to three layers of mid-epithelial cells; and a single layer of columnar basal cells. The functions of the corneal epithelium include providing a barrier between the environment and the stroma of the cornea and, through its interaction with the tear film, forming a smooth refractive surface on the cornea. The epithelial barrier is formed as the epithelial cells move from the basal layer to the surface of the cornea, progressively differentiating until the superficial cells form two layers of flattened cells encircled by tight junctions, which serve as a semi-permeable, high-resistance membrane. This barrier prevents the movement of fluid from the tears into the stroma and also protects the eye from the penetration of tear-borne pathogens. In physiology, the epithelium also has the ability to actively transport chloride ions from the stroma to the tears, which assists the corneal endothelium in the regulation of stromal hydration and thereby contributes to the maintenance of corneal transparency (Klyce & Russell 1979; Klyce & Crosson 1985).

The human corneal stroma is the middle connective tissue layer that is approximately 500 μm thick and forms the bulk of or about 90% of the corneal thickness. The stroma is arranged in three clearly defined layers of extracellular matrix. These include: bordering the epithelium, the 8–10 μm Bowman's layer; the middle lamellar stroma, which comprises by far the major portion of the stroma and adjacent to the endothelium, the 8–12 μm Descemet's membrane, the thickened basement membrane secreted by the corneal endothelium (Gipson 1994). The lamellar stroma is the major layer of the stroma and is comprised of 300–500 sheets formed from flattened bundles of collagen fibrils oriented in a parallel manner. Each lamella is about 2–3 mm broad and about 1.5–2.5 μm thick and is comprised of long collagen fibrils embedded in a matrix containing proteoglycans which is highly hydrated (Maurice 1988; Daxer *et al*. 1998). The mean diameters of the fibrils in the corneal stroma are in the range between 22 and 32 nm (Pinsky & Datye 1991; Doughty & Bergmanson 2005). The stroma meshes with the surrounding scleral connective tissue to form a rigid framework for maintaining intraocular pressure and, thus, alignment of the optic pathways, although it has very weak rigidity in its thickness direction and can swell to several times its normal hydration in the absence of active maintenance by the epithelium and endothelium.

The corneal endothelium is a thin monolayer of hexagonal cells that forms the posterior corneal surface. The thickness of the endothelium is about 4 μm, which is 10 times thinner than that of the epithelium. Compared to the epithelium, which is a semi-permeable, ‘tight’ barrier, the endothelium is a permeable, ‘leaky’ barrier that restricts the movement of water and solutes into the hydrophilic stroma. The primary function of the endothelium is corneal dehydration. For a normal cornea, hydration is maintained at 78% water (by weight). Since fluid and solutes are continuously driven into the stroma by stromal swelling and imbibition pressures, the maintenance of corneal thickness and transparency is dependent on the active removal of fluid that leaks into the stroma. It has been thought that this is due to the attribution of the corneal endothelium, that is the so-called pump-leak hypothesis. Although the exact mechanism of the endothelial pump is not completely understood, the effect of endothelial ion transport resulting in net fluxes of sodium and bicarbonate ions from stroma to aqueous humour has been confirmed in experiments (Edelhauser *et al*. 1994). Corneal endothelial carbonic anhydrase, which catalyses the conversion of carbon dioxide and water into bicarbonate and hydrogen ions, is believed to provide an important source of bicarbonate for the endothelial pump (Bonanno 2003).

## 3. Electrolytic transport model in corneal stroma

The transport model described below is similar to those proposed by Klyce & Russell (1979) and Li *et al*. (2004), which is a three-compartment, one-dimensional model used to predict the change of stromal thickness in response to the concentration changes applied to corneal epithelium and/or endothelium (see figure 1). However, since the solution considered here is a multi-component solution, the governing equations for both ions and ionic solution are derived based on the transport processes in electrolytic solutions, which is different from those early developed models.

The governing equations describing the transport of ionic solution and ionic species in corneal stroma can be derived in terms of the mass conservation (Li *et al*. 2004)(3.1)(3.2)where *W*_{f} in g g^{−1} is the mass of the solution in unit mass of dry tissue material, *t* in s is the time, *ρ*_{f} in g cm^{−3} is the density of the solution, *J*_{f} in cm^{3} (s cm^{2})^{−1} is the volume flux of the solution flowing in unit time through unit area in the thickness direction, *ξ* in cm is the thickness coordinate of the dry tissue material, *ρ*_{d} in g cm^{−3} is the density of the dry tissue material, *C*_{k} in mole cm^{−3} is the concentration of species *k*, *J*_{k} in mole (s cm^{2})^{−1} is the molar flux of species *k* flowing in unit time through unit area in the thickness direction. The index *f* stands for ionic solution and index *k* (*k*=1, 2, …, *M*) represents species *k* and *M* is the total number of species involved in the solution. Equations (3.1) and (3.2) represent the well-known mass conservation in the unit mass of dry tissue material for the solution and for species *k*, respectively.

The volume flux of the solution can be simply assumed to be the contribution of convection, which is given by(3.3)where *v* in cm s^{−1} is the velocity of convective flow. The molar flux of species *k* can be expressed in the form of the Nernst–Planck equation as follows (Newman 1973):(3.4)where *z*_{k} is the charge number of species *k*, *D*_{k} in cm^{2} s^{−1} is the diffusion coefficient of species *k*, *F*=96 485 C mole^{−1} is the Faraday constant, *R*=8.314 J (mole K)^{−1} is the universal gas constant, *T* in K is the absolute temperature, *ϕ* in V is the electrostatic potential, *x* in cm is the transient thickness coordinate, which has the following relationship with the thickness coordinate *ξ* in dry tissue material(3.5)

Note that for charged ions the molar fluxes of species should satisfy the conservation equation of current, which is expressed by(3.6)where *I* in A cm^{−2} is the current density. Substituting equation (3.4) into (3.6) and noting that because of the electro-neutrality, yields(3.7)

Substituting equations (3.5) and (3.7) into (3.4) yields(3.8)where is the new diffusion coefficient defined as follows:(3.9)

According to the Darcy law, the velocity of the convective flow in the corneal stroma is negatively proportional to the local gradient of hydrostatic pressure, i.e.(3.10)where *K*_{μ} in cm^{4} (s dyn)^{−1} is the flow conductivity coefficient and *P* in dyn cm^{−2} is the local hydrostatic pressure, which is calculated as follows (Fatt & Goldstick 1965)(3.11)where IOP=2.67×10^{4} dyn cm^{−2} is the intraocular pressure, *S* in dyn cm^{−2} is the stromal swelling pressure, *γ*=2.41×10^{6} dyn cm^{−2} is the empirical constant and *W*_{w} is the stromal hydration which is defined as the mass of water in unit mass of dry tissue material. Substituting equations (3.10) and (3.11) into (3.3) yields(3.12)

For dilute solutions for which it can be assumed that , the following equations can be obtained to determine *W*_{w} and *ρ*_{f}:(3.13)where *C*_{w} in mole cm^{−3} is the concentration of water (solvent), *υ*_{w}=18 cm^{3} mole^{−1} is the partial molar volume of water, *m*_{w}=18 g mole^{−1} is the molar mass of water, *m*_{j} in g mole^{−1} is the molar mass of species *j*, *ρ*_{w}=1 g cm^{−3} is the density of water.

For given initial conditions and boundary conditions, equations (3.1) and (3.2) can be used to determine *W*_{f} and *C*_{k}, in which *J*_{f} and *J*_{k} are defined by equations (3.8) and (3.12), *W*_{w} and *ρ*_{f} are defined by equation (3.13).

## 4. Electrolytic transport model at epithelial and endothelial membranes

Note that the flux expressions given by equations (3.8) and (3.12) are only applicable to the corneal stroma. For the epithelial and endothelial membranes, the flux expressions for ionic solution and ionic species can be obtained from the theory of irreversible thermodynamics (Kedem & Katchalsky 1958; Li 2004). By assuming that the epithelial and endothelial membranes are the infinitely thin uniform film, the following volume flux for ionic solution and molar flux for species *k* can be obtained(4.1)(4.2)where *L*_{p} in cm^{3} (dyn s)^{−1} is the hydraulic conductivity coefficient, Δ*P* in dyn cm^{−2} is the hydrostatic pressure difference on the two sides of the membrane, *σ*_{k} is the reflection coefficient of species *k*, Δ*C*_{j} in mole cm^{−3} is the concentration difference of species *j* on the two sides of the membrane, Δ*ϕ* in V is the electrostaic potential difference on the two sides of the membrane, *ω*_{k} in mole (dyn s)^{−1} is the permeability coefficient of species *k*, *J*_{ak} in mole (s cm^{2})^{−1} is the active pump rate of species *k*. The sign convention used in equations (4.1) and (4.2) is defined as follows. Net fluxes of solution and solutes from outside to inside are taken as positive. The difference in solute concentration across the membrane is calculated by the outside concentration minus the inside one. Differences of hydrostatic and osmotic pressures are defined in a similar manner.

Again, the fluxes of species across the membrane should satisfy equation (3.6) from which the electrostatic potential gradient across the membrane can be eliminated. Thus, equations (4.1) and (4.2) become (Li 2004)(4.3)(4.4)in whichwhere are the new membrane coefficients, and are the new volume flux of solution and new molar flux of species *k*, which are generated by the active pump rates. Note that these new coefficients are concentration dependent. The reason for this is because the fluxes can be generated by the hydrostatic pressure difference, the concentration difference or the electrostatic potential difference. Equations (4.3) and (4.4) are obtained by eliminating the electrostatic potential item using the condition of electro-neutrality. This is equivalent to absorb the electrostatic potential item into the membrane coefficients to form the new membrane coefficients. It is also interesting to note from equation (4.3) that the active pump rate may have direct contribution to the volume flux if . In past models, the ionic solution was assumed as an equal binary solution for which . Therefore, the counterbalance was actually provided by the concentration gradient, which, of course, is created by the active pump. However, in the present multi-component model, the active pumps not only provide the concentration gradient, but also directly contribute to the volume flow as is demonstrated in equation (4.3).

Equations (4.3) and (4.4) provide the flux boundary conditions of the mass transport equations (3.1) and (3.2) for the ionic solution and ionic species. If the thickness changes of the epithelial and endothelial membranes are negligible then the change of the corneal thickness can be calculated purely based on the change of the corneal stroma thickness as follows (Hedbys & Mishima 1966):(4.5)where Δ*h* in cm is the change of corneal thickness, *ξ*_{max} in cm is the thickness of the dry corneal stroma, *ρ*_{f0} in g cm^{−3} is the initial density of the solution, *W*_{f0} and *W*_{w0} are the initial masses of solution and water in the unit mass of dry tissue material, respectively.

## 5. Numerical results

For given initial conditions (such as initial ionic concentrations and initial hydration in stroma) and boundary conditions (ionic concentrations and hydrostatic pressures in tears and aqueous humour), one can solve partial differential equations (3.1) and (3.2) using finite element methods to obtain the hydration distribution, *W*_{w}, and then use equation (4.5) to obtain the thickness variation, Δ*h*. Figure 2 shows the comparison of the present model predictions and Fischbarg's experimental data (1973) obtained from the tests of rabbit corneas. Since in Fischbarg's experiments the bare anterior stroma (epithelium) was covered with silicone oil and the osmotic perturbations of NaCl solution were applied only to the endothelial surface, the fluxes in the numerical simulation defined at the epithelium for both solution and solutes were assumed as zero, whereas the fluxes at the endothelium were assumed to be described by equations (4.3) and (4.4). The parametric values employed in the numerical simulation for this equal binary solution are the same as those used in our previous paper (Li *et al*. 2004) and thus are not given again here. It can be seen from the figure that the predicted thickness response is in good agreement with the experimental data.

Under normal *in vivo* physiological conditions, there is no net fluid transport since the fluid into the cornea driven by the stromal swelling pressure is exactly counterbalanced by the active secretion of fluid out of the cornea. This is the so-called ‘pump-leak’ mechanism hypothesis for maintenance of corneal hydration and transparency (Maurice 1972). It has been several decades since the hypothesis was first made. However, so far it is still not clear about which species are responsible for generating such a fluid pump leak and which species are capable of providing active pump fluxes. In literature, it has been suggested that the endothelial pump leak involves complicated ionic exchanges (such as Na^{+}/H^{+} and Cl^{−}/HCO_{3}^{−}) and ion cotransport (such as Na^{+}/HCO_{3}^{−} and K^{+}/HCO_{3}^{−}, Na^{+}/K^{+}/2Cl^{−}) (Lane *et al*. 1997; Bonanno 2003). However, there is no quantitative model available to describe how these ionic species are interrelated and interact on each other (Kuang *et al*. 2004). Therefore, in the following numerical examples, the active pump rates are simply assumed.

The dominant ionic species found in tears and in aqueous humour are Na^{+}, K^{+}, Cl^{−} and HCO_{3}^{−}. It is probably reasonable to assume these species of constant concentrations, because the volumes of tears and aqueous humour are much greater than the stromal volume. The initial concentrations of Na^{+}, K^{+}, Cl^{−} and HCO_{3}^{−} are assumed to be uniform in the corneal stroma, the values of which together with other parametric values used in the following numerical study are given in table 1. The parametric values, which are either independent of the ionic species or only related to chloride ions, are taken from literature (Klyce & Russell 1979). The parametric values, which are related to other than chloride ions, are calculated based on the relative ratios of their diffusion coefficients to the chloride ion's diffusion coefficient.

Figure 3 shows the responses of the stromal thickness to the shocks of 15 mOsm NaCl, NaHCO_{3}, KCl and KHCO_{3} hypertonic solutions applied (after the steady state has been reached) to the endothelial surface for an hour duration, while the epithelial surface is either blocked (figure 3*a*) or left intact (figure 3*b*). The results show that the responses of the stromal thickness are different to the solutions of different solutes although in terms of the concentration they are actually identical. The greatest change in thickness is found to the solution of NaHCO_{3} for which the permeability coefficients are the smallest and the reflection coefficients are the largest. While the least change is found to the solution of KCl for which the permeability coefficients are the largest and the reflection coefficients are the smallest. The influence of the epithelial flow on the stromal thickness response to endothelial perfusion appears remarkable. When the epithelial surface is left intact and functioning, the thickness reverse becomes much quicker and also much greater.

Figure 4 shows the responses of the stromal thickness to the same shocks but applied at the epithelial surface, while the endothelial surface is either blocked (figure 4*a*) or left intact (figure 4*b*). The results show that the stromal thickness response to the epithelial perfusion is rather slow compared to those shown in figure 3. Significant reduction in stromal thickness changes is found when the endothelium is left intact and functioning. This indicates that the endothelium is able to alleviate the stromal thickness change when the corneal is subjected to a perfusion disturbance from the epithelial side.

The effect of the active pump rate provided by ionic species in the epithelial layer is examined by using zero pump rates. Figure 5 shows the results in which the ionic pump rate for Cl^{−} in the epithelial layer is set to zero. It is found that the stromal thickness response pattern in figure 4 is similar to that shown in figure 5. It is seen that when the endothelial surface is left intact there is almost no difference between figures 4*b* and 5*b*. The reason for this is probably partly because the value of the pump rate in the epithelial layer (figure 4) is small and partly because both the hydraulic conductivity and permeability coefficients of the epithelial membrane are very small compared to those of the endothelial membrane. The latter leads the stroma to have a slow response to the shock applied to the epithelial surface. This is why the stroma has different response patterns to the shocks applied to the epithelial and endothelial surfaces. However, when the endothelial surface is blocked, the pump rate does have some influence on the thickness responses, particularly on the reversal of the response curves. The comparison of figures 4*a* and 5*a* shows that when the active pump rate provided by Cl^{−} is considered more thickness reduction can be achieved.

In contrast, the effect of the active pump rate provided by ionic species in the endothelial layer appears very significant. Figure 6 shows the responses of the stromal thickness to the same shocks applied to the endothelial surface while the endothelial ionic pump rate is set to zero. The figure shows that, when the bicarbonate pump is tuned off, the stromal thickness responses are about double in magnitude when compared to the same shock with the pump on (see figure 3). Also, when the pump is tuned off, there is almost no recovery in the response curves after the reversal. The former is because the corneas have different states when applying shocks, although their initial states are the same (see table 1). Since the shocks are applied after the same certain time from the same initial state, the ionic concentrations and swelling pressure are different for the two cases because of the influence of the active pump. While the latter is simply due to the function of the swelling pressure that keeps water flowing into the stroma, which is demonstrated by the response curve that has the same slope at places before and after the shock.

Parametric study of the dynamic transport model on the cornea hydration has shown that, though reasonable in magnitude and in general shape, the model with using a standard set of fixed membrane transport coefficients cannot adequately describe the actual corneal response found in experiments, unless the membrane coefficients are allowed to change at the time of the perturbation and reversal (Ruberti & Klyce 2003). It is interesting to note from the present results, figure 3, for example, that the stromal responses are different to the shocks with different solutes, although each has the same concentration. Assume that there exist ionic exchanges between Na^{+} and H^{+} and between Cl^{−} and HCO_{3}^{−} in the endothelial layer as it is suggested (Bonanno 2003). This implies that the initial shock of NaCl at the outside surface of the layer could become NaHCO_{3} and HCl when it comes to the inside surface of the layer or *vice versa*. If this were indeed happened, the variation of the corneal thickness could be modulated by 30%.

## 6. Conclusions

A computational model of the corneal hydration has been proposed for simulating the response of the corneal thickness to the perfusion with hypotonic or hypertonic solutions. The transport equations for each ionic species and ionic solution within the corneal stroma are derived based on the transport processes developed for electrolytic solutions, whereas the transport across epithelial and endothelial membranes is modelled by using phenomenological equations derived from the thermodynamics of irreversible processes. The influence of the flow across the epithelium and endothelium on the stromal thickness response has been investigated. The role of the ionic active pumps in epithelial and endothelial layers has also been examined. The present numerical examples are based on a number of assumed parametric values which are not available. The encouraging results found with present model suggest that further development of theoretical models and experimental methods that can be used to determine these parametric values will be very useful. The present model provides a starting point to simulate the complicated transport process of the human cornea involving ionic cotransports and ionic exchanges taking place in the epithelial and/or endothelial layers.

## Footnotes

- Received May 5, 2005.
- Accepted September 5, 2005.

- © 2005 The Royal Society