## Abstract

The transparency of the human cornea depends on the regular lattice arrangement of the collagen fibrils and on the maintenance of an optimal hydration—the achievement of both depends on the presence of stromal proteoglycans (PGs) and their linear sidechains of negatively charged glycosaminoglycans (GAGs). Although the GAGs produce osmotic pressure by the Donnan effect, the means by which they exert positional control of the lattice is less clear. In this study, a theoretical model based on equilibrium thermodynamics is used to describe restoring force mechanisms that may control and maintain the fibril lattice and underlie corneal transparency. Electrostatic-based restoring forces that result from local charge density changes induced by fibril motion, and entropic elastic restoring forces that arise from duplexed GAG structures that bridge neighbouring fibrils, are described. The model allows for the possibility that fibrils have a GAG-dense coating that adds an additional fibril force mechanism preventing fibril aggregation. Swelling pressure predictions are used to validate the model with results showing excellent agreement with experimental data over a range of hydration from 30 to 200% of normal. The model suggests that the electrostatic restoring force is dominant, with the entropic forces from GAG duplexes being an order or more smaller. The effect of a random GAG organization, as observed in recent imaging, is considered in a dynamic model of the lattice that incorporates randomness in both the spatial distribution of GAG charge and the topology of the GAG duplexes. A striking result is that the electrostatic restoring forces alone are able to reproduce the image-based lattice distribution function for the human cornea, and thus dynamically maintain the short-range order of the lattice.

## 1. Introduction

The almost perfect transparency of the human cornea depends on the regular lattice arrangement of the collagen fibrils within the stromal lamellae and on the precise maintenance of an optimal hydration [1]. The achievement of both conditions depends on the presence of stromal proteoglycans (PGs) that associate with the collagen fibrils. The stroma, which comprises about 90% of the corneal thickness, is organized into lamellae (fibres organized in flat sheets) that contain collagen types I and V in the form of 25–30 nm diameter fibrils and small leucine-rich repeat PGs. The collagen fibrils within a lamella appear in parallel arrays following the direction of the lamella. In a lamella cross-section, the fibrils assume a quasi-regular packing arrangement with centre-to-centre interfibrillar distance of 40–60 nm; this arrangement is referred to as the fibril lattice. The short-range order of this lattice produces optical interference effects that render the tissue transparent [1–4].

The stromal PGs are sulfated and have linear sidechains of negatively charged disaccharide units called glycosaminoglycans (GAGs) that are bound covalently at one end to a PG core protein, which in turn is attached to a collagen fibril [5]. The charged GAG chains interact electrostatically with ions in the interfibrillar fluid to produce osmotic pressure. The resulting propensity of the stroma to swell by absorbing water is counteracted in the *in vivo* cornea by an active ion transport mechanism located in the endothelial cellular layer [6]. It is through the interplay of the GAG-based osmotic pressure and the active ion pumping mechanisms that hydration in the *in vivo* cornea is controlled. A measure of the overall swelling tendency of stroma *in vitro*, which includes osmotic and other effects, is provided by the tissue swelling pressure which has been measured by numerous investigators [7,8].

The GAGs also have an important but less-understood role in the maintenance of the fibril lattice. Keratan sulfate (KS), the predominant stromal GAG component, has been shown to be involved in modulating the fibril lattice by the knockout of *Chst5*^{1} in the mouse [9]. Scott [10] proposed that two or more GAG chains, originating at different core proteins on neighbouring fibrils, may form an antiparallel duplexed association that appears as a bridge-like structure spanning the interfibrillar distance in electron microscopy after staining, and that such structures may be capable of controlling the interfibrillar spacing. The open question of how negatively charged (and mutually repulsive) GAG chains might form durable antiparallel associations is discussed by Knupp *et al*. [11]. The idea of an entropic elastic interconnected fibril–GAG network was employed in the theoretical study on lattice forces and spatial order by Farrell & Hart [12] who proposed a model based on a sixfold symmetric elastic network. A sixfold symmetric structural model was also proposed by Müller *et al*. [13] in which GAG bridges (presumably formed by duplexing) were proposed to link next-nearest-neighbouring fibrils. However, recent three-dimensional electron tomography imaging by Lewis *et al*. [5] and Knupp *et al*. [11] shows that while GAG chains do appear to form bridging structures between two or more fibrils, these structures do not exhibit any approximation to a regular sixfold symmetric network or indeed any kind of regularity in spatial organization. Further, many GAG chains are not acting as bridges and these linear structures reach into the interfibrillar fluid at various orientations. All of this suggests that fibrillar bridges alone cannot mediate the lattice [5].

Maurice [2,14] speculated that repulsive forces must exist between adjacent fibrils to keep them separated and that the origin of these forces was related to the swelling force within the extracellular matrix. A conceptual model for the lattice forces has recently been described by Lewis *et al*. [5] in which an attractive force system between a fibril pair resulting from thermal fluctuations of antiparallel GAG duplexes is balanced by a repulsive force system which has an origin in osmotic pressure. In this model, the lattice dimensions will adjust to reflect a state of equilibration between the two opposing force systems, which depends on factors such as the number and type of GAG duplexes. The contractive entropic force from GAG duplexes remains an essential element of this model.

In the present work, a theoretical model for the forces that allow stable self-organization of the lattice is presented. An electrostatic-based fibril restoring force is introduced which results from a hypothesized local change in charge density induced by motion of a fibril relative to its neighbours in the lattice. It is also assumed that some GAG chains form antiparallel duplexes bridging between fibrils and these have been modelled by assuming that they possess the entropic elasticity of a worm-like chain. In developing the mathematical description, a conventional equilibrium thermodynamics approach [15] is adopted and Poisson–Boltzmann (PB) mean-field theory is used to describe interactions between GAG fixed charges and mobile ions. The tissue swelling pressure and lattice restoring forces are found from suitable derivatives of the free energy. An approximation for the electrostatic forces is also derived using Donnan theory and used in a coarse-grained dynamic model to describe the self-organization of the fibril lattice. The dynamic system, solved using a molecular dynamics technique, is used to explore the effects of randomness in the GAG chain organization by comparing the computed lattice radial distribution functions with an image-based estimate [16].

In a theoretical study on corneal transparency, Twersky [17] conjectured that collagen fibrils must be centred in a ‘tough adherent transparent coating’ that prevents fibril aggregation, an idea that has experimental support in the work of Fratzl & Daxer [18] whose measurements suggest that collagen fibrils are surrounded by a cylindrical coating consisting of PGs. The current model is used to explore the possibility that the collagen fibrils have a PG-dense coating. It is shown that such a coating will relatively have small influence on the tissue osmotic pressure and fibrillar restoring forces, except when the coatings of adjacent fibrils become close, at which point an additional and rapidly increasing electrostatic repulsive force is generated.

## 2. Ionization and charge properties of the glycosaminoglycans in corneal stroma

### 2.1. Predominant glycosaminoglycans

The predominant corneal GAGs are KS, dermatan sulfate (DS) and chondroitin sulfate (CS). These are linear unbranched polymers that are composed by repeating disaccharide units. In order to characterize the fixed charge carried by the corneal GAGs, we must estimate the average charge per disaccharide unit. The molecular structure of each GAG chain is shown in figure 1. The KS disaccharides are sulfated at both Gal and GlcNAc at the 6-carbon position. The CS disaccharides are sulfated at either the 4- or 6-carbon of the GalNAc residue, and are designated CS4 or CS6, respectively (figure 1). A detailed FACE^{2} [19] analysis discloses that only a portion of the available GAG monosaccharides is sulfated in the human cornea. In KS, the fraction of sulfated Gal and GlcNAc is 54% and 95%, respectively. For CS, the sulfation fraction is significantly lower. For CS4, the fraction of sulfated GalNAc is 28%, whereas for CS6 the fraction is only 8% (table 1).

### 2.2. Estimation of glycosaminoglycan ionization fractions

As sulfate groups carry one negative charge per group, the average charge per disaccharide unit for KS can be found from table 1 as *z*_{KS} = 0.95 + 0.54 = 1.49. Similarly, for CS/DS^{3} the value is *z*_{CS/DS} = 1.0 + 0.28 + 0.08 = 1.36. Introducing the GAG ionization fraction *f*, defined as the ratio of the average charge per disaccharide unit to the number of groups that are ionizable, it is found that *f*_{KS} = *z*_{KS}/2 = 74.5% for KS, and *f*_{CS/DS} = *z*_{CS/DS}/2 = 68.0% for CS/DS. The average ionization fraction of all corneal GAGs may then be estimated as
2.1where *x*_{KS} and *x*_{CS/DS} are the molar fractions of KS and CS/DS, respectively. Given the percent weight ratio between the two GAG species, *w*_{KS} = 65% and *w*_{CS/DS} = 30% [19,20] and the molecular weight (g/mol per disaccharide unit), *M*_{KS} = 464 g mol^{−1} and *M*_{CS} = 513 g mol^{−1} [21], the molar ratio of KS : CS/DS is estimated as *x*_{KS} : *x*_{CS/DS} = *w*_{KS}/*M*_{KS} : *w*_{CS}/*M*_{CS} = 2.40, which leads to *x*_{KS} = 70.5% and *x*_{CS/DS} = 29.5%. Substitution of these values in equation (2.1) produces *f*_{avg} = 72.6%.

### 2.3. Stromal glycosaminoglycan fixed charge

If a representative volume within the normally 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 volume is 2.2where

*e*is the unit charge. The

*average*stromal fixed charge density is then 2.3where is a measure of the GAG number density. However, in the current study, the charge density distribution must be characterized at the nanometre scale where it will be spatially non-uniform.

X-ray scattering studies have been used to reveal structural transformation of collagen fibrils under varying tissue hydration [18,22]. Fratzl & Daxer [18] performed x-ray scattering studies of collagen fibrils under varying hydration produced by drying the tissue. Their results suggest that the stromal PGs appear in relatively high density in the vicinity of the fibrils where they may form a charge-rich and water-binding PG-coating surrounding each fibril. The radius of the coating was measured to be *r*_{c} = 18.25 nm. They showed that the coating radius is insensitive to hydration over a wide range, suggesting that the water bound by the coating is released only at very low levels of tissue hydration. The existence of such a surface ultrastructure on the collagen fibrils is corroborated by image studies from Miyagawa *et al*. [23] and, to some extent, by Müller *et al*. [13]. In particular, Miyagawa *et al*. [23] had suggested that the PG-coating is primarily composed of KS. A theoretical study by Twersky [17] 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.

To accommodate the possibility that the collagen fibrils have coatings consisting of PGs, two classes of GAGs are introduced for modelling purposes (figure 2*c*): (i) an *interstitial* class, which spans the interfibrillar electrolyte volume and (ii) a fibril *coating* class, which spans a fibril coating volume. The stromal fixed charge *Q* is partitioned with a parameter such that the fixed charge *Q*_{inter} = *λQ* is associated with the interstitial GAGs, and the charge *Q*_{coat} = (1 − *λ*)*Q* is associated with the coating GAGs. The parameter *λ* facilitates modelling of the case with no coating GAGs, *λ* = 1, coating and interstitial GAGs, 0 < *λ* < 1, and the extreme (non-physiological) case of coating GAGs only, *λ* = 0.

The coating GAGs are not explicitly identified with KS or CS/DS. However, such identification can provide a lower bound on the physiological value of *λ*. The charge fraction contributed by KS and CS/DS is given as *Q*_{KS}/*Q*_{CS/DS} = *z*_{KS}*x*_{KS}/*z*_{CS/DS}*x*_{CS/DS} = 2.62. This implies that *Q*_{KS}/*Q* = 72.3% and *Q*_{CS/DS}/*Q* = 27.7%, which provides an estimation for the lower bound of *λ*_{l} = 0.28 based on the assumption that all the KS is contained in the PG coating [23].

## 3. Unit cell model

### 3.1. Thermodynamic equilibrium

The simplest possible geometric representation of a unit cell within a stromal lamella is that of an equilateral triangular prism with vertices at the centre of three fibrils (figure 2*c,d*). The cell contains a parallel array of collagen fibrils aligned with the lamella direction and arranged with hexagonal packing. The centre-to-centre interfibrillar spacing at normal hydration is denoted by , the fibril radius by *r*_{f}, and the PG-coating radius by *r*_{c} (figure 2*d*). The antiparallel GAG duplexes, which are formed by two or more GAG chains, bridge the neighbouring fibrils and are assumed to repeat continuously in the fibril axial direction with an average periodicity *h* (figure 3). The unit cell length is taken to coincide with *h*. Measurements of for the human cornea are summarized in table 2 and values for other unit cell parameters are summarized in table 3, which indicates that we have assumed and *h* = 25 nm as representative values from measurements.

For isolated corneal stroma immersed in an ionic bath, the total Gibbs free energy of the unit cell can be defined as
3.1where is the Helmholtz free energy, *P* and *V* are atmospheric pressure and total system volume (including the bath solution), respectively, and is the external work applied to the unit cell. At equilibrium, stationarity of the Gibbs energy implies:
3.2The *PV* term is essentially constant [30] and equation (3.2) reduces to . From the standard theory for electrolyte gels [15], the Helmholtz free energy has electrostatic, entropic elastic and mixing components, i.e.
3.3

For the cornea, it has been demonstrated that the free energy of mixing is very small [21,28,31] and this component is therefore neglected in the present study.

### 3.2. Electrostatic free energy

The unit cell domain is denoted and contains fibril and electrolyte regions The electrolyte domain is further partitioned into coating and interstitial regions (figure 2*d*). For a binary electrolyte (e.g. NaCl), the electrostatic potential *φ* over *Ω*_{electrolyte} is governed by the PB equation [32]
3.4where *ρ*_{GAG} is the non-uniform GAG fixed charge density, *ε*_{w} is the dielectric permittivity of the electrolyte solvent, *F* is the Faraday constant, *C*_{0} is the bath concentration at which the electrical potential is zero, *R* and *T* are the ideal gas constant and temperature, respectively. When *λ* ≠ 1 the charge density *ρ*_{GAG}(**x**) will be discontinuous at the *Ω*_{interstitial} − *Ω*_{coating} interface, whereas the potential *φ* and electric field **E** = −∇*φ* are required to be continuous across the interface.^{4}

The mean-field approximation of the free energy for a binary electrolyte can be expressed as a functional of *φ*, which is stationary at the solution to the PB equation [33] in the form:
3.5The first and second terms in equation (3.5) measure the electrostatic free energy of the GAG fixed charges and ionic charges, respectively, and the last term measures the contribution from the electric field.

### 3.3. Entropic elastic free energy

As noted in §1, Scott [10] proposed that two or more GAG chains, each attached to a fibril through its PG core protein, may form a duplexed association in which the linear chains are 180° apart, creating a bridge-like structure spanning the two fibrils (figure 3). Other imaging studies have supported this proposition, e.g. [5,13,34]. It is assumed that the antiparallel, duplexed GAG structures possess the free energy of a worm-like chain (WLC) [35], providing an entropic elastic behaviour arising from thermal fluctuations that will induce a contractive force between the bridged collagen fibrils. In this case, the free energy takes the form:
3.6where *N*_{b} is the number of GAG duplexes in *Ω*_{cell}, *k*_{b} is the Boltzmann constant, *L*_{p} is the persistence length of the GAG duplex, *L*_{c} is the contour length representing the length of the fully stretched duplexed chains, and *L _{i}* is the end-to-end distance of the

*i*th chain.

In order to characterize *L _{i}*, it is necessary to specify the locations of the PG binding sites where each protein core attaches to each fibril. Although it is known that the binding sites in the fibril axial direction occur at the a, c and d/e bands within the fibril D-period [36], the current model assumes that the antiparallel GAG duplexes are perpendicular to the collagen fibrils [36] with an axial periodicity

*h*(figure 3

*b*). However, the circumferential location of each of the PG binding sites on the two bridged fibrils can produce significantly different end-to-end bridging distances (figure 3

*a*). By assuming that each circumferential PG binding site is random for each fibril, the end-to-end distance

*L*may be replaced by the root mean square (r.m.s.) distance . For the nearest-neighbour fibril bridging 3.7where d

_{i}_{i}is the centre-to-centre distance of the bridged fibril pair. The distance function

*L*(

_{i}*θ*

_{1},

*θ*

_{2}) is illustrated in figure 3

*a*.

The r.m.s. end-to-end distance of a *free* WLC at normal hydration is denoted by *ℓ*^{0} and is assumed to be given by equation (3.7) with d* _{i}* replaced by the normal hydration distance Using the relationship between

*ℓ*

^{0}and

*L*

_{p}and

*L*

_{c}[37], it follows that 3.8This expression allows the contour length

*L*

_{c}to be determined in terms of

*L*

_{p},

*r*

_{f}, and . Table 3 gives the values used in this study. For the nearest-neighbour bridging, which leads to

*ℓ*

^{0}= 37.83 nm and

*L*

_{c}= 90.88 nm. This implies that duplexed GAG chains are approximately 42% stretched at normal hydration.

## 4. Stromal swelling pressure

### 4.1. Thermodynamic definition

Stromal swelling pressure *P*_{s} is defined as the pressure applied on a sample of the tissue in an ionic bath that is required to prevent it from swelling. The dependence of *P*_{s} on hydration, bath ionic strength, temperature and pH have been measured in numerous species including human, and by various means. The application of the current model to the corneal swelling pressure provides an important validation case. The results for *P*_{s} are also used directly in the lattice force model presented in §§5 and 6.

Consider a representative volume of stromal tissue at normal hydration. This volume will contain keratocyte cells and is distinct from the unit cell *Ω*_{cell} previously introduced. Let the tissue volume be subject to a dilation where *V*_{stroma} will be referred to as the current volume. If it is assumed that the total collagen, keratocyte and PG content are preserved under the dilation while the incompressible electrolyte fluid is free to enter or leave the volume during dilation,^{5} then at fixed temperature and bath ionic concentration, the external work needed to achieve the volume *V*_{stroma} is
4.1At equilibrium, it follows from equation (3.2) that
4.2where and are the electrostatic and entropic elastic swelling pressure components, respectively. The swelling of the corneal stroma has been observed to be predominantly unidirectional in the thickness direction [31,38] and experimental measurements of the interfibrillar spacing of the swollen cornea suggest that the fibril lattice repacks itself bidirectionally [38]. In the swollen stroma it is therefore reasonable to assume that the fibrils will assume a lattice with equal interfibrillar spacing. Then the interfibrillar centre-to-centre distance *d*_{c} under volume dilation *J* can be determined as
4.3

### 4.2. Stromal charge density

In the stroma, collagen fibrils and keratocyte cells occupy volume that is excluded from the electrolyte fluid. Table 4 lists the approximate weight fraction of each constituent of the stroma [27]. Denoting the collagen and keratocyte volume fractions as *ϕ*_{col} and *ϕ*_{ker}, respectively, and denoting the volume fraction of the remaining constituents as *ϕ*_{electrolyte} = 1 − *ϕ*_{col} − *ϕ*_{ker} (primarily water), and using the data of Meek & Leonard [24], we have estimated *ϕ*_{col} = 24.9%, *ϕ*_{ker} = 11.2% and *ϕ*_{electrolyte} = 63.9%.

Assuming that keratocyte cells dilate with the tissue dilation [39] and assuming collagen fibrils are non-swelling [22,31,38], we find two charge densities associated with *Q*_{inter} and *Q*_{coat} under macroscopic dilation *J* to be as follows:
4.4and
4.5where is given by equation (2.3). The above expression has used the coating volume fraction *ϕ*_{coat} defined by For consistency with the findings of Fratzl & Daxer [18], the coating volume does not change with dilation *J* such that

The non-uniform (piecewise constant) GAG fixed charge density *ρ*_{GAG} appearing in the PB equation (3.4) can now be defined as
4.6where the integer-valued *n*(**x**) is given by
4.7

### 4.3. Poisson–Boltzmann-based numerical solution for *P*_{el}

A sequence of dilated unit cells, parametrized by the dilation *J*, is constructed with interfibrillar centre-to-centre distance given by equation (4.3). The numerical method proceeds by solving the nonlinear PB equation (3.4) with equation (4.6) for the electrostatic potential *φ*(*J*) by the finite-element method. The electrostatic free energy is then evaluated using equation (4.2) for each member in the unit cell sequence. The electrostatic contribution to the swelling pressure is determined for any dilation *J* from equation (4.2) and where is approximated numerically by central difference. As *J* is reduced, the coating domains will overlap and the charge density in the overlap region is doubled as indicated by *n*(**x**) above. Electrostatic potential contours for *J* ∈ [0.42,1.0] are illustrated in figure 9*a–d* for typical stromal parameters, including the case of coating overlap.

### 4.4. Donnan-based approximation for *P*_{el}

The Debye length is the characteristic distance in an electrolyte over which the potential of a charged surface decays by 1/*e*. For the human stroma, this is estimated to be
4.8Because this value is small compared to the unit cell dimensions, it is reasonable to consider an approximation based on Donnan potentials associated with the *ρ*_{inter} and *ρ*_{coat} charge densities. The piecewise constant Donnan potential is given by [32]
4.9where *ρ*_{GAG} is given by equation (4.6). Although satisfies the PB equation (3.4), it will not satisfy the necessary continuity conditions at the *Ω*_{coat} − *Ω*_{interstitial} interface. Introducing equation (4.9) into (3.5) and using equation (4.2), the Donnan-based electrostatic contribution to the swelling pressure is given as (appendix A)
4.10The second term on the right-hand side of equation (4.10) provides the correction for the effect of the coating. If the coating charge is set to zero (*λ* = 1) and the keratocyte volume is neglected (*ϕ*_{ker} = 0), the above swelling pressure reduces to the classical Donnan form with the collagen volume exclusion effect included [40]. On the other hand, if all the GAG fixed charge is in the coating (*λ* = 0), then *ρ*_{inter} = 0 and is zero for all dilation *J*, consistent with the assumption that the coating volume does not change with tissue dilation. The accuracy of equations (4.9) and (4.10) is examined in §7.

### 4.5. Entropic elastic contribution

The GAG duplex r.m.s. end-to-end distance *ℓ _{i}*, given by equation (3.7), may be expressed as a function of

*J*by substituting equation (4.3) in equation (3.7) with

*d*=

_{i}*d*

_{c}resulting in 4.11Use of the chain rule allows

*P*

_{str}(

*J*) to be determined from equations (4.2), (3.6) and (4.11) as 4.12where is a measure of duplexed GAG bridging density in the corneal stroma. For the parameter values given in table 3, the entropic elastic swelling pressure component

*P*

_{str}remains negative (contractive) for all

*J*> 0.146.

## 5. Interfibrillar lattice restoring forces in an ideal hexagonal cell

### 5.1. Rationale and thermodynamic definition

Consider an individual fibril displaced to a new position in the lattice while all other fibrils maintain their positions. Linear GAG disaccharide chains are associated with each fibril and attach to that fibril at one end only via the PG core protein. It is hypothesized that the attached GAG chains will displace with their supporting fibril. Because the longer chains reach into the interfibrillar fluid, a local change in GAG fixed charge density is created with an increase in advance of the displacement and a reduction behind (figure 4). In our model, the interstitial GAG class introduced in §2 are assumed to be involved in this process, whereas the shorter coating class GAGs will not, until the coatings become close. Additionally, GAG chains that have formed interfibrillar antiparallel, duplexed structures will experience changed entropic forces resulting from changes in end-to-end distances. For the purpose of fully characterizing the two restoring forces, the analysis in this section is based on an ideal hexagonal unit cell with perfect duplexed GAG bridging based on nearest-neighbour connectivity. In §6, the important effects of random GAG chain topology is considered.

The interfibrillar lattice restoring force **f*** _{l}*(

**r**) is defined as the force that must be imposed on a fibril to maintain its current position

**r**. The external work required to move the fibril from its lattice position

**r**

_{0}to displaced position

**r**, while fixing the position of other fibrils in the lattice is given by 5.1From equations (3.2) and (3.3),

**f**

*is given as 5.2where and are defined as the electrostatic and entropic elastic components of the restoring force.*

_{l}### 5.2. Electrostatic restoring force

The region between the target (central) fibril and its six neighbours in the lattice is divided into triangular prism subcells with initial volume The central fibril is then displaced to position **r** causing a change *Δ**V _{k}* in the volume of each subcell (figure 4

*b*). Let measure the

*k*th subcell change in volume due to the fibril displacement to

**r**and set

The GAG fixed charges within subcells are assumed to be conserved during fibril displacement. The fixed charge density *ρ*_{inter} given by equation (4.4) is modified to for the *k*th subcell as
5.3Since the volume of the hexagonal cell is fixed (*J* is independent of **r**), it follows that *Σ*_{k}*Δ**J _{K}* = 0. Figure 5 shows for each subcell when the central fibril is displaced vertically and indicates that charge density increases in advance of the displacement and reduces behind it. Note that

*ρ*

_{coat}and

*ρ*

_{GAG}and are still given by equations (4.5) and (4.6), respectively.

The electrostatic free energy of the cell can be evaluated by summing over subcells
5.4where takes the precise form of equation (3.5), but with the *k*th subcell volume as the integral domain. An energy sequence *i* = 1, … ,*n* is determined numerically using equation (3.5), after solving the PB equation on the subcells. The force **f**_{el}(**r**) is then obtained from equation (5.2) using central difference. This scheme considers all electrostatic effects, including coating overlap.

### 5.3. Donnan approximation

We make the approximation , and an analytical formulation of **f**_{el} is found by
5.5where is the *k*th subcell Donnan osmotic pressure component from equation (5.5). This result shows clearly how the osmotic pressure in the subcell serves as the source for the interfibrillar restoring force. The vector ∂*J _{k}*/∂

**r**is directed towards the centre of the lattice for all

*k*(see appendix B for details); thus equation (5.5) describes a set of centrally oriented forces with differing magnitudes depending on the subcell osmotic pressure. Clearly, when the fibril is in its lattice position, the osmotic pressure is uniform and the electrostatic force vanishes identically. Away from the lattice position, the electrostatic force is always a

*restoring*force. An interpretation of the electrostatic restoring force as given by equation (5.5) is provided in figure 4. Equation (5.5), which evaluates the electrostatic restoring force without solving the PB equation, is employed in a fibril dynamics model described in the next section.

### 5.4. Entropic elastic restoring force

Consider the centre of the displaced fibril at location **r**, and the centre of a neighbouring fibril in the hexagonal cell at location at **r*** _{i}*, with centre-to-centre distance Observing from equation (3.6) that a GAG duplex between two fibrils will have entropic free energy

The relationship between the bridging duplex end-to-end r.m.s. distance *ℓ _{i}* and

*d*is given by equation (3.7). Then for all six nearest-neighbour bridges, we have, from equation (5.2) 5.6In the above expression, −∂

_{i}*d*/∂

_{i}**r**

*=*

_{i}**v**

*, a unit vector pointing from*

_{i}**r**to

**r**

*, indicating the direction of the contractive force on the displaced fibril.*

_{i}## 6. Lattice dynamics model with an imperfect glycosaminoglycan topology

A dynamic model is introduced to describe fibril movements and lattice self-organization. Consider a realistic fibril lattice, with many collagen fibrils, in a transient state; in analogy with molecular dynamics theory, any imbalance in the interfibrillar forces on fibril *i* will induce a fibril displacement
6.1where *m _{i}* is the fibril mass and

**a**

*is the acceleration. To incorporate the force computations from the previous sections, the lattice is divided into equilateral triangular subcells defined by the nearest-neighbouring fibrils. In order to efficiently compute the electrostatic force, equation (5.5) is used to avoid solving the nonlinear PB equation at each time step. In addition, a soft-sphere potential [41] is introduced to prevent fibril–fibril interpenetration as fibrils have a finite size. The soft-sphere potential can be described by a generalized inverse power law (GIPL) 6.2where*

_{i}*r*refers to the distance between fibril pairs and σ

_{ij}*= 2*

_{r}*r*that corresponds to the fibril diameter. A force field based on a value of

_{f}*n*= 32 was used in this study and gives an extremely short-range force field such that the force magnitude drops by more than 95% after two fibrils are separated by 3 nm.

To describe a realistic (i.e. random) GAG chain system [5], two types of imperfection are introduced: (i) the stromal charges are not equally allocated to each subcell and (ii) not all bridging sites are occupied. In this approach, the reference charge densities in subcells have been randomly assigned which obeys the Gaussian distribution , where σ* _{p}* is the standard deviation. To describe a defective GAG bridging network, each bridging site has been assigned a probability

*p*

_{b}for the existence of GAG duplex.

The dynamic system has been solved using a typical molecular dynamics routine implemented in MATLAB. As is customary, all the physical parameters have been reduced to dimensionless (table 5). The velocity-Verlet algorithm is used for the time integrator, and the standard cell subdivision and neighbour-list algorithms [42] have been implemented to reduce the interaction computations. The results have been interpreted by computing the radial distribution function (RDF), which is a measure of the average number density of fibril centres at a given radius *r* from every other fibril normalized by the bulk fibril number density. A structure with a long-range order will have many sharp peaks in its RDF corresponding to the distances between the particles with high probability, whereas the RDF for a purely random structure will be constant with a value of unity. For structures with short-range order, such as the fibril lattice in the corneal stroma, the RDF has multiple peaks with decaying magnitude over a fixed interval and then quickly converges to unity thereafter [16]. The classical approach of evaluating the two-dimensional RDF is given as [42]
6.3where *A* is the system area and *h _{n}* is the number of fibril pairs (

*i*,

*j*) for which (

*n*− 1)

*Δ*

*r*≤

*r*≤

_{ij}*n*

*Δ*

*r*and

*r*= (

_{n}*n*− 1/2)

*Δ*

*r. N*is the total number of the fibrils in the system.

_{m}## 7. Results

### 7.1. Swelling pressure

#### 7.1.1. Dilation and hydration

For convenience in interpreting results, stromal hydration *H*, defined as water weight per unit dry weight, may be approximately^{6} related to the tissue volume dilation *J* by
7.1such that normal (physiological) hydration corresponds to *J* = 1 and *H*^{0} = 3.2 [4]. The linearity is consistent with measurements [8,43].

#### 7.1.2. Keratocyte and collagen electrolyte volume exclusion effects

Predicted electrostatic swelling pressure component *P*_{el}, given by equation (4.10), is plotted against volume dilation *J* in figure 6*a* for three cases: no volume exclusion, only collagen volume exclusion and collagen and keratocyte volume exclusion. The volume exclusion effects arising from these two stromal components are important for all *J* < 1 corresponding to hydration levels lower than normal. For dilations *J* < 0.5 (hydration less than 33% normal), the volume exclusion effect is significant, increasing *P*_{el} by almost an order. Collagen volume exclusion is the most important because the keratocyte volume is smaller and also because it dilates with the tissue.

#### 7.1.3. Collagen proteoglycan-coating effects

The effect of the coating charge partition *λ* on *P*_{el} is illustrated for three values of *J* in figure 6*b. P*_{el} is relatively insensitive to *λ* for 0.5 < *λ* ≤ 1, and reduces significantly in the range 0 ≤ *λ* < 0.5 as progressively more of the total charge is allocated to the coating region and less to the interfibrillar region. In the extreme case where all charge is located in the coating region, *λ* = 0, the predictions converge to zero swelling pressure, reflecting the condition that the coating volume does not change with hydration [18].

#### 7.1.4. Swelling pressure prediction versus experimental measurements

Predicted swelling pressure *P*_{s}, given by equation (4.2), is compared with experimental measurements [7,8] in figure 7. The average stromal charge density used in these calculations was found by requiring the predicted *P*_{s} to exactly match the measured swelling pressure at *J* = 1 (normal hydration) and *λ* = 0.65. Olsen & Sperling [8] provide a power-law fit of their data where *t* is the stromal thickness, and this was used to determine the physiological swelling pressure as 84.35 mmHg.

In figure 7*a*, computed and measured swelling pressure are compared over a range from swollen *J* = 1.5 (168% normal hydration) to compressed *J* = 0.55 (39% normal hydration) for three values of *λ*. The curve for *λ* = 0.65 agrees by design with the experimental data at *J* = 1, and shows excellent agreement with the data over the entire range of hydration. The other two curves corresponding to *λ* = 1 and 0.3 show the effect of changing the charge partition; the solution for *λ* = 0.3 deviates significantly from the experimental data. In figure 7*b*, the predicted swelling pressure is shown for three values of interfibrillar spacing (from 43 to 60 nm) with other unit cell parameters being fixed.

To illustrate the effect of coating overlap on swelling pressure, two coating radii *r*_{c} = 18.25 nm and 21.56 nm [17,18] have been modelled at very low hydration and the results are summarized in figure 7*c*,*d*. For both radii, the swelling pressure exhibits a very rapid increase with reducing hydration, owing to fibril coating overlap, with values that are in good agreement with the experimental data. The range of hydration for each coating is terminated when the coating reaches the collagen core.

In figure 8, the ratio between the electrostatic *P*_{el} and entropic *P*_{str} components of *P*_{s} is given and indicates that at normal hydration the magnitude of *P*_{el} is at least an order greater than the magnitude of *P*_{str}. At lower hydration, the ratio increases quickly, exceeding two orders at 60% normal hydration. When swollen to about 300% normal hydration, they become comparable.

#### 7.1.5. Assessment of the Donnan potential-based approximate solution

Contour plots of the electrostatic potential *φ*, obtained by finite-element solution of the PB equation (3.4) are shown in figure 9*a*–*d* for four hydration values. In order to compare *φ* with the Donnan approximation given by equation (3.4), the profile of the two potentials is plotted along the fibril centre-to-centre line in figure 9*e*–*h*. At normal hydration, figure 9*e* shows that agrees reasonably well with *φ*. As the hydration is reduced, the coating domains approach figure 9*f*,*g* and then overlap figure 9*h*. As can be anticipated, the four representative hydration values illustrated in figure 9*e*–*h* clearly indicate that the Donnan approximation cannot accurately represent *φ* when the coatings are close or overlapping. A comparison of the value of *P*_{el} computed from the PB solution *φ* using equation (4.2) and from the Donnan potential using (4.10) is shown in figure 10*a*,*b*. The Donnan solution, which takes into account the coatings but not their interactions, is very accurate until the coatings become close, at which point the Donnan approximation does not share the rapid increase in swelling pressure predicted by the PB solution.

### 7.2. Interfibrillar restoring forces in an ideal hexagonal cell

In order to characterize and compare the electrostatic and entropic restoring forces **f**_{el} and **f**_{str}, developed in §5, we consider an ideal hexagonal cell with fibrils at each vertex (figure 4). The hydration is taken to be normal with *J* = 1, and GAG fixed charges are allocated equally in each subcell. The resulting electrostatic restoring force magnitude *f*_{el} is computed from equation (5.2) and depicted in figure 11*a*,*b* for two fibril coating radius values *r*_{c} of 18.25 and 21.56 nm. The force–displacement relationship exhibits two quasi-linear regimes with a transition where the coatings begin to interact. This result quantifies the effect of the PG-coatings on the fibrils as a mechanism preventing fibril aggregation. For example, with *r*_{c} = 18.25 nm and *λ* = 0.65, the restoring force on the fibril at a displacement of 10 nm is less than 1 mN m^{−1}; when the displacement is extended to 20 nm, the force increases to approximately 3 mN m^{−1}.

In order to characterize the isotropy of this force, the central fibril was displaced around a circle while the vertex fibrils remained fixed. The results for the magnitude of the restoring force are shown in figure 11*c*,*d* with two circles selected for each *r*_{c}. The smaller circle produces almost perfectly isotropic results as indicated by the flat force–angle curves in figure 11*c*,*d*. On the extended circle, the coatings interact and produce a degree of anisotropy.

For the purpose of characterizing the entropic restoring force associated with the antiparallel GAG duplexes, it is assumed in this ideal situation that all six nearest-neighbour bridging sites are occupied producing a perfect network with hexagonal symmetry. The results, shown in figure 12, indicate an essentially linear force–displacement relationship. At a fibril displacement of 10 nm, the entropic elastic restoring force is approximately 0.03 mN m^{−1} and is thus two orders lower than the electrostatic restoring force. Note, however, that at lattice location **r** = **r**_{0}, we have **f**_{el} = −**f**_{str} = 0, which reflects the ideal symmetry of this case.

Table 6 summarizes the electrostatic restoring stiffness, defined as the fibril restoring force per metre of fibril divided by the fibril displacement. The two quasi-linear regimes were considered and the electrostatic restoring stiffness in each regime was computed by least-squares regression. For example, with *r*_{c} = 18.25 nm and *λ* = 0.65, the electrostatic stiffness was 69 kPa for a fibril displacement less than 10 nm and 483 kPa for a fibril displacement greater than 16 nm (a sevenfold increase). In contrast, the entropic restoring stiffness is only 3 kPa.

### 7.3. Fibril dynamics model for lattice self-organization

The proposed lattice restoring force mechanisms are studied further using the fibril dynamics model which features an imperfect GAG chain topology modelled by introduction of random variability in both fixed charge density and GAG chain elasticity. The fibril dynamics model is based on a two-dimensional periodic box containing 256 (16 × 16) collagen fibrils. The system is initialized with a perfect hexagonal lattice and zero initial velocities. The hydration is taken to be normal (*J* = 1) and charge partition *λ* = 1. The simulation has been run long enough to ensure equilibrium and the results are found independent to the system size and choice of time step. The RDF *g*(*r*) is obtained from equation (6.3) at each time step and averaged with respect to time.

Figure 13*a* shows the result of *g*(*r*) for the case where only the electrostatic restoring force **f**_{el} is considered, using three standard deviation values of 0.1 and 0.2. The charge samplings in the subcells for these three cases are summarized in the bar histogram in figure 13*b*. When and 0.1, the computed *g*(*r*) closely matches measurements on electron micrographs of human corneas [16] and indicates short-range order in the lattice such that correlations between any two fibrils becomes nearly random after *r* = 200 nm. As *σ*_{ρ} increases to , a more irregular lattice with a short-range order up to 80 nm results. The average interfibrillar distance, which corresponds to the location of the first peak in *g*(*r*), is maintained at in all three cases.

The influence of duplexed GAG bridging elasticity on the fibril lattice order has been studied by computing *g*(*r*) for four imperfect networks with **f**_{el} excluded. The results are summarized in figure 13*c*. For a nearly perfect network (i.e. *p*_{b} = 0.9), the average interfibrillar distance is maintained as but the peaks of the computed *g*(*r*) do not decay as fast as the measurement with increasing *r*, indicating that a long-range order lattice is predicted in this case [12]. For the networks with more than 20% bridges missing (i.e *p*_{b} < 0.8), the location of the first peak in *g*(*r*) is immediately reduced to 31 nm, which is close to the value of fibril diameter and implies that fibril aggregation is occurring in the lattice. The results indicate that the duplexed GAG structures alone cannot maintain the interfibrillar distance when the network is imperfect.

Figure 13*d* shows the computed *g*(*r*) when both **f**_{el} and **f**_{str} are present, with and *p*_{b} varying from 0 to 0.9. The peak locations are maintained in all these cases and clearly the lattice structure becomes *insensitive* to the duplexed GAGs bridging topology when **f**_{el} is present.

## 8. Discussion

The proposed collagen–PG interaction model governing the fibrillar lattice was validated by predicting the stromal swelling pressure. The results, which showed good agreement with experimental data over the entire tested range of hydration, were used to establish a reasonable value for the coating charge parameter of *λ* = 0.65, although there is not a high sensitivity to this parameter with regard to swelling pressure. This value lies midway between the limit of no coating charge, *λ* = 1 and the lower limit *λ* = 0.28 established in §2. It may be possible to explore this parameter more directly by using three-dimensional electron tomography images such as those in [5,11]. Experimental results for swelling pressure at extremely low hydration, which are on the order of 2500 mmHg at a hydration of around 20% of normal, cannot be matched by classical Donnan theory. The proposed model modifies the Donnan result with two additional inputs: the electrolyte volume exclusion provided by the collagen and keratocyte stromal constituents, and the effect of fibril PG-coating interaction.

In the present work, we have used a theoretical model to describe mechanisms by which movement of a fibril out of its lattice position produces electrostatic and entropic elastic restoring forces. In this model, the GAG chains generate the underlying tissue osmotic pressure and they are also responsible for a local perturbation in that pressure when a fibril, along with its GAG cargo, is displaced from its lattice position. All GAG chains can contribute to such a perturbation, regardless of whether they are individual chains or are part of an antiparallel, duplexed association bridging between adjacent fibrils. It has been shown that in the dynamic situation where a fibril is out of its lattice position, the electrostatic and entropic elastic restoring forces: (i) are both restoring forces, (ii) essentially operate independently and not as an antagonistic mechanism in which they are required to oppose each other and (iii) the electrostatic forces alone can restore the lattice. In addition, according to our calculations, the resulting electrostatic-based restoring forces subordinate the entropic elastic forces of the antiparallel GAG duplexes and are relatively insensitive to the random organization of the GAGs.

As described in §4, stromal swelling pressure has both osmotic (positive) and GAG entropic elastic (negative) components; it depends strongly on hydration and has a considerable physiological value of approximately 60 mmHg. However, common experience informs us that healthy living corneas do not swell; the hydration (and thus swelling) of the stroma is precisely controlled to be 78% water weight by means of an active ion transport mechanism located in the bounding endothelial layer (figure 2*a*) [14]. Hydration values exceeding optimal have been shown to result in corneal opacity. The active mechanism that counters the swelling pressure in the living cornea is implicit in the presented model because the stromal volume (for any hydration level) is maintained by forces that precisely balance the swelling pressure. Further, the fibril lattice theory presented is valid for the living system when the local perturbations in the osmotic pressure are based on the physiological hydration (as was done for the lattice dynamics model in §§6 and 7).

It has also been shown that, according to the model, if the fibrils do have a GAG-dense coating, and our results do not preclude that, then when fibril pairs approach the point where the coatings are close or overlapping, an additional restoring force mechanism is activated and the restoring stiffness (force per unit displacement) increases sixfold over the value before the coatings engaged. In this way, the model gives a quantitative assessment of the coating effectiveness with regard to preventing fibril aggregation.

In our modelling of randomness of the GAG organization, it is always assumed that the duplexed GAG chains are perpendicular to the fibril axis [13]. In the light of recent findings, this is not realistic because many are tilted with respect to the fibril axis [5]. Such considerations could be included in the theory, but presumably would tend to further reduce force estimates because only a component of the force vector will act as a restoring force. In regard to the relatively small bridging forces, it is interesting to note that in recent studies of other connective tissues such as tendon, cartilage and ligament, interfibrillar PG-based filaments have been found to be unlikely to have any mechanical role [44].

The dynamic analysis was based on a lattice with a random GAG organization. The randomness of the GAGs was reflected in both the spatial distribution of charge and the topology of the bridges. A striking result is that the electrostatic restoring force alone was able to reproduce the radial distribution function of the fibrils for the human cornea, confirming the short-range order of the lattice; no input from the GAG duplexes was needed. The converse did not hold. The dynamic case with bridging entropic forces alone actually predicted fibril aggregation on a defective network because the entropic force is essentially contractive.

The lattice dynamics study presented was performed at normal hydration. At a swollen state, the cornea can become opaque and manifest stromal ‘lakes’ [1,27], which are regions in which the lattice is profoundly disturbed with reduced numbers of collagen fibrils. Although results have not been presented here, the proposed model is applicable to modelling these conditions. In the current model, PG/GAGs have been assumed to be distributed uniformly throughout the stroma. Biochemical studies have shown that the anterior corneal stroma has a higher ratio of DS to KS [45]. The proposed theory could be adapted to model how the corneal hydration profile may be affected by this distribution. Another potential application concerns the loss of KS which occurs in corneas with macular corneal dystrophy [5,26] or to genetically mutated corneas [9], which leads to a reduction in interfibrillar distance and a more chaotic lattice due to the lack of PG regulation.

## Funding statement

This research was supported by the Stanford University Bio-X Interdisciplinary Initiatives Program (IIP), which is gratefully acknowledged.

## Acknowledgements

X.C. received a Stanford Graduate Fellowship (SGF) from Stanford University.

## Appendix A. Derivation of the Donnan-based electrostatic swelling pressure component *P̂*_{el}

The piecewise constant Donnan potential is written by combining equations (4.6)–(4.8) as

A1

where the charge densities *ρ*_{inter}(*λ,J*) and *ρ*_{coat}(*λ*) associated with the two GAG classes are given by equations (4.4) and (4.5), respectively. Substitute equation (A 1) into equation (3.5), the Donnan electrostatic free energy approximation is given as
A2where *I*_{inter} and *I*_{coat} are the two integrands in the interstitial and coating regions:
A3and
A4

The Donnan-based electrostatic contribution to the swelling pressure is then given by equation (4.2) with replaced by the Donann approximation , equation (A 2). Using equations (A 1), (A 3) and (A 4), it is observed that A5and A6therefore can be derived as A7

## Appendix B. Derivation of the Donnan-based electrostatic interfibrillar restoring force

From equation (5.5), the electrostatic interfibrillar force **f**_{el} is expressed under Donnan approximation as
B1The volume dilation on the *k*th subcell *J _{k}* is expressed as
B2and
B3where is the current volume of the subcell and measures the subcell cross-section area and is a function of

**r**,

**r**

_{k}_{1},

**r**

_{k2}, which are the coordinates of fibrils forming the

*k*th subcell. The analytical form of and its derivative are given as B4and B5where

**r**= (

*x*,

*y*) and refers to the sign of the term in the absolute operator of equation (B 4). Combining equations (B 1), (B 3) and (B 5) gives B6

## Footnotes

↵1 Carbohydrate sulfotransferase 5 enzyme.

↵2 Fluorophore-assisted carbohydrate electrophoresis.

↵3 It is assumed that the GlcUA/IdoUA moiety always carries one negative charge.

↵4 Generally it should be the dielectric displacement

**D**=*ε***E**that is continuous across the interface. Here we assume that the PG-coating domain*Ω*_{coating}is ion-penetrable and has the same dielectric permittivity as the interstitial region*Ω*_{interstitial}.↵5 As is the case in swelling pressure experiments, e.g. [7,8].

↵6 Here we ignore the amount of bound water in the collagen fibrils [38].

- Received June 11, 2013.
- Accepted July 8, 2013.

- © 2013 The Author(s) Published by the Royal Society. All rights reserved.