A discrete-to-continuum model for the human cornea with application to keratoconus

Abstract

We introduce a discrete mathematical model for the mechanical behaviour of a planar slice of human corneal tissue, in equilibrium under the action of physiological intraocular pressure (IOP). The model considers a regular (two-dimensional) network of structural elements mimicking a discrete number of parallel collagen lamellae connected by proteoglycan-based chemical bonds (crosslinks). Since the thickness of each collagen lamella is small compared to the overall corneal thickness, we upscale the discrete force balance into a continuum system of partial differential equations and deduce the corresponding macroscopic stress tensor and strain energy function for the micro-structured corneal tissue. We demonstrate that, for physiological values of the IOP, the predictions of the discrete model converge to those of the continuum model. We use the continuum model to simulate the progression of the degenerative disease known as keratoconus, characterized by a localized bulging of the corneal shell. We assign a spatial distribution of damage (i.e. reduction of the stiffness) to the mechanical properties of the structural elements and predict the resulting macroscopic shape of the cornea, showing that a large reduction in the element stiffness results in substantial corneal thinning and a significant increase in the curvature of both the anterior and posterior surfaces.

Keywords: collagen lamellae; corneal mechanics; discrete-to-continuum asymptotics; keratoconus; multiscale modelling; proteoglycan matrix.

Figure 1.

The unloaded geometry. Zoom onto the macroscale geometry of an idealized two-dimensional corneal slice (created with BioRender.com). The unloaded-configuration cornea is shown in blue, the key parameters are indicated in red and the green dotted lines (circles) depict the curves of constant Φ ($\tilde{R}$). The unloaded configuration is then discretized (N = 4, M = 80, γ = 20). The panel on the bottom right presents a zoom onto the unit cell with the dimensional lengths of the elements (denoted with tilde) dependent on the radial position (index i). Terms in red represent corresponding quantities in the continuum limit (N → ∞) of the dimensionless model. Note that while the unit cell for a finite N and M forms an isosceles trapezoid, in the continuum limit this becomes a rectangle.

Figure 2.

Healthy cornea loaded with IOP. The central panel shows the loaded configuration ($\widetilde{\hspace{0.17em}p}=2$ kPa) of the discrete model with N = 4 and M = 80. Note that the colour represents the axial force in the elements with positive (negative) values indicating tension (compression). Schematic in the left panel depicts how continuum IOP is transformed into discrete forces at posterior hinges. Panels on the right present zoomed-in steady-state force distributions through the corneal thickness at the apex (top) and the limbus (bottom).

Figure 3.

For a healthy cornea under physiological IOP (2 kPa), node displacements are small (a) and the axial forces in individual elements decrease as N → ∞ (b–d). The discrete axial forces in the lamellae can be divided by an appropriate cross-sectional area to get equivalent lamellar stresses. These stresses for large enough N (e–g) agree well with the circumferential stress ${\tilde{T}}_{\mathit{\Phi }\mathit{\Phi }}$ in the continuum limit (h), because the lamellar segments are the stiffest elements. Note that while ${\tilde{T}}_{\mathit{\Phi }\mathit{\Phi }}$ in the continuum model exceeds the displayed range in small regions near domain corners, we for ease of comparison restrict the colourbar to the same range as for the discrete model.

Figure 4.

Details of posterior (a) and anterior (b) surfaces near the apex for a healthy cornea inflated to 2 kPa IOP. Using M = 20N, the discrete profiles for N = 2 (blue squares), N = 4 (red pluses), N = 8 (green diamonds), N = 16 (cyan circles) and N = 32 (purple crosses) converge to the continuum predictions (solid black curves). Note that in order to make visible the otherwise small differences between the discrete and continuum solutions, we used different scales for X and Y axes (all indicated length units are still in mm). Panels on the right depict the predictions for apex displacement from discrete and continuum models. Panel (c) shows that the prediction of the discrete model for increasing N (blue squares) approaches that of continuum model (horizontal black line). Panel (d) documents that the absolute error, defined as the difference in apex displacement between the discrete and the continuum model, decreases to 0 as O(1/N).

Figure 5.

Predicted diseased corneal shapes using the continuum model with the parabolic (ξ = 2 in (4.2)); (a–d), quartic (ξ = 4; e–h) and sextic (ξ = 6; i–l) damage profiles and varying value of central damage ${\mathcal{D}}_{\max }=0.8$ (first column), 0.9 (second column), 0.95 (third column), 0.99 (fourth column). Panels on the right show the circumferential stress ${\tilde{T}}_{\mathit{\Phi }\mathit{\Phi }}$ at ${\mathcal{D}}_{\max }=0.99$ for the three considered values of ξ—this can be compared with the healthy profile in figure 3h.

Figure 6.

Key experimental metrics, apex displacement (AD; a) and central corneal thickness (CCT; b), plotted for varying damage parameters, ${\mathcal{D}}_{\max }$ and ξ.

Figure 7.

Comparison of reference, healthy and diseased (using ${\mathcal{D}}_{\max }=0.99$ and ξ = 4 in (4.2)) curvatures. Panels (a,b) show curvatures of the posterior and anterior surface, respectively.

Figure 8.

Comparison of key metrics of diseased corneas (using ${\mathcal{D}}_{\max }=0.99$ and ξ = 4 in (4.2)) for varying γ, K⁽²⁾ and K⁽³⁾. In left, central and right panels, we fix at their baseline values γ, K⁽²⁾ and K⁽³⁾, respectively, and vary the remaining parameters. We plot the apex displacement (AD; a–c), the central corneal thickness (CCT; d–f), the maximum anterior curvature (MAC; g–i) and the maximum posterior curvature (MPC; j–k). For comparison, we note that the reference (unloaded cornea) values are AD=0 mm, CCT=0.62 mm, $\mathrm{MAC}=1/{\tilde{R}}_A\approx 0.128\hspace{0.17em}{\mathrm{mm}}^{-1}$, $\mathrm{MPC}=1/{\tilde{R}}_P\approx 0.139\hspace{0.17em}{\mathrm{mm}}^{-1}$ and the values for healthy cornea loaded with physiological IOP are AD=0.04 mm, CCT=0.617 mm, MAC = 0.131 mm⁻¹ and MPC = 0.143 mm⁻¹. Black circles at bottom-left corners of the panels indicate the case where all parameters are fixed at their baseline values (γ = 20, K⁽²⁾ = 0.001 and K⁽³⁾ = 0.07).