## Abstract

A structural model of the *in vivo* cornea, which accounts for tissue swelling behaviour, for the three-dimensional organization of stromal fibres and for collagen-swelling interaction, is proposed. Modelled as a binary electrolyte gel in thermodynamic equilibrium, the stromal electrostatic free energy is based on the mean-field approximation. To account for active endothelial ionic transport in the *in vivo* cornea, which modulates osmotic pressure and hydration, stromal mobile ions are shown to satisfy a modified Boltzmann distribution. The elasticity of the stromal collagen network is modelled based on three-dimensional collagen orientation probability distributions for every point in the stroma obtained by synthesizing X-ray diffraction data for azimuthal angle distributions and second harmonic-generated image processing for inclination angle distributions. The model is implemented in a finite-element framework and employed to predict free and confined swelling of stroma in an ionic bath. For the *in vivo* cornea, the model is used to predict corneal swelling due to increasing intraocular pressure (IOP) and is adapted to model swelling in Fuchs' corneal dystrophy. The biomechanical response of the *in vivo* cornea to a typical LASIK surgery for myopia is analysed, including tissue fluid pressure and swelling responses. The model provides a new interpretation of the corneal active hydration control (pump-leak) mechanism based on osmotic pressure modulation. The results also illustrate the structural necessity of fibre inclination in stabilizing the corneal refractive surface with respect to changes in tissue hydration and IOP.

## 1. Introduction

In the simplest terms, the cornea is a fibre-reinforced fluid membrane, which resists the intraocular pressure (IOP) applied on its internal boundary. Its tensile strength derives from the three-dimensional organization of collagen fibres and its bulk properties derive from the interfibrillar fluid pressure. The fluid pressure depends on both the IOP and the tissue osmotic pressure. The osmotic pressure causes the cornea to have a tendency to swell by imbibing water from the adjacent anterior chamber of the eye [1,2]. In the *in vivo* cornea, the osmotic pressure is modulated by active ionic transport processes as a means to control the level of tissue hydration, which is important for transparency and also affects the mechanical behaviour of the tissue (cf. bending of a swollen and non-swollen cornea). Ideally, the tissue fluid pressure, swelling effects (including active modulation processes) and the local interaction of collagen fibres with the swelling tissue, should be accounted for in any structural model based on first principles. Some important modelling concepts towards this end have already been introduced, as described below. However, no fully three-dimensional and comprehensive model has yet been presented. It is noteworthy that current finite-element-based models for structural analysis of the cornea treat the interfibrillar fluid as an incompressible or nearly incompressible elastic solid [3–6]. While this approach is convenient, it cannot describe swelling behaviour or (steady-state) bulk compressibility. The goal of this work is to present a mathematical description of the corneal stroma as an electrolyte gel, characterizing the fluid pressure, swelling behaviour and collagen-swelling interaction. We show that such a model can extend predictive accuracy in corneal biomechanics.

The cornea has five principal layers which, from anterior to posterior, are the epithelium, Bowman's, stroma, Descemet's and endothelium. Bowman's and Descemet's layers are basement membranes for the cellular epithelium and endothelium layers, respectively. The stroma occupies 90% of the cornea's thickness and like other soft, highly hydrated and charged collagenous tissues such as cartilage and intervertebral disc, is a polyelectrolyte gel consisting of a mixture of interacting fluid, solid and ionic phases. Water is the principal component of stroma which saturates the collagen solid phase, solvates the ionic phase and accounts for about 78% of the cornea by weight [7]. A small portion of the water is cellular or bound to the stromal collagen [8], but much of the water is free to flow within the stroma in response to gradients in fluid pressure. The solid phase of corneal stroma comprises a flexible collagen network, organized as fibrils within lamellae and associated proteoglycans (PGs). Stromal PGs have sulfated linear side chains of negatively charged disaccharide units called glycosaminoglycans (GAGs) that are covalently bound at one end to the PG core protein [9,10]. At normal pH, the stromal fixed charge is almost entirely due to GAG ionization [1]. The ionic phase includes dissolved salts, primarily Na^{+} and Cl^{−}, and metabolites such as (lactate ion) and (bicarbonate ion).

The stromal mobile ions interact electrostatically with the GAG fixed charges and form cloud-like distributions giving rise to osmotic pressure. The aqueous humour filling the anterior chamber produces the IOP, acts as an ionic bath for the stroma and provides a reservoir of water that is available for exchange with the stroma. Water transport across the permeable endothelial layer is driven by the osmotic pressure difference between the stroma and the aqueous humour. When the stromal ionic concentration exceeds the ionic concentration in the aqueous humour, a positive osmotic pressure difference exists and the stroma will swell by inflow of water from the anterior chamber until the ionic concentrations are equalized or the expansion is restrained by some agent. By contrast, the anterior stroma is sealed by the epithelial cellular layer, which is nearly impermeable to water and ions (although permeable to O_{2} and CO_{2}) [11]. In the *in vivo* cornea, the endothelial layer supports active molecular mechanisms that produce an outward flux of ions from the stroma into the aqueous humour. Our analysis herein confirms that such active pumping reduces the stromal ionic concentrations and hence lowers the osmotic pressure [1]. The regulatory system and details of the molecular mechanisms responsible for active ion transport, which requires Na^{+}, K^{+}, ATPase and carbonic anhydrase activity to transport Cl^{−} and possibly Na^{+}, are not yet fully understood [12].

From a macroscopic perspective, the control of hydration has been explained by a ‘pump-leak’ theory [12–14]. In the pioneering work of [15], based on the phenomenological membrane transport theory of Kedem–Katchalsky (KK) [16], non-equilibrium fluid and ionic fluxes through the corneal stroma and endothelium were modelled and time-dependent swelling resulting from osmotic perturbations were predicted. This work, and the related study by Ruberti & Klyce [17], provides a theoretical confirmation that active endothelial ion transport produces an osmotic gradient across the endothelium which can modulate water flow and swelling. Li & Tighe [18] later proposed an extended model which considers multiple ionic species and Leung *et al.* [19] have further extended the theory to include metabolites (e.g. ) and metabolic reactions in a steady-state model. These models employ a one-dimensional representation of the cornea, do not consider stromal fixed charge or collagen interaction with swelling, and the swelling (i.e. osmotic) pressure is not derived, but assumed as an empirical function of tissue hydration based on measurements by Hedbys & Dohlman [20]. In the important extension by Bryant & McDonnell [21], a spherically symmetric (one-dimensional) steady-state model based on the triphasic theory of Lai *et al*. [22] was combined with KK theory for passive and active transport across the endothelium. This model describes ionic interaction with fixed charges, fluid interaction with a solid elastic phase and ideal Donnan osmotic pressure. While greatly simplifying the stromal elasticity, this model does capture the full range of interactions in an elegant triphasic framework. The above models are complex, and it is not surprising that they do not agree for all predictions, for example, on the sign of the stromal fluid pressure. In this work, we build on these foundational efforts but take an alternative modelling approach.

In the proposed model, conditions of thermodynamic equilibrium are assumed to hold, rendering water and ionic fluxes time-independent. Such an approach has utility when the long-time, steady-state response of the cornea to disease processes or surgical alteration is desired. While modelling based on steady-state conditions does not allow the time course of swelling to be described, it does avoid some of the high complexity of non-equilibrium multiphasic theory and the analytical simplicity of the resulting theory leads to a practical and effective theory for finite-element-based general structural analysis of the living cornea.

Modelled as an electrolyte gel, the free energy is taken to be additively decomposed into various components which characterize the behaviour of the tissue under general deformations [23,24]. A key ingredient in this approach is the statement of electrostatic free energy. For the *ex vivo* cornea with no active ionic transport, we use the mean-field approximation of the Helmholtz free energy for a binary electrolyte which measures the energy of the GAG-based fixed charge, the osmotic energy of mobile ions in a Boltzmann distribution and the dielectric free energy [25]. To account for active endothelial ion transport in the *in vivo* cornea, we consider thermodynamic equilibrium expressed in terms of ionic and water electrochemical potentials and show that the stromal mobile ions satisfy a modified Boltzmann distribution that depends on the active ionic fluxes and endothelial ionic permeability. The theory easily takes into account the important volume exclusion effects associated with collagen and keratocyte populations [23] and the electrostatic free energy corresponding to the modified Boltzmann ionic distribution is readily found and used to predict the osmotic pressure in the *in vivo* cornea. Under certain assumptions on the electrostatic potential, we are able to find analytical approximations for the electrostatic free energy density and osmotic pressure for the *in vivo* cornea. The electrostatic free energy functional is convex in the volume dilation and provides the key variational ingredient for the finite-element formulation. The approach can be easily implemented within a standard finite-element framework using only the displacement field.

Turning to collagen-swelling interaction, the influence of collagen architecture on tissue swelling is readily demonstrated in the *ex vivo* cornea. When a sample of human corneal stroma is immersed in deionized water, the tissue undergoes extreme swelling to many times its original thickness by absorbing water from the surrounding bath. However, the observed swelling is far from uniform, reflecting the interaction of local swelling with the collagen network. The swelling is observed to occur primarily in the thickness direction and is most pronounced in the posterior region of the stroma [26–29]. An analogous situation also occurs in the *in vivo* swollen cornea [30]. This observed rigidity of the anterior cornea with respect to swelling under extreme hydration states has been ascribed to the specific architecture of the collagen fibres [29,31].

The stroma consists of 200–500 superposed sheet-like lamellae. Meek *et al*. [32] quantified the orientation of lamellae when viewed in the plane perpendicular to the optical axis using X-ray diffraction. That study, and numerous others since, described superior–inferior and nasal–temporal preferential collagen orientations in the central cornea and a circumferential orientation near the limbus. These data have been employed directly or indirectly in a number of continuum mechanics-based models to characterize the elastic anisotropy induced by this ‘in-plane’ description of collagen architecture, e.g. [3–6]. Abahussin *et al.* [33] have also performed X-ray diffraction studies on isolated third-thickness central stromal samples and showed that the percentage of lamellae exhibiting preferential orientation is greatest in the posterior stroma and reduces towards the anterior.

To investigate the inclination of lamellae as seen in corneal cross-sections, Morishige *et al*. [34] used second harmonic-generated (SHG) image stacks to reveal lamellae exiting Bowman's layer at steep angles (average 23° maximum 43°) and extending up to 120 µm into the stromal depth. It has long been thought that lamella inclination and insertion into Bowman's layer provides mechanical stabilization of the refractive surface [26–29]. High-resolution macroscopic imaging, which combines many SHG image stacks to produce entire limbus-to-limbus stromal three-dimensional cross-sections, has allowed some lamellae to be traced over millimetres, even with branching and led to the discovery of ‘bow string’ fibres that enter and exit Bowman's layer [35,36]. Recently, Winkler *et al*. [37] processed SHG images to quantify inclination angle distributions in the anterior stroma and found that lamellae inclination angles satisfy a Gaussian distribution with inclination angles being maximum at the anterior stroma and reducing linearly towards the posterior.

In this work, we have synthesized X-ray scattering intensity data and SHG image processing to describe the spatially varying, three-dimensional orientation of lamellae at every point in the human stroma, including through the depth. The raw data are still incomplete in some regions and extrapolation has been necessary, as described in §5. Nevertheless, the most significant features of collagen architecture are now most probably accounted for, allowing a comprehensive description of the elastic anisotropy of the stromal collagen network and regional collagen-swelling interaction.

The proposed model is implemented in a finite-element framework and employed to predict free and confined swelling of stroma in an ionic bath. For the *in vivo* cornea, the model is used to simulate corneal hydration changes which occur as a result of increasing IOP, such as occurs in glaucoma. Corneal swelling in Fuchs' corneal dystrophy is employed as an illustration of how the model can be adapted to describe certain pathological conditions. Whenever the cornea is altered by the surgical removal of tissue, the tissue fluid pressure will change and a swelling response can be anticipated. We analyse the biomechanical response of the *in vivo* cornea to a typical LASIK surgery for myopia and predict, for the first time, how the fluid pressure and tissue swelling will respond. As a validation of the model, the predicted biomechanical response of the cornea (based on elastic deformation and swelling) is compared to extensive clinical data.

## 2. Electrolyte gel in thermodynamic equilibrium

### 2.1. Electrostatic free energy

Consider a binary electrolyte gel, consisting of a mixture of solid, fluid and ionic phases, immersed in an ionic bath at constant electrostatic potential hydrostatic pressure *p*_{0,} and ionic concentration The stroma can be modelled as such a binary electrolyte of Na^{+} and Cl^{−}. This is justified by the fact that the molar concentration of these ions accounts for about 85% of the total molarity of the aqueous humour [18,19], which acts as an ionic bath. Noting that all ionic species in the stroma have the same valence value of −1, we set the ionic concentration of Na^{+} and Cl^{−} in the aqueous humour to and thereby account for the total molarity of all solutes in the aqueous humour, which is approximately 300 mM [18,19,23].

The reference state for the electrostatic free energy is taken to be the normo-hydrated cornea (i.e. the physiological state) which can be routinely measured. A unit cell reference configuration is denoted *Ω*_{0}, with volume *V*_{0,} and with a fixed charge concentration distribution modelled by the function with The motion maps *Ω*_{0} into the current configuration *Ω* with deformation gradient **F**, considered uniform over the unit cell by virtue of the small unit cell size. The volume of *Ω* is *V* = *JV*_{0} with The fixed charge concentration distribution in the current configuration is modelled by the function *C*_{f}(**x**, *J*) with The fixed charge concentration functions and *C*_{f}(**x**, *J*) are developed in §3.1.

The reference configuration *Ω*_{0} will generally not be a stress-free configuration for the solid (collagen) phase^{1} and it is necessary to introduce a stress-free reference configuration for the solid phase, denoted [3]. The deformation gradient **F**′ from to the current configuration *Ω* may be multiplicatively decomposed such that where **F**_{0} is the deformation gradient from to *Ω*_{0}. Algorithms for finding and **F**_{0} have been described, for example, by Pinsky *et al*. [3] and Gee *et al*. [38].

The Gibbs free energy density, measured per unit reference volume, of the electrolyte gel in *Ω* may be additively decomposed [24] into four components
2.1where is the electrostatic potential and and are right Cauchy–Green deformation tensors. The elastic strain energy density of the solid phase is associated with fibril stretching and matrix shearing and must be based on the three-dimensional collagen organization within the stroma. A detailed description of this term is deferred to §5.

The mean-field approximation for the Helmholtz free energy density of a binary electrolyte in *Ω* is given by [25]
2.2

where *C*_{f}(**x**, *J*) is the fixed charge concentration in *Ω* and measured in moles per litre, *z*_{f} = −1 is the fixed charge valence value, and *R*, *T*, *F* and are the gas constant, temperature, Faraday constant and dielectric permittivity of the electrolyte solvent, respectively. The first term in the integrand measures the free energy of the fixed charge. The second term includes the excess mean concentration of the mobile ions at any point in the electrolyte compared to the bath and, after multiplication by *RT*, may be interpreted as the osmotic work of introducing excess ions into the neighbourhood of the fixed charges. The last term where is the electric field and is the electric displacement, measures the dielectric free energy. Observe that is defined to be the *average* electrostatic free energy density over the unit cell. The average free energy density of the bath is
2.3

Finally, is the potential of the external loads, including body forces and boundary tractions. The list of terms appearing in the free energy expression (2.1) is not exhaustive; for example, the free energies of mixing and disassociation [39,40] can also be considered. However, osmotic energy is dominant in highly hydrated gels [2,39] and additional terms are neglected in the current theory.

### 2.2. Thermodynamic equilibrium

The condition corresponds to electrostatic equilibrium and the Euler–Lagrange equation is the Poisson–Boltzmann equation for the electrostatic potential on *Ω*,
2.4with ionic concentrations satisfying the Boltzmann distribution
2.5

The condition for mechanical equilibrium is found by setting The variation of the elastic free energy is taken with respect to **C**′ to incorporate the solid phase prestress, as discussed above. Noting that and depend on **C** only through *J*, it follows that:
2.6The Euler–Lagrange equation of this variational form on *Ω* is found to be
2.7where the effective Cauchy stress ** σ** is given in the standard manner by
2.8

in which is the second-order push-forward operator from to *Ω*, with The unit cell osmotic pressure *P*_{os} in *Ω* is given by
2.9

It is also seen that,
2.10so that the electrolyte fluid pressure *p* is given by
2.11

### 2.3. Recovery of Donnan osmotic pressure

To illustrate the application of functional (2.2), consider the special case in which the reference fixed charge density is uniform over the unit cell. The fixed charge density in the deformed configuration is then simply 2.12

For this problem, the (Donnan) electrostatic potential, denoted is uniform over the unit cell and be may be solved from (2.4) as 2.13

Observing that and introducing (2.13) into (2.2) and with use of (2.12), the Helmholtz electrostatic free energy density reduces to 2.14

Employing the above equation in (2.9) yields 2.15

which recovers the classical Donnan osmotic pressure [41,42].

## 3. Application to the human cornea

### 3.1. Fixed charge concentration

The predominant corneal GAGs are keratan sulfate, dermatan sulfate and chondroitin sulfate [9,10]. A detailed FACE^{2} [43] analysis reveals that only a portion of the available GAG monosaccharides is sulfated in the human cornea. It has been estimated [23] that the average ionization fraction (i.e. ratio of the charge per disaccharide unit to the number of groups that are ionizable) over all corneal GAGs is *f*_{avg} = 72.6%. If the unit cell within the normo-hydrated stroma contains *N*_{g} GAG chains with average contour length *L*_{c} and if *b* is the length of a monosaccharide unit, the total GAG fixed charge in the unit cell is
3.1

where *e* is the unit charge. The average fixed charge concentration, measured in mM units, is then
3.2

where is the average GAG number density and *F* is the Faraday constant. *C*_{avg} has been measured indirectly [1,2], and the value used in this study is taken from [23] and given in table 1.

X-ray scattering studies under varying tissue hydration [45] suggest that some portion of stromal PGs form a charge-rich and water-binding PG-coating surrounding each collagen fibril. The radius of the coating has been measured to be *r*_{c} = 18.25 nm and this radius is insensitive to hydration over a wide range. The existence of such a surface ultrastructure on the collagen fibrils is corroborated by image studies from [46,47]. In addition, a theoretical study on transparency by Twersky [48] proposed that collagen fibrils must be centred in a transparent coating and the coated fibrils occupy approximately 60% of the matrix volume, giving *r*_{c} = 21.56 nm. Recent three-dimensional electron microscopy reconstructions of corneal collagen and GAG chains [9] also suggests that some GAG chains are in close association with the fibrils while the remaining GAG chains have random orientation in the interfibrillar fluid.

Given the regular arrangement of fibrils within a lamella, the nanoscale unit cell is selected to be of the form shown in figure 1*c*. It is a rectangular prism with cross-section as shown in figure 1*c* and with axial axis aligned with the collagen fibrils. Each fibril has a GAG-coating region. As described below, fixed charge concentrations in the unit cell will depend on various volume fractions, including the collagen and coating regions.

Let the unit cell current configuration *Ω* be partitioned into collagen fibril, coating and inter-coating domains (figure 1*c*) such that and where We seek expressions for the fixed charge density over the subregions of *Ω* as a function of dilation *J*. A parameter *λ* is used to partition the total unit cell GAG fixed charge *Q* into interfibrillar GAG fixed charge *λQ* and coating GAG fixed charge (1 − *λ*)*Q*. The value of *λ* = 0.65 has been estimated from swelling pressure measurements [23]. Note that because all GAG chains bind at one end to a collagen fibril, the interfibrillar GAG chains must pass through the GAG-coating region and their charges therefore combine in that region. The coating and interfibrillar GAGs are assumed to produce uniform charge concentration distributions over their respective domains.

In determining the fixed charge concentrations as a function of macroscopic dilation *J*, account must be taken of keratocyte volumes which are assumed to exclude GAG fixed charge. Keratocytes are morphologically flattened fibroblasts that occur between lamellae [11]. They are considerably larger than the unit cell and their effect is therefore introduced in an average way (figure 1*b*). Stromal collagen and keratocyte volume fractions are denoted *ϕ*_{col} and *ϕ*_{ker}, respectively. Assuming that keratocyte cells dilate with the tissue dilation [49] and assuming collagen fibrils are non-swelling [1,50,51], the charge concentrations in the current unit cell configuration *Ω* are then given by *C*_{coat} and *C*_{inter} in the coating and inter-coating regions, respectively, as
3.3and
3.4where *ϕ*_{coat} is the volume fraction of the coating region. Note that the coating volume is independent of dilation *J* as required by the experimental observations of Fratzl & Daxer [45]. In (3.3) and (3.4), *C*_{avg} = *Q*/*FV*_{0} is the average charge concentration in the reference configuration *Ω*_{0} given by (3.2). It is easily verified that (3.3) and (3.4) conserve unit cell total fixed charge under any volume dilation *J*. Volume fraction values used in this study are provided in table 1.

Finally, assuming that keratocytes have the same mass density as water, the unit cell volume dilation *J* can be related to the stromal hydration *H*_{w}, defined as water weight per unit dry (collagen) weight, as follows:
3.5where *ρ*_{w} = 1 g cm^{−3} and *ρ*_{col} = 1.36 g cm^{−3} are the mass density of the water and dry collagen fibrils [51], respectively, and *ϕ*_{r} is the relative volume fraction of collagen molecules within the fibrils (which excludes the intrafibrillar bound water) [1,50]. The value of *ϕ*_{r} is estimated to be 0.75 such that the stromal hydration *H*_{w} is 3.2 at the normal condition *J* = 1 [1]. Therefore, equation (3.5) is rewritten as
3.6

The linearity is consistent with measurements [52,53]. It should also be noted that the estimated value *f*_{r} = 0.75 is in very close agreement with the value 0.77 reported by Goodfellow *et al*. [50].

### 3.2. Active endothelial ion transport

The corneal endothelium is a 5 µm thin cellular monolayer located on the posterior surface of the cornea. It is permeable to water, metabolic species, including glucose and lactate ion, and other salt ions [11]. The endothelium regulates stromal hydration by providing active ion transport across the endothelium, from the stroma to the aqueous humour [1,12]. To allow sufficient generality, we model steady and independent anion and cation active transport and measure their effects (singularly or combined) on osmotic pressure by deriving a modified Boltzmann ion distribution for the stroma and the corresponding electrostatic free energy density for the *in vivo* cornea.

The endothelial layer is here modelled as an ideal membrane in which the action of the ion pumps is represented by independent active steady anion flux *J*_{a−} and cation flux *J*_{a+}, transporting ions out of the stroma and into the aqueous humour [21]. Because the endothelium is thin compared with the stroma, we will assume transport across the endothelium is one-dimensional. A coordinate *z* is introduced with *z* = *z*_{0} at the endothelium–aqueous interface and *z* = *z** at the endothelium–stroma interface. At equilibrium, the net ionic fluxes *J*_{−} and *J*_{+}, which result from both active and passive transport, must vanish so that
3.7where *L*_{i} is the membrane permeability of ionic species *i* and *J*_{ai} > 0 is the steady active ionic flux across the endothelium, from the stroma and into the aqueous humour. The ionic electrochemical potential *μ*_{i} is given by
3.8where *M*_{i}(*i* = +,−) are the cation and anion atomic weights, respectively, are the ionic chemical potentials at the reference state. Using (3.8) in (3.7) gives
3.9and adding both equations in (3.9) results in
3.10

Integrating (3.10) across the endothelium leads to
3.11where is the thickness of the endothelial layer, are the ionic concentrations at and *C*_{0} is the bath (i.e. aqueous humour) ionic concentration. Observe that the right-hand side of (3.11) is dimensionless as required. Defining,
3.12

it follows that and from (3.11), 3.13

Finally, integrating (3.9) across the endothelium, and taking the electrostatic potential in the aqueous gives
3.14where is the electrostatic potential at *z* = *z**.

At thermodynamic equilibrium, the gradients of the ion electrochemical potentials vanish and where *μ*_{i} is the electrochemical potential at any point in the stroma, and assuming that ionic concentrations are continuous across the endothelial–stroma interface, this condition implies that
3.15

Introducing (3.14) into the above equation allows the stromal ionic concentrations *C*_{i} to be expressed as
3.16

which generalizes the Boltzmann distribution (2.5). From this result, it may be seen that at any point in the stroma, the ionic concentrations *C*_{+} and *C*_{−} must satisfy the condition
3.17

### 3.3. Modified electrostatic free energy

Based on (3.16), a modified Poisson–Boltzmann equation for the electrostatic potential in the stroma is introduced as 3.18

and the corresponding modified electrostatic free energy density is 3.19

It may be verified that this functional is stationary at the solution of the modified Poisson–Boltzmann equation (3.18), i.e. It may also be observed that if the active fluxes *J*_{a−} and *J*_{a+} increase, the osmotic free energy, measured by the second term in the integrand, reduces.

### 3.4. Analytical approximation for osmotic pressure

Recognizing that the fixed charge density *C*_{f}(**x**, *J*) is piecewise constant over the unit cell, and noting that the unit cell is large compared to the Debye length [23], we employ a piecewise constant Donnan potential as an approximation for the exact electrostatic potential *φ* which is required to satisfy the Poisson–Boltzmann equation (3.18). By assumption and equation (3.18) may be used to solve for the Donnan potential. Recalling that the fixed charge density for and for the Donnan potentials over *Ω*_{coat} and *Ω*_{inter} are given by
3.20where the non-dimensional fixed charge concentrations and are defined by
3.21and where the fixed charge densities *C*_{coat}(*J*) and *C*_{inter}(*J*) are given by (3.3) and (3.4), respectively.

Employing potentials (3.20), the Helmholtz free energy density (3.19) reduces to
3.22where and It can be shown that is convex in *J*. The unit cell osmotic pressure then follows from (2.9) as giving
3.23

in which a modified inverse hyperbolic sine function is introduced to simplify the expression 3.24

The osmotic compressibility, defined as is the core function in the stroma electrolyte tensor required for the finite-element tangent operator (see the electronic supplementary material for details) and is given by
3.25Equations (3.23)–(3.25) show that the osmotic pressure and compressibility are dictated by the combined ionic transport parameters *Q*_{+}*Q*_{−}. These parameters must be estimated indirectly since direct measurement is not currently feasible. This is considered in the following section. In the *ex vivo* cornea, no active transport occurs and *Q*_{+} = *Q*_{−} = 1. This condition is employed in §6.1 to model *in vitro* free swelling and swelling pressure experiments on isolated stroma samples in ionic baths to provide evidence for the accuracy of the proposed model.

## 4. Model calibration from imbibition pressure

The presented model requires the determination of *Q*_{+}*Q*_{−}. From (3.12), it is found that
4.1the above equation implies that any increase in permeability, described by *L*_{i}, that is not balanced with a compensatory increase in ionic pump function, described by *J*_{ai}, will result in stromal swelling and eventually clinical oedema [54]. The permeability and pumping rates have been estimated via transient swelling experiments [55,56], but these are not relevant for equilibrium conditions. Instead, we shall infer the value of *Q*_{+}*Q*_{−} from measurement of the imbibition pressure in rabbit.

Imbibition pressure is measured in the living cornea by inserting a saline-filled cannula into the stroma [57]. At equilibrium, a stable suction pressure is developed in the cannula and is referred to as the imbibition pressure *P*_{imb}. As the saline solution is in contact with the stroma in front of the cannula tip, Donnan equilibrium is reestablished locally [21] and equilibrium of water requires
4.2where and are the water electrochemical potentials of the stroma in front of the tip and of the saline solution in the cannula, respectively. The stromal fluid pressure is *p* and *P*_{imb} is the hydrostatic (imbibition) pressure measured in the cannula. Setting the local osmotic pressure difference equation (4.2) implies
4.3

As noted by Bryant & McDonnell [21], *П* also corresponds to the *ex vivo* osmotic pressure because the canula and ionic bath have the same ionic concentration *C*_{0}. In this case, equation (4.3) states that imbibition pressure is the difference between the stromal fluid pressure and the *ex vivo* osmotic pressure. However, direct experimental *in vivo* measurement of imbibition pressure in rabbit [57] indicated that the following relationship holds:
4.4where is the IOP and is the independently measured swelling pressure in an ionic bath of concentration . As discussed in §1, for the stroma it can be safely assumed that swelling and osmotic pressure coincide and therefore equations (4.3) and (4.4) allow the conclusion that in rabbit the stromal fluid pressure is equal to the IOP, i.e. *p* = IOP at normal hydration. Furthermore, the measured *P*_{imb} was found to have insignificant variation throughout the cornea [57]. It is therefore reasonable to conclude that the fluid pressure is *p* = IOP everywhere in the normal *in vivo* rabbit cornea. Lacking any alternative, we shall take this condition to also hold in the human cornea.^{3}

Recalling that equilibrium fluid pressure in the cornea is given by (2.11), and setting *p*_{0} = IOP, it follows that:
4.5

However, by the above argument based on imbibition pressure measurements, we must set *p* = IOP. In this case, equation (4.5) implies that in the normal *in vivo* cornea,

The condition is used to estimate *Q*_{+}*Q*_{−} as follows. Figure 2 depicts the variation of *P*_{os} with *Q*_{+}*Q*_{−} at normal hydration, determined from (3.23). A nearly linear relationship between *P*_{os} and *Q*_{+}*Q*_{−} is found. The curve indicates that stromal osmotic pressure reduces monotonically with increasing active ionic flux values. A calibrated value of *Q*_{+}*Q*_{−} = 0.965 is found for the normo-hydrated *in vivo* cornea. As will be discussed in §6.2.1, the positive slope of figure 2 is a necessary requirement for hydration regulation.

## 5. Collagen organization and stromal elasticity

### 5.1. Angular averaging

In this section, we review the modelling of the stromal fibres which contribute to the elastic strain energy density appearing in (2.1). It is assumed that is additively decomposed into an anisotropic part describing the strain energy density of the three-dimensional collagen fibre network and an isotropic part that gives a simple phenomenological description of the small background shear stiffness of the extracellular matrix such that 5.1

Recall from §2.1, that describes the deformation measured from the stress-free configuration , and **F**′ is the deformation gradient from the stress-free configuration to the current configuration *Ω*. The invariant is used to describe an isotropic contribution based on the unimodular right Cauchy–Green deformation tensor and [58].

The invariant which measures the stretch along a fibre with direction **A**, is used to describe elastic anisotropy. It is important to note that depends on the dilation through its dependence on **C**′, allowing the description of fibre stretching resulting from electrolyte swelling. The form of is based on the use of angular integration [3,59] such that
5.2
5.3where ** ω** is the unit sphere,

*ρ*(

**X**,

**A**) is the distribution of fibre directions, is the normalization and the fibre stretching energy is described below.

The function *ρ*(**X**, **A**) may be based on fibre orientation information obtained from X-ray diffraction and SHG imaging as described in [6]. A fibre direction **A** can be expressed in spherical coordinates with azimuthal angle and inclination angle As the anterior and posterior surfaces have different curvatures, a local coordinate system is introduced that varies pointwise through the corneal thickness and matches both the anterior and posterior curvatures at those limits; details can be found in [6]. Every material point within the cornea has a unique angular distribution with no perfect symmetry over quadrants. We set where is a point on the cornea anterior surface, is the non-dimensional corneal depth measured from the anterior surface, and assume,
5.4where is based on X-ray diffraction data and is based on SHG imaging. The current model generalizes the approach in [6] in two respects: incorporates depth-dependence of preferred fibre directions, and incorporates direct statistical measurement of inclination data. These are briefly summarized below.

### 5.2. Azimuthal distribution

In [3], analytical distribution functions were introduced to represent data from X-ray diffraction experiments by Meek *et al.* [60]. Following [6], we have eliminated the need to introduce analytical distribution functions and directly employ the raw X-ray data. These data provide scattering intensities versus orientation on a discrete grid of points over the cornea's anterior surface. At any point in the cornea, the orientation distribution is obtained by interpolating the X-ray diffraction measurements taken at the nearest four grid points using bilinear Lagrange functions. This procedure ensures that any subtle variations in distributions found in the X-ray data are faithfully reproduced in the model. At any scan point the X-ray scattering data which corresponds to a measurement through the entire thickness, can be analysed in the following way. The total scatter *η*_{total} is the area under the curve,
5.5This can be additively decomposed into aligned and uniform parts where and (figure 3*a*).

The percentage of all aligned fibres is then 5.6

The depth-dependence of aligned fibres may be introduced based on the work of Abahussin *et al*. [33], who performed X-ray diffraction studies on isolated third-thickness stromal samples and characterized the variation of scattering intensity with depth at the corneal centre. They found the average percentage of total fibres exhibiting preferred azimuthal directions *v*_{align} was 22%, 31% and 42% in the anterior, central and posterior thirds, respectively. However, depth-dependent alignment was not measured at the cornea's periphery where there is significantly more alignment. A linear best-fit of alignment as a function of depth for the data is 0.3*s*+1/6 with *R*^{2} = 0.997. Based on the slope of this function, a first approximation for the depth-dependence of aligned fibres is For the central full-thickness X-ray scans in the current database of implemented X-ray diffraction data (four donor corneas, including a right and left cornea from the same donor), we find for all subjects, which corresponds well to the mid-thickness value measured by Abahussin *et al*. [33]. Thus, to evaluate *δ* we set which gives Then,
5.7defines the fraction of aligned fibres based on in-plane position and depth *s* (figure 4).

In order to generate a depth-dependent distribution, the full thickness distribution *η*_{total} is augmented by Δ(*s*) such that (figure 3*b*). Solving for Δ(*s*) gives
5.8

The final distribution will then be,

Incorporating the X-ray diffraction data in this manner describes the depth-dependent anisotropy resulting from the well-documented S–I and N–T preferred directions of lamellae in the vicinity of the corneal vertex as well as the circumferentially preferred orientations at the limbus [44].

### 5.3. Inclination distribution

The presence of inclined fibres has been well-documented with macroscale SHG imaging [34–36]. However, because the angular distribution of inclined fibres was missing, an approximation was proposed by Petsche & Pinsky [6] based on inspection of the macroscale images and calibration based on the depth-dependence of shear properties. Recently, Winkler *et al*. [37] have processed SHG images to detect and quantify the spatial distribution of inclined lamellae and their distribution is employed in this study. Winkler *et al*. [37] analysed half cross-sectional image stacks to measure the inclination angle of every lamella in the anterior half of the stroma. Lamellae were binned based on depth, corneal quadrant and radial position. For each bin, a histogram of fibre inclination was fit to a Gaussian distribution for the full-width at half-maximum (FWHM), with the mean fixed at 0° from the plane tangent to the anterior surface (figure 5).

They found no significant differences between quadrants and radial position but measured a significant variation with depth. On the anterior half of the cornea, where analysis was performed, a linear decrease of FWHM was found. In order to complete the data through full thickness, it is assumed that FWHM extrapolates linearly to zero at the posterior surface. A best-fit line that forces no inclined lamellae at *s* = 1 (*R*^{2} = 0.87) was used to find FWHM as a function of depth (figure 6).

Using the one-to-one relationship between FWHM and the standard deviation *σ* this gives
5.9

The inclination distribution therefore takes the form of a Gaussian 5.10

where the mean is *π*/2 because *ϕ* is the zenith angle [44].

### 5.4. Fibril and shear elasticity

The term in (5.3) represents the strain energy density resulting from fibril stretching. In tension, collagenous tissues are characterized by a J-shaped stress–strain curve resulting from unkinking of collagen molecules at low strains followed by the stretching of stiffer covalent bonds at high strains. However, they exhibit negligible compressive stiffness due to buckling effects. To capture this behaviour, Markert *et al*. [61] introduced a strain energy function of the form
5.11

This model has been applied to the cornea by Studer *et al*. [5]. The model parameters *α* and *β* have been identified using corneal inflation experimental data and optimal values found to be *α* = 30 kPa and *β* = 75 [44]. These values satisfy the constraints on *α* and *β* given by Markert *et al*. [61] and render convex in

The matrix shear term is modelled as simply as possible and is assumed to be
5.12where *μ* = 5 kPa is the matrix shear modulus which has been measured experimentally by Petsche *et al*. [62] and Sondergaard *et al*. [63].

Finally, the elasticity tensor and electrolyte tensor needed for the finite-element tangent operator, are described in the electronic supplementary material.

### 5.5. Boundary conditions

The *in vivo* cornea is loaded by the IOP at the endothelial–aqueous interface (figure 1*a*). The anterior stroma is sealed by the epithelium, which is a stratified squamous cellular membrane that is composed of four to six cellular layers and approximately 50 µm in thickness. In contrast to the endothelium, epithelial cells exhibit tight junctions and the epithelium is essentially impermeable to water and ions under normal conditions [11]. As the mechanical stiffness of both epithelium and endothelium is insignificant compared to the stroma [11,64], the mechanical roles of these bounding layers are reasonably ignored. Therefore, for modelling purposes, the IOP is assumed to act directly on the stroma and the posterior boundary condition is written as
5.13where ** σ** is the effective Cauchy stress,

*p*is the stromal fluid pressure,

*p*

_{0}= IOP is the hydrostatic pressure in the aqueous humour and

**n**is the normal unit vector. At the epithelium, no external pressure is applied and then 5.14

A full interpretation of the above boundary conditions is extremely important since it provides insight into specific corneal structural function, and in particular, the implied necessity of fibre inclination in the anterior (such that needed to balance the stromal fluid pressure at the anterior corneal surface. A detailed discussion is presented in §7.

## 6. Results

### 6.1. Swelling in the *ex vivo* stroma

#### 6.1.1. Confined and unconfined swelling pressure

Finite-element predictions for *ex vivo* swelling pressure were compared to experimental measurements on a sample of isolated corneal stroma immersed in a bath solution (figure 7*a*). Stromal swelling pressure is the mechanical pressure required to maintain a specified thickness of the tissue at a specified bath ionic concentration. *Ex vivo* conditions (no active ion transport) are modelled by setting the ionic transport function *Q*_{+}*Q*_{−} = 1 in equation (3.23). The stroma sample was modelled as a 7 mm diameter cylinder with an initial thickness of 0.5 mm and immersed in an ionic bath of concentration *C*_{0} = 150 mM [20,52]. The finite-element mesh employed 2500, 27-node (tri-quadratic) hexahedral elements. Fixed displacement boundary conditions were applied on the bottom surface and prescribed transverse displacements were applied to the top surface to model changes in sample thickness. The equilibrium swelling pressure at any sample thickness *t* is computed as the total reaction force on the top surface divided by the surface area. In order to model confined [20] and unconfined [52] experimental measurements, edge boundary conditions were taken as zero radial displacement or free, respectively (figure 7*a*). Predictions are compared with experimental measurements [20,52] in figure 7*b*. The results show good agreement for the full range of experimental measurements with *t* varying from 0.75 mm (swollen) to 0.2 mm (highly compressed), with the corresponding swelling pressure values varying over two orders of magnitude. Swelling pressure predictions for the confined and unconfined cases are nearly identical; this is expected because lateral expansion of the sample is restricted by the collagen fibrils and, as a consequence, stromal volume change in these tests is predicted to be associated with changes in sample thickness [1].

#### 6.1.2. Free swelling of a stromal sample

The free swelling of the above described cylindrical sample of *ex vivo* stroma is modelled to investigate the role of the collagen organization in mediating the swelling response. Results for edge-confined free swelling of the stroma sample are shown in figure 8*a* (solid curve), which depicts the ratio of swollen thickness to original thickness (swelling ratio) versus bath ionic concentration *C*_{0}. For dilute bath solutions, the model predicts that the stroma sample will swell to approximately four times its original thickness, which may be compared to experimental measurements in de-ionized water [28,29] in which human corneas were observed to swell to approximately three times their original thickness. For concentrated bath solutions, the abundance of ions results in ionic shielding of the fixed charges, reduction in osmotic pressure and minimal swelling. Both limiting states are captured by the swelling predictions.

It has been observed in free swelling experiments that the anterior third of the stroma remains virtually unswollen with most swelling taking place in the deeper stroma [28,29]. We used the model to investigate the depth-dependence of stromal swelling at two bath ionic concentrations *C*_{0} = 150 mM and 150 × 10^{−6} mM, corresponding to the physiological state and de-ionized water, respectively. The results are presented in figure 8*b* which depicts the profile of local volume dilation *J* across the depth of the stroma (*J* = 1 corresponding to the normo-hydrated state). At *C*_{0} = 150 mM, the local volume dilation varies from 1.5 in the anterior to 1.8 in the posterior. Massive swelling occurs when *C*_{0} = 150 × 10^{−6} mM and volume dilation varies from 1.8 in the anterior to 10.0 in the posterior. The highly non-uniform swelling across the corneal thickness agrees qualitatively with experimental observations [28,29]. Comparison of the two results in figure 8*b* indicates that the stromal anterior third thickness is predicted to maintain its thickness for any bath concentration, indicating the rigidity of this region with respect to extreme hydration changes.

A theoretical study was also undertaken to confirm the importance of lamellae inclination with respect to anterior stromal rigidity. The above study was repeated while constraining lamellae to have zero inclination. This is accomplished by replacing the Gaussian inclination distribution given by (5.10) with where *δ* is the Dirac delta. In this case, all lamellae throughout the corneal depth will have zero inclination. As shown in figure 8*a* (dashed curve), for dilute bath solutions the stroma sample now swells to almost 10 times its original thickness. In fact, the only mechanical constraint preventing infinite swelling is the work done by the matrix shear stiffness which bounds the swelling. As expected, the model still predicts little swelling when the bath ionic concentration is high.

### 6.2. Swelling in the *in vivo* cornea

#### 6.2.1. Effect of active endothelial ion transport

Recall that the active endothelial ionic transport term (equation (4.1)) satisfies with *Q*_{+}*Q*_{−} = 1 corresponding to no active transport and *Q*_{+}*Q*_{−} = 0.965 corresponding to normal active transport in the *in vivo* cornea. Predicted values of osmotic pressure *P*_{os} (equation (3.23)) are plotted in figure 9*a* for these two values of the ionic transport function and also for *Q*_{+}*Q*_{−} = 0.930 which corresponds to hyper-active ionic transport. Osmotic pressure is reduced with reducing values of *Q*_{+}*Q*_{−}. For *Q*_{+}*Q*_{−} = 1, the osmotic pressure is positive for all dilation *J* and the tissue will tend to swell. For the normal *in vivo* cornea with *Q*_{+}*Q*_{−} = 0.965, it is seen that a reduction in dilation from the physiological condition (*J* = 1) will result in positive stromal osmotic pressure which will produce a tendency to swell. On the other hand, an increase in dilation will result in negative stromal osmotic pressure and produce a tendency to de-swell. Consider the curve for the hyper-active *Q*_{+}*Q*_{−} = 0.930. In this case, positive or negative deviations in dilation from *J* = 1 will always result in negative osmotic pressure and produce a tendency to de-swell. These conditions illustrate the hydration regulation mechanism of active endothelial ion transport. Increasing active transport shifts the osmotic pressure–dilation curve downwards in figure 9*a*. Predicted values of osmotic compressibility *K*_{os} (equation (3.25)) versus dilation *J* are shown in figure 9*b* for no active transport. Values closely match measurements reported in [52].

#### 6.2.2. Swelling of a cornea with Fuchs' dystrophy

Fuchs' dystrophy is usually characterized by morphological changes in endothelial cells or by an accelerated loss of endothelial cells [11,54]. In this situation, the cornea will swell due to increasing endothelial permeability, decreasing active ion flux, or both mechanisms simultaneously. The effect of these pathological changes on stromal ionic concentrations is described by the ionic transport function *Q*_{+}*Q*_{−} (equation (4.1)). Because the absolute values of endothelial ionic permeability *L*_{+} and *L*_{−} and active ionic fluxes *J*_{a+} and *J*_{a−} have limited clinical significance, swelling effects due to relative changes in these parameters are considered. For simplicity and without loss of generality, we consider active anion transport only and set *Q*_{+} = 1. Let and be the anion membrane permeability and active anion flux rate of the normal cornea, respectively, and consider pathological increases in ionic permeability such that and reductions active flux such that For any value of and the value of can be found through equation (4.1)
6.1where is the calibrated value for normal cornea as determined in §4.

The cornea was modelled geometrically with a central thickness of 520 µm and the anterior and posterior surfaces defined as spherical surfaces with radii of 7.87 and 6.7 mm, respectively. Part of the sclera was also modelled with anterior and posterior surfaces defined with radii of 12.00 and 11.38 mm, respectively. Because the sclera swells much less than the cornea [8], its bulk behaviour was modelled as a (non-swelling) compressible neo-Hookean material, as described by Petsche & Pinsky [6], and the limbal corneal collagen fibre elasticity was extended to the scleral tissue. At the junction of the cornea and sclera, there is an abrupt interface of swelling and non-swelling tissue, which should be replaced by use of a gradual transition (although we have not done so). Homogeneous Dirichlet boundary conditions are applied to fix the periphery of the scleral section. The solution proceeds by first solving for the normal cornea^{4} with and then solving a sequence of problems in which is systematically varied to simulate changes in endothelial permeability and active anionic flux.

Figure 10*a* shows the predicted swollen central corneal thickness (CCT) for endothelial ionic permeability increased with and with normal active ionic flux Figure 10*b* shows the predicted swollen CCT for active ionic flux rate reduced up to a factor of 0.2, so that with normal anionic permeability In both cases, the cornea swells. The predicted maximum swollen CCT of approximately 690 µm lies in the range of clinical observation (500–1000 m) measured for patients with Fuchs' dystrophy [30].

The model predicts that at any value of the swollen CCT, the anterior surface deforms much less than the posterior surface, with swelling concentrated in the posterior region. This prediction agrees with the clinical observation of Brunette *et al.* [30] that the anterior surface of the cornea is nearly normal among patients with Fuchs' dystrophy, whereas the posterior surface shows significant change. Figure 11 provides a fringe plot of axial (vertical) displacements due to swelling resulting from a 78% reduction in active ion transport, The results again suggest the stability of the anterior surface with respect to swelling resulting from the presence of inclined lamellae.

#### 6.2.3. Swelling due to changes in intraocular pressure

An example of corneal oedema with an intact and functional endothelium occurs in acute glaucoma. In this case, elevated IOP combined with normal stromal osmotic pressure can create an increase in corneal thickness [65]. Likewise, corneal thickening is reported in patients with ocular hypertension [66,67]. We emphasize that we are considering the equilibrium solution that is achieved after long-time exposure to increased IOP. It is very likely that the short-time non-equilibrium response of the cornea is that of elastic thinning. In fact, if we model an increase in IOP and do not allow the fluid pressure to equilibrate, the model does predict small corneal thinning. However, persistent elevated IOP produces an increased pressure gradient which drives more fluid across the endothelium, creating oedema of the stroma (and epithelium), as seen in acute angle-closure glaucoma [65,68]. This is undoubtedly a simplified view of a complex problem but it is supported by predictions from the current model.

There are two competing mechanisms with regard to corneal swelling. The increased (mechanical) IOP loading at the endothelial–aqueous humour interface acts to thin the cornea, but an increase in stromal fluid pressure produces a countering ‘inflation’ effect. At Bowman's layer, the stromal fluid pressure must be balanced by forces in the inclined lamellae, because there is no external pressure loading at that boundary which can balance the stromal fluid pressure (figure 12). As shown in figure 13 (solid curve), at 44.5 mmHg, the CCT is predicted to swell by approximately 15 µm, which is in agreement with clinical data [65] in which the CCT was measured to be 575 and 593 µm for a group of patients who had normal IOP (average 15.8 mmHg) in one eye and high IOP (average 44.5 mmHg) in the other eye, giving a thickness difference of 18 µm. In the simulation of increasing IOP, the active ionic pumping was maintained at a normal level throughout. It is important to note that other measurements, performed on eyes with retinal detachment or non-human species, or performed within a relatively short timescale, have concluded an opposite trend [69–71].

To illustrate the importance of lamella inclination in resisting the IOP, the cornea was reanalysed with lamella inclination constrained to be zero, as described in §6.1.2. In this case, the response shown in figure 13 (dashed curve) includes a significant further increase in swelling.

#### 6.2.4. Biomechanical behaviour after LASIK surgery

LASIK (laser *in-situ* keratomileusis) is commonly used to correct refractive errors including myopia and hyperopia. An anterior circular hinged flap of about 150 µm thickness is first created using a femtosecond laser. The flap is then folded back to expose the stromal bed which is reshaped by excimer laser ablation. The procedure is completed by returning the flap to conform with the ablated profile. Although laser ablation creates a new central corneal curvature with the goal of achieving emmetropia, the biomechanical response of the residual (significantly thinned) stromal tissue is an important but challenging factor that should be considered [72]. In this study, we assume that active endothelium ionic pumping is functioning normally and use the proposed model to investigate the biomechanical response of the cornea, including stromal fluid pressure and local swelling, resulting from a typical and representative myopic correction.

The geometry shown in figure 14 is for the pre-operative cornea and defines the shape of the normo-hydrated cornea and surgical profile while the cornea is supporting an IOP of 15 mmHg. The ablation depth profile was determined using theory presented by Munnerlyn [73] for an ablation radius of curvature of 8.847 mm (corresponding to an intended correction of approx. −5 dioptres). The CCT is 520 µm, and the maximum ablation depth of 87 µm combines with the 150 µm flap to produce a residual central stromal thickness of 283 µm. With this ablation profile, the residual stromal thickness is 54.4% at the cornea centre and gradually increases to 71.2% at the periphery of the flap. The flap is returned after ablation (where negative osmotic pressure holds it in place by suction) but it is no longer mechanically integrated with the stroma and can be easily separated at any time after surgery [72]. LASIK therefore effectively removes that part of the stroma in which lamella inclination is most pronounced (the anterior third).

Because the collagen fibril distribution is not symmetric about any axis, a full 360° model is required. The eye is supported in the eye socket by fatty tissue and stabilized by six muscles responsible for eye movement. However, it is not necessary to extend the sclera to the posterior pole; numerical experiments were conducted and it was determined that extending the scleral tissue for an arc length of 2 mm was sufficient, at which point the tissue can be ‘clamped’. The finite-element model consists of 7008 27-node hexahedral elements. The solution was obtained in two phases. In the presurgical phase, the normal cornea was loaded by an IOP of 15 mmHg and a global fixed-point algorithm was employed to obtain the collagen prestress. With this prestress, the normal cornea is in equilibrium under the IOP and is in the pre-operative geometry of figure 14. In the postsurgical phase, the elements corresponding to the flap and ablation region are removed from the mesh while maintaining the IOP and normal ionic pumping. The model is no longer in equilibrium and the new deformed state is obtained by Newton's method.

We first compare the predicted biomechanical response to clinical measurements. Ortiz *et al.* [74] made measurements of changes in corneal anterior curvature and correlated these to the surgical radius of curvature in the ablation profile. The case study analysed 85 eyes with myopia or myopic astigmatism and subdivided them into two groups: 44 with flaps created by microkeratome (‘mechanical’) and 41 with flaps created by femtosecond laser (‘laser’). The biomechanical response to LASIK was evaluated by means of the coefficient of curvature change defined by
where *R*_{post} is the postsurgical corneal radius at one month and *R*_{S} is the surgical radius imposed on the stromal bed by ablation. The postsurgical corneal radius was measured by topography over a central 5 mm diameter central zone. Likewise, the postsurgical corneal radius of the finite-element model used least-squares fitting over a 5 mm diameter zone.^{5} The mean value of *C* found by Ortiz *et al*. [74] is compared to the predicted value found by the present model in table 2. Ortiz *et al.* [74] also performed a statistical analysis for a linear fit function between *R*_{post} and *R*_{S} and found for the ‘mechanical’ group and for the ‘laser’ group. Using the model value the fit functions for *R*_{post} are compared to the predicted value in table 2. Finally, the spherical equivalent measures the power change in dioptres achieved by changing the presurgical radius *R*_{pre} to the postsurgical radius *R*_{post}. Clinical values for SE with *R*_{s} = 8.847 are compared to the predicted value^{6} in table 2. It may be seen that the finite-element model predictions agree closely with the clinical data.

We analysed the fluid pressure and swelling response of the postsurgical cornea. The fluid pressure is given by (equation (2.11)), where *P*_{os} is the osmotic pressure and *p*_{0} is the bath pressure (IOP) of 15 mmHg. Because most inclined lamellae have been removed by surgery, it may be anticipated that the fluid pressure must approach zero in the anterior region of the residual stroma. Since *p*_{0} is constant, this implies that the tissue will swell and the osmotic pressure will become negative (see also figure 9*a*). Fluid pressure fringes are shown in the cutaway half-mesh of figure 15*a* and local volume dilation fringes are shown in 15*b*. Figure 15*a* shows that the fluid pressure reduces to approximately half of its normal value in the ablation region.^{7} The local swelling shown in figure 15*b* is necessary to produce a negative osmotic pressure. It is concentrated in the upper region of the residual stroma, particularly where the ablation is deepest and lamella inclination is lost. It is observed that the swelling is modest, around 3.5% maximum and decreasing with depth. This is so because only small swelling is needed to produce negative osmotic pressure of the order of 15 mmHg (figure 9*a*). This will add some contribution to the biomechanical response in LASIK, but it is probably secondary to the elastic response. However, if IOP increases, then the swelling increases significantly, as shown in table 3.

## 7. Discussion

The proposed model follows the general theory for hydrophilic and highly swollen polyelectrolyte gels introduced by Katchalsky & Michaeli [24] combined with the membrane transport theory of Kedem & Katchalsky [16]. The main work was the characterization of the stromal electrolyte and elastic free energies, including the modelling of modifications due to active endothelial ionic transport. The resulting model predicts the *ex vivo* (no active endothelial ionic transport) osmotic pressure and osmotic compressibility with good accuracy compared to experimental measurements (as described in §6.1.1) and provides some confidence in the general framework. Previous equilibrium modelling of corneal swelling by Bryant & McDonnell [21] was based on numerical solution of the spherically symmetric (one-dimensional) fluid and ionic continuity equations in a multiphasic theory [22,75]. The current energy-based approach is equivalent (see the electronic supplementary material) in that it implies satisfaction of the same governing equations but it offers some distinct advantages.

The energy approach provides an average (or homogenized) energy-equivalent value of osmotic pressure over a unit cell (see equation (2.9)). This is useful for consistently embedding nanoscale effects, for example, the distribution of GAG charge into collagen coating and non-coating regions, in a macroscopic theory. By contrast, the triphasic theory produces a pointwise varying expression for the osmotic pressure. The energy approach is also convenient for describing the deviation of the stromal electrolyte from ideal conditions, including electrolyte volume exclusion effects associated with the collagen fibrils and stromal cell (keratocyte) population. Also noteworthy is the fact that the energy formulation is a natural approach for finite-element approximation which can be based directly on the variational forms provided by the free energy terms. The proposed formulation only requires solving for the displacement field, which gives practical convenience in finite-element implementation. The model contains numerous physical constants and parameters (see reference values in table 1), but most are well characterized and all simulations presented used these reference values—no problem-dependent parameters appear in the theory.

Modelling the influence of active endothelial ionic fluxes on the biophysics of corneal tissue hydration in a fully three-dimensional context presents a challenge. In this approach, we have analytically integrated the phenomenological KK [16] membrane transport equations to describe the membrane osmotic gradient as a function of the active ionic fluxes and endothelial ionic permeability.^{8} This leads to a modified Boltzmann ionic distribution in the stroma and corresponding electrostatic free energy. The effects are contained in the introduced ionic transport parameters *Q*_{+} and *Q*_{−} which correspond to cation and anion transport, respectively. Interestingly, for a binary electrolyte the influence of active ionic transport appears as a function of the combined parameter *Q*_{+}*Q*_{−}. The model suggests that active endothelial ionic transport reduces osmotic pressure and hence swelling pressure (these are essentially equivalent in the case of the highly hydrated cornea) and the pump-leak mechanism can be explained as active pumping producing a downward shift of the swelling pressure–hydration curve as discussed in §6.2.1 and illustrated in figure 9*a*.

Precise quantification of the fluid and osmotic pressures remains an open question for the human *in vivo* cornea, as direct measurements are not currently feasible. Conventional analysis [15,18,19] predicts a negative fluid pressure throughout the cornea. In this study, osmotic calibration based on imbibition pressure measurements for rabbit indicate positive fluid pressure at normal hydration, which is in agreement with [21]. The calibration suggests that stromal fluid pressure must be close to IOP and osmotic pressure must be close to zero. It is interesting to analyse how the forces acting at the corneal bounding surfaces are balanced under these conditions. The posterior stroma boundary condition, equation (5.13), may be rewritten in the direction normal to the corneal posterior surface as where *σ*_{el} is the collagen effective stress component and *p*_{0} denotes the IOP. As the posterior collagen fibres are parallel to the corneal surface, they can provide negligible force in the transverse direction and so Thus, it may be concluded that the *in vivo* osmotic pressure *P*_{os} must be zero at the posterior stroma—in agreement with the calibration.

The situation is quite different at the epithelial layer. Here there is no external pressure loading to balance the stromal fluid pressure (figure 12). Inclined lamellae insert into Bowman's layer [35] where they must act to resist the stromal fluid pressure applied to the epithelium. The mechanical balance equation across the epithelium, equation (5.14), written in the normal direction is Thus, the corneal anterior boundary could not be stably maintained without the contribution from inclined fibres and which leads to the following hypothetical question—if lamella inclination was absent, how would the boundary condition then be satisfied? This in essence is what occurs in LASIK where the inclined lamellae are largely removed by flap creation and ablation. In this case, the effective stress will become zero in the normal direction and the *in vivo* cornea must swell until a negative osmotic pressure is generated so that This is precisely what the model predicts in LASIK as described in §6.2.4 (although only approx. 3% swelling is needed to achieve ).

We have calibrated this model for the normal living cornea by assuming that GAG fixed charge is not depth-dependent; this leads to a hydration state that is also uniform with depth. It has been observed by Castoro *et al*. [76] that the anterior stroma may be slightly less hydrated than the posterior. The water content of these two stromal layers is measured to be 75% and 79% for rabbit, respectively [77]. Furthermore, biochemical studies on bovine corneas have shown that posterior stroma has a higher ratio of keratan sulfate to dermatan sulfate [76]. For human corneas, the depth-dependence of GAG types is unknown. Given that the hydration gradient is small through the corneal depth (less than 5%), it seems justifiable to employ a uniform fixed charge distribution as a first approximation. However, this can be easily modified as new information becomes available. Likewise, flexible GAG chains, the linear macromolecular structures formed from disaccharide units, may contribute an entropic elastic energy [23]. Each GAG chain is covalently bound at one end to the PG core protein and so is anchored to a collagen fibril (and thereby providing *fixed* charges). It has been proposed [9,10] that two or more GAG chains can link to form antiparallel duplex structures bridging between neighbouring collagen fibrils and can transmit entropic-based forces deriving from thermal motions. Such contributions can be modelled through appropriate free energy expressions but are considered to be a second-order effect and neglected in this study.

The stromal electrolyte free energy (equation (3.22)) is convex in the dilation *J* and represents a rational basis for modelling the osmotic behaviour of the tissue. Indeed we have shown that the model is capable of simulating massive swelling and the associated hydration change in the case of stroma immersed in an ionic bath. The predictions also show that swelling interaction with the collagen network leads to the observed anisotropy of swelling (primarily through the thickness direction) and inhomogeneity (most swelling occurs in the posterior cornea) which strongly supports the idea that the anterior cornea exhibits rigidity with respect to changes in hydration [29,31]. Measurements for *in vivo* swelling are not precise, but model predictions appear to agree qualitatively in all presented cases, and quantitatively for the rather subtle problem of swelling under increasing IOP. The model has provided predictions for swelling under Fuchs' dystrophy and should also be adaptable to modelling other diseases of the cornea such as macular corneal dystrophy, where the GAG concentration is found to be abnormally low.

The model supports the necessity of lamellae inclination for stability of the corneal anterior (refractive) surface [29,31]. When lamella inclination is removed synthetically by adjusting the fibre distributions, the model predicts an enormous increase in swelling in a stroma sample immersed in de-ionized water (§6.1.2). Theoretically, such swelling would become infinite, but in this model is limited by the small matrix shear stiffness. The structural role of lamella inclination has been examined here for the first time and the results suggest that this aspect of the collagen architecture is extremely important for understanding swelling behaviour and how the stromal fluid pressure is equilibrated.

Finally, it is remarked that clinical swelling (oedema), which results from hydration changes in the pathological cornea, can occur under a wide variety of circumstances. For example, during hypoxia due to contact lens wear, lactate ion production increases and results in oedema [19,78]. The underlying mechanism has been identified as the accumulation of lactate acid due to increasing anaerobic respiration [78]. Extension of the model to include these multiple species and their interactions through metabolic reactions for aerobic and anaerobic respiration is currently under development.

## Authors' contributions

X.C. and P.M.P. designed the research and developed the mathematical modelling of the electrolyte; S.J.P. and P.M.P. developed the mathematical modelling of the stromal elasticity; X.C. and P.M.P. performed the numerical analysis and interpretation; all authors wrote the paper jointly.

## Competing interests

We declare we have no competing interests.

## Funding

This research was supported by the Stanford University Bio-X Interdisciplinary Initiatives Program (IIP). X.C. received a Stanford Graduate Fellowship (SGF) from Stanford University. S.J.P. received a Bio-X Graduate Student Fellowship from Stanford University.

## Footnotes

↵1 For example, in application to the

*in vivo*cornea, collagen fibres in*Ω*_{0}are prestressed by the action of the IOP.↵2 Fluorophore-assisted carbohydrate electrophoresis.

↵3 If future studies produce a different conclusion regarding human stromal fluid pressure, including variation through depth, the value of

*Q*_{+}*Q*_{−}can be re-calibrated without any change to the general theory.↵4 This requires iterative solution for the configuration and collagen prestress that equilibrates the IOP—see §2.

↵5 The model analysed the biomechanical (elastic and swelling) response of the residual stroma. The flap was not analysed and the postsurgical anterior radius of the cornea was taken to be that of the postsurgical ablation profile radius, found by least-squares fitting, extended by a constant flap thickness of 150 µm.

↵6 The mean clinical and finite-element model values of

*R*_{pre}were 7.8 and 7.87 mm, respectively—the value*R*_{pre}= 7.8 mm was used in calculating SE for consistency.↵7 Because

*p*_{0}is constant at 15 mmHg (0.002 MPa), we have*P*_{os}=*p*−0.002, and the fringe plot can be interpreted in terms of osmotic pressure.↵8 Water permeability does not appear because at equilibrium water flux will be zero.

- Received March 18, 2015.
- Accepted June 17, 2015.

- © 2015 The Author(s)

Published by the Royal Society. All rights reserved.