A mathematical model for the corneal transparency problem

Understanding the physical basis of corneal transparency has been a subject of interest amongst physicists, basic scientists and ophthalmologists. Impairment of corneal clarity is a significant cause of visual morbidity worldwide. Thus, it is essential to understand the mechanisms behind corneal transparency and how the alterations due to corneal pathologies affect vision. We use Maxwell’s equations to model light propagation in ocular tissues and a nodal discontinuous Galerkin method combined with an explicit Runge-Kutta method to simulate light propagation in normal and pathological corneas. Our simulation results illustrate that an increase in the diameter of some fibres causes an increase in backscattering. Thus, these may represent some of the physical changes in the cornea that might result in loss of transparency and visual morbidity.


Introduction
The eye is a slightly asymmetrical globe composed of various tissues (see Fig. 1), about an inch in diameter, which possesses a quite complex and sensitive mechanism.
The cornea is the transparent, clear tissue at the front of the eye. It is a differentiated tissue to allow refraction and transmission of light. Its characteristics are highly specialised in maintaining transparency and resisting external adverse factors like dehydration, microbial invasion, and trauma. Optically, it is a convex-concave lens with an average refractive power of 43 diopters, which corresponds to about 70% of the total refractive power of the optical ocular system [8].
The cornea comprises six layers: the epithelium, basement membrane, Bowman's layer, the stroma, Descemet's membrane, and the endothelium (see Fig. 1). Many corneal pathologies are caused by changes in at least one of these layers, leading to increased light scattering and consequent loss of corneal transparency. The stroma is the layer that gives the eye essential corneal strength and constitutes 90% of the cornea's thickness [3,22].
Corneal transparency is achieved and maintained by combining several features: the specific structural organisation of the collagen fibres in the stroma, the absence of blood vessels in the corneal layers, the lack of myelin sheath in the corneal nerves and the ability to preserve a tightly controlled hydration state. These factors contribute to the cancellation of light scattering in all directions except forward and to the low absorption in the spectral range from 300 nm to 1400 nm (UV-A, visible and IR-A light).
The corneal stroma consists of collagens, keratocytes and other extracellular matrix components (EMC). The collagen molecules are organised into fibrils 1 with uniform diameters between 25 nm and 35 nm, further gathered into collagen lamellae. There are from 200 to 300 lamellae of various thicknesses in the stroma lying parallel and interwoven together. Lamellae are quite interwoven in the anterior region, whereas parts near Descemet's membrane seems to be less interwoven, allowing the cornea there to swell easier. Getting closer to the sclera, the collagen lamellae get disorganised and therefore less transparent [22].
Explaining which factors lead to corneal transparency deterioration, the respective extent and where the responsible changes occur within the cornea, are still open questions that this works aims to address and shed light on.
The paper is organised as follows. Section 2 describes the mechanisms that explain corneal transparency and the causes of light scattering in corneal pathologies. In Sect. 3, we propose a computational model that mimics light propagation through the eye's tissues. We use Maxwell's equations to model electromagnetic radiation and, as discretisation methods, we use a nodal discontinuous Galerkin method for space discretisation and a low-storage explicit Runge-Kutta method for time discretisation. Section 4 presents the results obtained for the computational simulation of light propagation in a normal and pathological corneas. Section 5 discusses the results.

Corneal transparency
Corneal transparency depends strongly on the structural organisation of the stromal collagen fibres, namely on the uniformity of their diameter fibres [15] and on the regularity of their spatial arrangement. Many different ideas have emerged on this topic over time [4,5,9,10,12,13,15,19,[21][22][23]. In particular, to ensure corneal transparency, it is sufficient that the distance between adjacent fibres is constrained. This structural organisation ensures that while light is scattered by the collagen fibres, the scattered light interferes destructively in all directions except the forward direction. This way, the cornea acts as a waveguide, allowing light to travel forward through the cornea to reach the retina [25]. Doutch et al., in 2008, summarised parameters that are believed to be crucial in maintaining corneal transparency [7]. These are: (i) density of collagen fibres; (ii) diameter of collagen fibres; (iii) refractive index differences between fibres and interfibrillar space or EMC; (iv) thickness of stroma, and (v) order in the spacing of the fibre network. Thus, every model trying to explain corneal transparency should consider the shape and size of corneal stroma and its elements such as collagen fibres and proteoglycans, and their refractive indexes because these factors influence the amount and direction of scattered light. The average refractive index of the stroma has been calculated to be 1.375. It is based on the refractive index of collagen fibres and interfibrillar fluid and on the volume of collagen fibres and fluid [23].
Another important requirement for maintaining corneal transparency is that collagen fibres do not touch each other. This is achieved due to the negative charge of the proteoglycan-rich coating around collagen fibres, which prevents their aggregation [5]. The space between the collagen fibres depends on the water content in the matrix. The transparency is maintained by the balance between the proteoglycan matrix and the collagen fibres [17]. The physiological hydration in the cornea is tightly maintained at about 78% in humans through the pump and barrier functions of the endothelium. When the cornea swells, the distance between collagen fibres within the lamellae increases, and voids form within the lamellae [4], increasing light scattering and decreasing the cornea transparency.
The uniformity of the diameters of the collagen fibres and the restriction of the range of distances between the collagen fibres are two of the most important factors for corneal transparency. Several corneal pathologies that increase light scatter and result in corneal opacity are associated with phenomena that affect those two factors.

Maxwell's equations
We use Maxwell's equations to model the electromagnetic wave's propagation through the eye's tissue [16,20]. Maxwell's equations are the fundamental set of equations that describe how the electromagnetic field propagates in free space and in any media. Four-vector fields describe the electromagnetic field in space and time: E, H, D, B : × R + → R 3 , where E represents the electric field, H the magnetic field, D the electric flux density and B the magnetic flux density. The SI units of these fields are volt per meter (V /m), ampere per meter (A/m), coulomb per square meter (C/m 2 ) and weber per square meter (Wb/m 2 ), respectively. We consider that is a bounded domain in R 3 .
Assuming that there is no electric current density nor electric charge density and using the constitutive relations D = εE, B = μH, where ε is the medium's electric permittivity and μ is the medium's magnetic permeability, the Maxwell's equations can be written as: for heterogeneous anisotropic media with no source in three-dimensional spaces.
We decompose the electromagnetic wave in transverse electric (TE) and transverse magnetic (TM) modes, reducing the number of equations implemented in our model. TE or TM modes describe the propagation when the electric field lays in the plane of propagation or when perpendicular to it, respectively. Those two modes can be decoupled since they don't contain any common field vector components. These assumptions are appropriate when studying light propagation in two-dimensional spaces. From now on, we will analyse the time-domain Maxwell's equations in TE mode where the only non-vanishing components of the electromagnetic fields are E x , E y and H z . Thus, we have the following equations: where E = (E x , E y ) and H = (H z ). These equations are solved on the bounded domain ⊂ R 2 . We use the following notation for the vector and scalar curl operators: We consider that the electric permittivity ε is an anisotropic tensor, and that the magnetic permeability μ is isotropic. Both, μ and ε, are functions in space. We complete the model equations with the initial conditions E(x, y, 0) = 0 and H(x, y, 0) = 0 in and with first-order Silver-Müller absorbing boundary conditions [1,14] ν is the speed with which a wave travels along the direction of the outward pointing unit normal ν. These boundary conditions are relevant since we need to consider a truncated computational domain to mimic an unbounded propagation space. In this case, when the electromagnetic wave hits the boundary, it should not be reflected but absorbed.

Discretisation methods
There has been a great interest in solving Maxwell's equations accurately and efficiently in realistic applications in the last decades because of their relevance in many different areas. The first -and most well-known -method for solving Maxwell's equations is the so-called finite difference time-domain (FDTD) scheme proposed by Yee in 1966 [27]. Despite its success, FDTD, like all finite difference methods, is difficult to generalise to irregular domains and unstructured grids.
With the growing need to solve geometrically complex large-scale problems, there has been an interest in the flexibility offered by finite volumes or finite element schemes [16]. More recently, the development of discontinuous Galerkin (DG) time-domain methods to solve Maxwell's equations gained relevance on the simulation of electromagnetic waves propagations. DG methods gather many desirable features such as being able to achieve high-order accuracy and easily handle complex geometries. The advantages of using DG methods when compared with classical FDTD methods, finite volume time-domain methods or finite element time-domain methods, have been reported by several authors (see e.g. [11] and the references therein). Our computational algorithm is based on a nodal DG method for space discretisation and an explicit Runge-Kutta method for time discretisation.
We start by the conservation form of the TE-mode Maxwell's equations in 2D equivalent to (2a)-(2b), with

Space discretisation
Let the domain be a bounded set with a polygonal boundary ∂ and T h = {T k } K k=1 be a conformal triangulation of . On each element T k , the solution fields are approximated by polynomials of degree less or equal to N . The global solution u(x, y, t) is assumed to be approximated by piecewise N order polynomials.
defined as the direct sum of the K local polynomial solutions u k h . The fields are expanded in terms of interpolating Lagrange polynomials L i (x, y), where (x i , y i ), i = 1, . . . , N p , states for interpolation points whose number N p is related with the polynomial order N by N p = (N + 1)(N + 2)/2. We start by multiplying equation (4) by test functions v (which are usually the Lagrange polynomials) and integrate over each element T k . Next, we employ one integration by parts and substitute in the resulting contour integral the flux F by a numerical flux F * . Integrating by parts again yields where ν = (ν x , ν y ) is the outward pointing unit normal vector of the contour.

Introducing the notation for jumps of fields values [u
where "+" refers to the neighbouring element and "-" refers to the local cell, we define the following upwind flux where Z ± = μ ± /ε eff ± is the local impedance and Y ± = (Z ± ) -1 is the local conductance. The effective permittivity is defined by ε eff = det ε/ν T εν.
For Silver-Müller absorbing boundary conditions, we consider at the outer boundary which is equivalent to set

Time discretisation
The discretisation in space leads us to a global semi-discrete scheme, which can be written as A low storage explicit Runge-Kutta (LSERK) method is applied to these equations to solve the problem. One of the most celebrated LSERK methods is a fourth-order version with five stages, LSERK(5, 4), introduced by Carpenter & Kennedy in [6]. When compared to the classical fourth-order explicit Runge-Kutta methods, the application of LSERK(5, 4) diminishes the memory requirements related to arrays storage without increasing the computational cost. One alternative to LSERK (5,4) is the version that uses 14 stages instead of five, proposed by Niegemann et al. in [24], the LSERK (14,4). Unlike LSERK (5,4), this method comes with an extra cost of ten additional functions due to the ten more stages. However, LSERK(14, 4) allows a large stable time-step [2].

Numerical experiments
We consider a two-dimensional model of backscattered light intensity in two different (healthy and pathological) scenarios that correspond to the organisation of the fibres represented in Fig. 2. Our simulation settings are composed of 38 randomly distributed collagen fibres denoted by circles in the mentioned figure. Figure 2 (left) represents a healthy human cornea with fibres of diameter 31 nm and where the distance between every two adjacent fibres is not less than its doubled diameter [22]. In Fig. 2 (right), we mimic a pathologic situation where positions of fibres are kept and 20% of those (8 fibres) are randomly chosen to have doubled diameter [22]. On the computational domain = [-1, 1] 2 , we represent by F the union of circles that model healthy corneal collagen fibres' positions. \ F represents the EMC stromal components. The union of circles F represents a pathologic cornea where the diameter of some of the fibres is doubled.
Since the eye's tissues are non-magnetic media, the relative magnetic permeability is μ = 1 [26]. The refractive indices n of collagen fibres and EMC are well-established in the literature. Taking those values from [18] we define In our experiments, we consider ε = I, where = n 2 and I is the identity matrix. This yields that the relative permittivity of the medium is given by 0 = 1.365 2 .
To simulate the scattering of light in the cornea, we use the scattered field formulation of the TE-mode. We begin by separating the electromagnetic fields u = (E x , E y , H z ) T into two fields: the incident fields u i = (E i   Assuming that the incident field is also a solution of Maxwell's equations (4), with coefficients 0 and μ 0 = μ being the relative permittivity and permeability of the medium in which the incident field propagates in the absence of scatterers, we insert (6) into (4) and we obtain the scattered field formulation We want to solve for the scattered fields with initial conditions u s (x, y, 0) = 0 and the incident fields given by the planar wave u i (x, y, t) = (0, cos (10(xt)), 0) T . For the discretisation, we use the DG method with the LSERK (14,4) with the spatial meshes illustrated in Fig. 3. The evolution in time of the scattered electric field intensity Fig. 4, where the solution is approximated by a polynomial of order N = 4 and plotted for different values of simulation time T. The results in Fig. 4 show an increase in backscattering caused by the increase in the diameter of some fibres, which is compatible with previous findings (e.g. [22]).
All simulations were done using the MATLAB software package on a computer with one 11th Gen Intel Core i7-1165G7 @ 2.80 GHz × 8 processor and 8 GB of RAM. The meshes used in the paper were designed using the software package FreeFem++ and the computational time ranged from 26 seconds, for T = 0.25, and 105 seconds, for T = 1, in both scenarios.

Conclusion
This study assessed the impact of a small percentage of stromal fibres being modified on their properties, specifically, the increase of their diameter and the consequences of that Time evolution of scattered electric field intensity I s . Left: healthy cornea; Right: cornea with some fibres whose diameter is doubled change in corneal transparency. We, therefore, simulated a corneal pathology resorting to computational methods to find that an increase in backscattering occurs being compatible with previous findings in [22]. We may also note that the planar characteristic of the wavefront is significantly lost in the situation where we double the diameter of 20% stromal fibres. Our results, depicted in Fig. 4 for these conditions, demonstrate the increase in the background scattering of light that results from the aggregation of adjacent fibres and the resulting loss of its transparency. Our mathematical model can be used to simulate other pathological scenarios, for instance, considering other spatial arrangements in the increased fibre diameters or allowing regions devoid of fibres ( [22]).