An Improved Monte Carlo Method for Quantitative Analysis of Transparency Degradation Caused by Corneal Edema

Abstract

Purpose: To address the impact of corneal edema on transparency, we proposed an improved Monte Carlo (MC) method that uses a cost-effective and easily applicable approach to extract fibril features. This method aims to enhance computational accuracy and elucidate microscopic mechanisms of corneal transparency.

Methods: Fresh ex vivo porcine corneas were immersed in deionized water to simulate various states of corneal edema. For each edema state, we carefully analyzed its transparency and fibrous conditions. The diameter and spatial distribution of corneal fibrils were measured using image recognition algorithms. An improved MC method, analyzing the multiple scattering of each incident light, was proposed to simulate the propagation of incident light in the cornea.

Results: Higher edema reduces transparency. After 10 hours of hydration, corneal water content increased by 7% and transparency decreased by over 30%. Cross-sectional observations indicated that the radius increases from the original 15.1 nm to 16.8 nm, the fibril-to-base material area ratio decreases from the initial 15% to roughly 7%, and the minimum gap between fibrils expanded by 35%. The proposed MC method provides more accurate transparency estimates compared to the direct summation of field theory.

Conclusions: Corneal transparency decreases as a result of multiple factors, including fibril and cornea thickening, variations in the refractive index, and changes in the fibril distribution.

Translational Relevance: The MC method is cost-effective, easily applicable, and accurate in predicting the transparency of corneas with varying water contents. Moreover, this method intuitively demonstrates the scattering of light during the propagation process.

Introduction

Cornea possesses excellent mechanical properties and optical transparency. As the outer layer of the eye, it protects internal structures while allowing most visible light to pass through. When the pump function of the corneal endothelium is impaired,¹ aqueous humor will flow uncontrollably into the stroma, causing the cornea to swell. Corneal edema leads to increased corneal thickness, whitening, and turbidity. This is a common finding of corneal injury,^2,3 inflammation,^4,5 keratoconus,^6–8 and damage to the corneal endothelial layer.^9–11 Studying the micro-physiological structure of corneal transparency is essential for elucidating how cornea generates and maintains its optical and biomechanical functions. These studies provide a reliable basis for biological modeling, enabling clinicians in grasping pathologies and predicting the onset of diseases.^12,13 Additionally, they also support the development of artificial cornea materials.¹⁴ 

The cornea is constituted of stacked fibrils and extracellular matrix, creating a complex structure rather than a uniform material. There is a refractive index difference between the fibrils and the extracellular matrix.^12,15 The diameter of corneal fibrils (about 30 nm) is much smaller than the wavelength of visible light (390–760 nm).¹⁶ The fibrils in adjacent layers are orthogonal to each other, as well as parallel to the corneal surface,¹⁷ which is essential for transparency and mechanical properties of corneas. Fibrils can continuously scatter the incident light, leading to those light travels through the cornea in a path filled with obstacles. Since the mid-20th century, researchers have attempted to explain how the cornea maintains such high transparency. Maurice^18,19 suggested that the superposition of scattered fields would propagate posteriorly through the cornea by arranging corneal fibrils in a regular hexagonal lattice. The assumption of the regular lattice is not well consistent with observation results by electron microscopy studies of corneal fibrils. Thereby a subsequent proposition of a paracrystalline model was proposed for fibril distribution.^20–22 The direct summation of field (DSF) method assumes the scattering from a single fibril has a negligible effect on the incident wave and ignores the wave's attenuation.²³ This allows the use of the Born approximation to solve the Schrödinger equation, which describes the motion of the incident wave.²⁴ Scholars have predicted corneal transparency using the DSF theory.²⁵ The prediction results were not well consistent with the experimental results at the positions away from the center of the cornea. This may be because the DSF theory simplifies calculations by neglecting the effects of multiple scattering. In the practical scenario of light transmission, the fibrils will inevitably be excited by the scattered fields of other fibrils. This will produce additional vibrations and induce multiple scattering.¹⁴ When the position of fibrils is far away from the center of the cornea, this effect is strong.²⁵ 

Considering the effects of multiple scattering, the MC method, also known as the statistical simulation method, establishes a probabilistic model that matches the parameters or numerical characteristics of the problem and gives approximate solutions. This method has been widely used in fields such as computational physics and biomedicine.^26–29 However, in practice, extracting appropriate and accurate physical parameters is still challenging. This work proposes a low-cost and easily applicable approach for fibril parameter extraction and an improved MC method for accurate prediction of corneal transparency. 

Methods

Sample Preparation

Twenty-four porcine eyeballs were obtained from a local slaughterhouse in Guangzhou for this study. All eyes were intact and undamaged. They were stored in physiological saline solution at 4°C and delivered to the laboratory within four hours postmortem. Corneal tissues were dissected from the eyeballs using a scalpel, and the corneal epithelium was scraped off with a blade. A total of 16 corneas were placed in deionized water at 4°C to induce stromal edema. Among them, eight of the corneas were hydrated for five hours, whereas the other eight were hydrated for 10 hours. The remaining eight normal corneas, served as the control group, were stored in physiological saline solution at 4°C. 

As shown in Figure 1, preparing samples for transmission electron microscopy (TEM) involves six stages: fixation, dehydration, embedding, curing, sectioning, and staining. The goal of fixation is to preserve the microstructure of the cornea and prevent deformation, as well as degradation. The fixation was carried out at 4°C. First, the corneas were fixed in a 2.5% glutaraldehyde solution for four hours, followed by a secondary fixation in 1% osmium acid solution for two hours. After each fixation step, the corneas were rinsed with 1% phosphate buffered saline solution for 15 minutes. Second, the corneas were dehydrated in a gradient with ethanol solutions at concentrations of 30%, 50%, 70%, 90%, 95%, and 100%, respectively. The 100% concentration ethanol dehydration was done twice. The duration of each dehydration was 20 minutes. And then a treatment in pure acetone for 20 minutes was done at the end of the dehydration phase. Third, embedding was conducted at room temperature. The corneas were infiltrated with acetone and embedding agent (Epon812 resin) mixtures at volume ratios of 2:1 for four hours, followed by a ratio 1:2 for another four hours. They were then immersed in pure embedding agent for eight hours. Fourth, during the curing stage, place the samples in ovens at 35°C, 45°C, and 60°C for 12 hours, 12 hours, and 48 hours, respectively. After curing, the samples were sectioned into 70 nm thick slices using an ultramicrotome (EM KMR3; Leica, Wetzlar, Germany), ensuring that the slices were thin enough to allow electron beam penetration. The sections were then stained with 2% uranyl acetate and lead citrate solutions for 15 minutes each to enhance image contrast and resolution. Finally, the sections were observed using a transmission electron microscope (HT7800; Hitachi, Tokyo, Japan). 

Figure 1.

 

Six stages of sample preparation for TEM: (a) fixation; (b) dehydration; (c) embedding; (d) curing; (e) sectioning; (f) staining. And (g) physical photo of TEM.

Transparency Measurement

Transparency is determined by measuring the loss of light intensity after it passes through the sample. An optical lens transmittance tester (LS108H; Linshang Technology, Shenzhen, China) was used to test the paths of visible light wavelengths of 400 nm, 500 nm, and 600 nm. The tester has a beam diameter of 1 mm and measurement resolution of 0.1%. Place the fresh and edematous corneas into the optical lens transmittance tester, ensuring the geometric center of the cornea to align with the optical center of the tester along the same vertical line, consistent with the location of TEM. The laser is directed perpendicularly onto the center of the cornea surface to minimize the impact of corneal curvature. Once the optical path stabilizes, the ratio of the transmitted light intensity to the incident light intensity indicates the transparency of the cornea. 

Transparency Calculation

The transparency of the cornea primarily depends on the ordered arrangement of fibrils. The theoretical calculation of transparency hinges on the determination of the scattering cross-section of the scatterer.^18,20,24 The scattering cross-section is defined as the ratio of scattered light intensity to incident intensity in a certain scattering angle, which quantifies the directional scattering probability. The value of the scattering cross-section is influenced by factors such as the wavelength of the incident light, the size of the scattering particle, and the refractive index of the scattering particle.^20,30 

Natural light typically enters the cornea perpendicularly. As shown in Figure 2(a), the corneal fibril layers oriented parallel to the corneal surface. This alignment allows us to simplify the scattering of incident light using a two-dimensional approach,¹⁴ shown in Figure 2(b). Thus the differential scattering cross-section can be briefly described as  

\begin{eqnarray}d\sigma = \frac{{dn}}{{NdS}}\end{eqnarray}

(1)

where N is the total number of incident particles, dS is the area of target screen, and dn is the number of particles scattered into the target screen. When a corneal fibril behaves as a single scattering particle, its scattering cross-section can be given by²⁴ 

\begin{eqnarray}\sigma (\theta ) = \frac{{n_{\rm{s}}^3{{{(\pi a)}}^4}{{{({{m}^2} - 1)}}^2}}}{{2{{\lambda }^3}}}\left[1 + {{\left(\frac{{2\cos \theta }}{{{{m}^2} + 1}}\right)}^2}\right]\end{eqnarray}

(2)

where θ is the plane scattering angle, n_s is the refractive index of the corneal stroma, a is the fibril diameter, m is the ratio of the refractive index of the fibril to that of the stroma, and λ is the wavelength of the incident wave. 

Figure 2.

 

(a) The cornea is the first gateway for natural light to enter the human eye. Light enters the eyeball perpendicular to the corneal surface and corneal fibrils. (b) The scattering path of incident light by a single fibril. The scattering cross-section is described as the ratio of the light intensity within the scattering angle to the incident light intensity, reflecting the probability of scattering in different directions.

During a complete scattering process, the attenuation of the light intensity is given by  

\begin{eqnarray}dI = - I\rho {{\sigma }_t}d\delta \end{eqnarray}

(3)

where I is the intensity of the transmitted light, dδ is the thickness of medium, ρ is the density of scattering particles, σ_t is the integral sum of the differential scattering cross-section of all fibrils, also known as the total scattering cross-section. After passing through the cornea with a total thickness of δ, Equation 3 integrates in the domain of integration I₀ (the intensity of the incident light) through I. Thus the intensity of the transmitted light is given by  

\begin{eqnarray}I = {{I}_0}\exp \, ( - \rho {{\sigma }_t}\delta )\end{eqnarray}

(4)

 

The transparency is defined as the ratio of transmitted light intensity to the incident light intensity:  

\begin{eqnarray}T = \frac{I}{{{{I}_0}}} = \exp \, ( - \rho {{\sigma }_t}\delta )\end{eqnarray}

(5)

 

It is known that I is directly proportional to the square of the amplitude of the scattered wave defined as ψ. For a single fibril with a plane wave incidence, we can describe it using the time-independent stationary Schrödinger equation:  

\begin{eqnarray}{{\nabla }^2}\psi + ({{k}^2} + U(r)) \, \psi = 0\end{eqnarray}

(6)

where U(r) is the potential function and k is the wave number. Its solution is a linear superposition of the attenuated incident wave and the scattered spherical waves:  

\begin{eqnarray}\psi = A \exp \, (ikz) + \frac{{f(\theta )}}{r}\exp \, (ikr)\end{eqnarray}

(7)

where A is the attenuation coefficient of the incident wave, and f(θ) is the scattering amplitude. When the incident light passes through multiple-layered stacking corneal fibrils, (Equation 7) will change to be:  

\begin{eqnarray}\psi &=& {{A}_n}{{A}_{n - 1}} \cdot \cdot \cdot {{A}_2}{{A}_1}\exp (ikz) \nonumber \\ && + \sum\limits_j^{{{N}_n}} {\sum\limits_j^{{{N}_{n - 1}}} { \cdot \cdot \cdot \sum\limits_j^{{{N}_2}} {\sum\limits_j^{{{N}_1}} {\frac{{{{f}_j}(\theta )}}{{{{r}_j}}}\exp (ik{{r}_j})} } } } \end{eqnarray}

(8)

where A_n is the attenuation coefficient of the incident light passing through the n^th layer, and N_n is the number of fibrils in the n^th layer. In reality, the intensity of the incident wave attenuates rapidly,^14,18 and the fibrils at different positions receive an incident wave that is a complex superposition of multiple scattered waves from the fibrils in front of them. Therefore A_n and f_j(θ) are obtained with difficulty; as a result, the theoretical calculation of the multiple scattering wave is extremely challenging. 

MC Method for Multiple Scattering

To address this challenge, the MC method uses an approach tracking the movement of incident particles. The first and most important step is to obtain the radius and distribution characteristics of the fibrils from the TEM images. At the center of the cornea, the fibrils are arranged orthogonally along the superior-inferior (S-I) and nasal-temporal (N-T) directions, respectively. As shown in Figure 3a, TEM imaging is done at the center of the cross-section in the S-I direction to clearly present the fibrils. And then, we developed an image recognition program to accurately detect the thickness and position of the fibrils. As shown in Figures 3b and 3c, the images are binarized to enhance the contrast between the fibrils and the extracellular matrix. Next, the best-fitting circle for each fibril bundle's contour is identified using edge detection methods. Finally, the radius and their centers, representing the corresponding fibril's radius and position, are recorded. As shown in Figure 3d, the radial distribution function (RDF), a statistical method, is used to describe the spatial distribution of particles. The RDF is the ratio between the local density of fibrils around a specific fibril and the total density, which effectively illustrates fibril aggregation patterns. As distance increases to long range, the value of the RDF approaches 1, reflecting that the distribution of corneal fibrils becomes uniform. This characteristic provides the statistical foundation that enables the MC method to simulate multiple scattering in corneal fibrils. 

Figure 3.

 

The process of obtaining the radius and the RDF of corneal fibrils. (a) Use trephine to obtain cornea from fresh porcine eyeball. Cut cornea open alone S-I direction to expose the position of TEM. (b) Corneal fibrils within a layer. (c) Detecting corneal fibrils through the most approximate circles to obtain the fibril radius and center position. (d) The RDF describes the ratio of local fibril density to total fibril density as a function of distance from a fibril center. Calculate the average radial distribution function to reflect the characteristics of fibril distribution.

The second step is to identify key parameters in the statistical model. These parameters involve the scattering cross-section of the fibrils and the free path of the incident particles. The free path refers to the average distance a particle travels before encountering another fibril after scattering, which depends on the fibril distribution and falls between the minimum and maximum distances between fibrils. From Equation 2, we can calculate the scattering cross-section. 

The final step is to prepare a large number of incident particles (more than 10⁶) with defined initial positions and wavelengths. After each scattering event, the position and energy of the incident particles will be updated. This updated process will continue until reaching the stopping criteria, such as detected by a receiver, exiting the boundary or reaching a maximum number of scatterings. The activation probability of these criteria is affected by the fibril distribution. As shown in Figure 4a, the normal distribution of corneal fibrils helps keep the scattering angle range and the free path of incident light reasonable. Most of the light can reach the receiving area behind the cornea. As shown in Figure 4b, when the free path increases because of sparse fibril distribution, incident light is more likely to pass through the boundary. When the scattering angle increases, the light will collide more with the fibril and lose more energy, as shown in Figure 4c. Because the intensity of the light wave is proportional to the number of photons, according to Equation 5, the ratio of particles received by the detector to the number of incident particles indicates the transparency of the cornea. 

Figure 4.

 

Simulating the motion of different incident light in the MC method. (a) The appropriate scattering angle and free path facilitate the forward traversal of incident light through the corneal fibrils. (b) The wider spacing between fibrils increases the free path, allowing incident light to exit the boundary more easily. (c) Increased scattering angles boost the light's scattering frequency from the fibrils. In the simulation, incident light exceeding the maximum scattering event stops moving due to energy depletion.

Results

The parameters used in the MC method are listed in the Table. The mass and the central thickness of corneas were measured in our laboratory. After 10 hours of hydration, corneal water content increased by 7%, and both corneal mass and thickness increase significantly. The hydration level was defined as the ratio of the water weight to the dry weight of the cornea. The fibril volume fraction was calculated as the ratio of the total fibril area to the area of the detection region. The fibril radius increases from the original 15.1 nm to 16.8 nm, the fibril-to-base material area ratio decreases from the initial 15% to roughly 7%. These parameters are the mean values of each group with eight corneas, respectively. As shown in Figure 5, the corneal fibrils after hydration are illustrated through TEM. From these images, the radius and volume fraction of corneal fibrils were determined using our developed image recognition program. There are about 10,000 fibrils detected in each group. The radius of the fibrils increases slightly, yet the fibril volume fraction decreases noticeably. 

Table.

 

Material Parameters Utilized in the MC Method of Multiple Scattering

 

Figure 5.

 

TEM Images of corneal fibrils under the same scale. (a) normal cornea; (b) hydrated for 5 hours; (c) hydrated for 10 hours. The images within the same group correspond to the upper, middle, and lower layers of the cornea, respectively. As the hydration time increases, it can be qualitatively observed that the fibril distribution transitions from dense to sparse, and neat to messy.

As shown in Figures 6a and 6b, the distribution of corneal fibrils at various edema levels can be described visually by the RDF method. The value of zero means there are no other fibrils within that distance from a given fibril. When the RDF moves from zero to non-zero, the corresponding distance represents the nearest distance. When the RDF moves to the highest value, where the local fibril density is much higher than the average, the corresponding distance represents the densest distance. Both the nearest and densest distances reflect the characteristic of the short-range order of corneal fibrils. As the distance increases from the short-range zone into the long-range zone, the RDF value approaches 1. The corresponding fibril distribution becomes uniform, which is essential for maintaining the cornea's stable biological properties. As shown in Figure 6c, with increasing corneal edema, the swelling of the stroma leads to greater nearest and densest distances between fibrils. The short-range effects of fibril distribution become more pronounced, and compared to the more uniform long-range distribution, the short-range distribution appears more disordered at a larger scale. 

Figure 6.

 

RDF of corneal fibrils: (a) normal cornea; (b) edematous cornea. (c) The relationship between the short-range characteristics of corneal fibril distribution and hydration. The densest distance is the area with the highest local density of fibrils, where a fibril has the highest probability of discovering other fibrils. The nearest distance represents the minimum gap between fibrils, with no other fibrils present within this distance.

Figure 7 shows the results of experimental measurements and transparency calculations. These calculations were used the DSF and MC method, respectively. In a non-edematous cornea, the calculation results of the DSF and MC method slightly deviate from the experimental data. As corneal edema increases, the calculation results of the DSF method become larger, whereas those of the MC method remain consistent with the experimental data. Meanwhile, the cornea also becomes hazy and whitish, which is well in agreement with clinical observations. It is noted that these corneas exhibit lower transmittance for lower wavelength light. 

Figure 7.

 

Comparison of the DSF theoretical values, MC theoretical values, and experimental measurements: (a) normal cornea; (b) hydrated for 5 hours; (c) hydrated for 10 hours.

Discussion

As mentioned above, the calculation results of the MC method are closer to the experimental data rather than those calculated by the DSF method. When applying DSF,^20,24,25,30 researchers first assess the scattering effect of a single fibril, then combine the effects of all fibrils. The DSF method is simple and intuitive but overlooks the continuous attenuation of incident light and multiple scattering among fibrils. Our improved MC method tracks incident particles as they undergo complex scattering among the fibrils, simulating realistic light intensity attenuation during transmission. Using the developed image recognition program, we accurately capture the size and distribution of corneal fibrils. The obtained fibril radius in porcine corneas is well agreed with the results measured by X-ray scattering.³¹ The existing research suggests the main causes of corneal transparency loss due to edema are (1) increases in corneal thickness and fibril radius; (2) greater differences in refractive indices between fibrils and the extracellular matrix; and (3) increased short-range disorder in fibril distribution, leading to a more uneven arrangement.^14,25 Figure 8 illustrates microscopic mechanisms of the main causes obtained from the improved MC method. After edema occurs, the radius of the corneal fibrils as scatterers slightly increases. Because the scattering ability of the fibrils grows with the fourth power of their radius, their scattering capacity is enhanced significantly. Additionally, as the cornea hydrates, its thickness increases, leading to a longer path for transmitted light and less transparency. This phenomenon can be described by Equation 5. Both fibrils and the extracellular matrix absorb water. According to observation of TEM images, the volume increase of the extracellular matrix is much greater than that of the fibrils, indicating that more water is absorbed by the extracellular matrix. The difference in water-binding capacities disrupts the already existing refractive index imbalance, making light transmission less efficient, because hydration reduces the refractive indices of both the fibrils and the extracellular matrix.¹⁵ The short-range arrangement of corneal fibrils is crucial for enhancing light transmission.¹² Over a hydration period of 0 to 10 hours, the short-range zone of the fibril distribution expands significantly, over 30%, increasing the free path of scattered particles. Incident particles are more likely to escape the boundaries rather than being received. 

Figure 8.

 

Mechanism of corneal transparency deterioration due to corneal edema. (a) Normal water content cornea has high transparency. (b) Most of the incident light can pass through the normal cornea. (c) The edematous cornea will swell and deteriorate in transparency. (d) Corneal edema leads to more light intensity loss. (e) Microscopic mechanisms of light intensity loss due to edema.

The common DSF theory uses the Born approximation to solve the Schrödinger equation for scattering, which is only accurate when the scattering is a perturbation relative to the incident light. When the wavelength of the incident light is short, according to Equation 2, the scattering cross-section of fibrils increases sharply. Moreover, when highly edematous cornea occurs, scattering cross-section of fibrils also increases. These will induce that Born approximation is not satisfied, leading to a large deviation between the calculation result and the experimental result. Compared to the DSF theory, our improved MC method has merits: (1) tracking the trajectories of incident light particles rather than using the Born approximation; (2) considering multiple scattering among fibrils. As shown in Figure 9, prediction errors of the MC method are much lower than those of the DSF method, especially in the conditions of high hydration and short wavelengths. The predicted errors of the MC method are less than 4%, even in the conditions of low hydration and longer wavelengths. 

Figure 9.

 

Relative errors of theoretical values compared with experimental measurements under different working conditions: (a) the DSF method; (b) the MC method. The arrow line shows the trend of the value.

In summary, this article presents an improved MC method for predicting corneal transparency. This method tracks the trajectories of incident light particles and considers multiple scattering among fibrils. It can elucidate the microscopic mechanisms underlying the reduced transparency of edematous corneas and well predict the corneal transparency with errors less than 4%. We also developed an image recognition program to accurately detect the characteristic of the fibrils. These results help understanding of scattering abnormalities in edematous corneas and effective supporting advancements in biomedical imaging. 

Acknowledgments

Supported by National Natural Science Foundation of China (Grant Nos. 12432008, 11932007, 12472180, 12472179, 12372181, 12072116, and 12072115), Guangdong Provincial Natural Science Foundation (Grant Nos. 2024A1515011076 and 2023A1515012942), and State Key Laboratory of Subtropical Building and Urban Science (Grant No. 2022KA05). 

Disclosure: S. Li, None; Y. Liu, None; B. Yang, None; S. Lv, None; X. Zou, None; L. Tang, None; L. Zhou, None; Z. Jiang, None; Z. Liu, None