Simulation of time-fractional oxygen diffusion in cornea coated by contact-lens.

In this paper, the time-fractional oxygen diffusion has been simulated in a one-dimensional (1D) corneal-contact lens (CL) system. Different CLs have been employed as Balafilcon, thin- and thick-Polymacon. It is assumed that homogeneous and isotropic porous mediums of cornea and CL is saturated with compressible oxygen. The computations of the time-fractional derivations are done based on the Caputo method. The obtained results show that the fractional derivative order (FDO) severely affects pressure distribution in cornea and CL. Consequently, the magnitudes of post-lens-tear-film (PoLTF) pressure change due to diverse FDOs. Particularly, maximum changes have been observed in the results gained from the CLs with thicknesses more than 100 μm. The agreement of the results obtained from the time-fractional modeling with the experimental data compared to the standard diffusion modeling has been improved by more than 36%. Finally, it has been demonstrated that high-thickness CLs can cause exist anomalous diffusion process in cornea tissue.

Introduction

It is generally accepted that the standard topology of fluid particles diffused in a space has a linear relationship to time (Zayernouri 2015). However, this point of mathematical view is sometimes unable to find suitable responses for experimental tasks. For the sake of this incapacitation, it is preferred on some issues to take into account the anomalous diffusion process. In this regard, the displacements of fluid particles can have a non-linear relationship to time. From the point of mathematical view, first-order derivations of time replace with fractional-time derivations. The application of the time-fractional calculus is nowadays increasing in computational analysis of diseases behavior such as COVID-19 transmission patterns (Awais et al. 2020), epidemic childhood diseases (Baleanu et al. 2020), damages in skin tissue caused by laser irradiation (Hobiny et al. 2020), growth of tumor cells in the presence of chemotherapy (Kumar et al. 2021) and Alzheimer medical phenomena (Mohammad and Trounev 2021). Fractional modeling of these medical problems allows finding appropriate curves fitted by the unpredicted distribution of experimental data. In this paper, the fractional modeling concept is extended to achieve a better understanding of oxygen diffusion in cornea tissue coated by CLs.

Numerous studies have been accomplished related to the estimation of corneal oxygen diffusion. Fatt et al. (1974) studied pressure distribution for a corneal-CL system assuming steady-state conditions. The CLs included Softlens, Bionite, and Silicon rubber lenses. Harvitt and Bonanno (1999) considered the effect of acidosis on pressure distribution in cornea worn by the polymethyl methacrylate, Oxyflow f30, Fluorex 700, and Oxyflow 151 lenses. Brennan (2005) improved Harvitt and Bonanno’s solution to eliminate negative pressure detections. In these three reviewed studies, the simplified governing equations have been solved employing the solution algorithms of ordinary differential equations. Alvord et al. (2007) simulated oxygen diffusion in a two-dimensional cornea geometry with average measures using the finite-element method. It is reported that the actual shape of any individual cornea can significantly have a deviation from the used dimensions. Although there is a layered pressure distribution same as 1D geometry in most cornea volumes, the peripheral part of the cornea contains spatial variation, especially in the closed-eye situation. Chhabra et al. (2009b) determined oxygen diffusion numerically for the cornea-CL system based on the second Fick’s law. To this end, the reaction term is defined by the Monod kinetics model. Different values of oxygen solubility coefficient and oxygen permeability have been computed for a certain cornea when coated by the Balafilcon and thin-Polymacon CLs. Larrea and Büchler (2009) offered a new numerical simulation for transient diffusion of oxygen in light of the Michaelis-Menten-type relationship. Chhabra et al. (2009b) have developed a more complete model accounting for the role of reactive species such as glucose, bicarbonate, lactate, hydrogen, and carbon dioxide along with oxygen to define an index as oxygen deficiency factor in oxygenation of the cornea-CL system. Takatori and Radke (2012) prepared a quasi-two-dimensional model for respiration of the cornea to gauge corneal hypoxia. Del Castillo et al. (2015) repeated the methodology suggested by Chhabra et al., assuming constant cornea characteristics. It has been claimed that the maximum oxygen consumption rate (OCR) is unfixed when the oxygen pressure diminishes at the interface CL. Compañ et al. (2017) and Moreno et al. (2021) have been analyzed the application of the generalized Monod kinetics model to describe the human corneal OCR during various soft CLs wear. It has been observed that maximum OCR increases to oxygen pressure of 105 mmHg and then decreases. These variations can be related to limitations in all the OCR models. Therefore, a generalization of them should be performed to acquire a better description of the behavior of the cornea in humans. The common denominator of all the recent studies reviewed is that the standard form of Fick’s law has been used for almost all simulations performed. Besides, the corneal-CL system has been taken into account as a 1D geometry. The cornea is thin compared to its area, so the effects of lateral diffusion are neglected (Harvitt and Bonanno 1999).

From the points of material and skeleton views, a certain cornea tissue must have constant characteristics such as oxygen solubility coefficient, oxygen diffusion coefficient, and oxygen permeability. The employment of different CLs can solely affect the value of the maximum OCR. This idea is utilized in the numerical simulations by Del Castillo et al. (2015). However, the final outputs differed from the laboratory results by Bonanno et al. (2002). To this end, anomalous diffusion simulations are presented again, considering the constant characteristics of corneal oxygen. In other words, the main focus of the current paper is on the development of the oxygen diffusion analysis in cornea based on the time-fractional derivative. Most of the fractional-time derivatives are modeled using the Riemann-Liouville (Fei et al. 2021; Wu et al. 2021; Zhou et al. 2021) and Caputo (Caputo 1967; Almeida 2017; Wu et al. 2020; Liu et al. 2021a,b) approaches. Unfortunately, the Riemann-Liouville approach leads to initial conditions containing the limit values of the Riemann-Liouville fractional derivatives at the lower terminal (Podlubny 1998). The main advantage of Caputo’s approach is that the initial conditions for fractional differential equations with the Caputo derivatives take on the same form as for integer-order differential equations (Podlubny 1998). Therefore, the fractional-time derivative is modeled here employing the Caputo-based derivative. It is deserved to say that the current study explores the effect of various FDOs on the pressure distribution and PoLTF pressure of CL-equipped cornea as a novel experience. The obtained results can be used in the geometrical design of CLs. The presented simulations include a 2-layer system made of the cornea, , and CL, , mediums. The CLs are in the types of Balafilcon, thin- and thick-Polymacon materials. Mathematical modeling of the problem, as well as the results and discussions of time-fractional simulations, is described in the following sections with more details.

Mathematical modeling and assumptions

Figure 1 depicts a corneal-CL system. The cornea geometry with a thickness of L, which its interior surface contains endothelial cells. Epithelial cells of the cornea are coated by an isothermal post lens tear film. The CL geometry has a thickness of L whose exterior surface is exposed to atmospheric air pressure at the open-eye situation and palpebral conjunctive at the closed-eye situation. Since the CL is porous, oxygen particles can penetrate and reach the cornea. Practical assumptions given below are made to simplify the present study:

Fig. 1

Schematic representation of a 1D cornea-CL system

The human cornea consists of three different porous mediums as Bowman’s layer, corneal Storma, and Descemet’s membrane (Mobaraki et al. 2019), which transport metabolic species (oxygen, lactate ion, glucose, hydrogen ion, bicarbonate ion, and carbon dioxide) (Chhabra et al. 2009b). Porous mediums of cornea and CLs are homogeneous, isotropic, and saturated with oxygen.

Characteristics of cornea and CLs are invariable.

Cornea and CL mediums have chemical equilibrium.

Cornea thickness is much less than its surface area, and thus, lateral diffusions have been neglected.

The spatial variation of oxygen pressure is 1D.

Thicknesses of endothelium, epithelium, and stroma are integrated into cornea thickness.

Thicknesses of pre- and post-lens tear films (∼3 μm) are ignored.

Oxygen is considered a compressible fluid.

Standard governing equations

The standard governing equations are presented here to facilitate the creation of the time-fractional equations represented in Sect. 2.2. In an open-eye situation, the governing equation describing anomalous oxygen diffusion in the cornea-CL system in Fig. 1 consists of the 1D and time-fractional Fick’s law expressed as (Chhabra et al. 2009a) where represents the distribution of oxygen pressure at time in -direction. and () are oxygen solubility coefficient and oxygen permeability for each medium, respectively. denotes the OCR, which can be computed based on the Monod kinetics model as shown below (Chhabra et al. 2009a). where is maximum OCR. denotes the Monod dissociation equilibrium constant. Initial conditions for the governing equation can be obtained from a closed-eye situation using the steady-state solution of Eq. (1) as (Chhabra et al. 2009a) Both time-dependent and steady-state governing equations can be solved under boundary conditions of the Dirichlet type summarized in Table 1 on the basis of the pressures of the anterior chamber, ac, air, air, and palpebral conjunctiva, pc. In addition, the employed characteristics of the cornea medium are represented in Table 2 based on the Del Castillo et al. (2015) data.

Table 1

Information of boundary conditions (Chhabra et al. 2009a)

Pressure symbol	Value*
Pac	24
Ppc	61.5
Pair	155

*All values are in mmHg

Table 2

Characteristics of cornea medium (Del Castillo et al. 2015)

kΩc×10−5 (mL/cm3/mmHg)	(Dk)Ωc (Barrer*)	Lc (μm)
3.1	93	532

*Barrer = 10−11 (cm2/s)(mL O2(STp)/cm3/mmHg)

To complete the required boundary conditions for each of open- and closed-eyes situations, two more boundary conditions are defined at the location of post lens tear film based on continuity of pressure distribution and oxygen flux continuity as (Chhabra et al. 2009a):

Time-fractional modelling

Time-fractional modeling of the diffusion process in the cornea-CL system can be represented as shown below. where is the order of fractional derivative, and / is the fractional (non-integer) time derivative. The standard oxygen diffusion equation obtains considering = 1. Therefore, Eq. (6) is the same as Eq. (1) when = 1. Anomalous diffusion process can be observed with 0 < < 1. Using the theory of operators, the time derivative term of Eq. (6) can be rewritten as below: where denotes fractional derivative operator.

Discretization method

For the numerical solution of the problem, the spaces of the cornea and CL have been decomposed with and points, respectively. The time scale is partitioned with levels. Discretized form the time-fractional term presented in Eq. (7) can be computed on the basis of the Caputo methodology as (Caputo and Cametti 2008): where is spatial point indicator. denotes the Gamma function. The left-hand side of Eq. (8) can be computed numerically in the following manner (Toubaei et al. 2019): where is the time level indicator. refers to time step, and Discretizing the diffusion term in Eq. (6) based on the central finite-difference approximation as well as the utilization of the result of Eq. (8), the implicit discretized time-fractional equation at th time level can be presented as where is the spatial step. The non-linear partial differential equations of Eqs. (3) and (11) can be solved by the Newton-Raphson method (Press et al. 1989).

Numerical results and discussion

The simulation results of the time-fractional oxygen diffusion in a cornea coated by distinct CLs are first studied in this section and then discussed compared to standard oxygen transport. The utilized CLs are in the types of Balafilcon, thin- and thick-Polymacon; their characteristics are been given in Table 3. It deserves to be mentioned that the lens characteristics are taken from the Chhabra et al. (2009b) data.

Table 3

Characteristics of the CLs medium (Chhabra et al. 2009a)

CL	Qc,max×10−4 (mL(STp) cm−3 s−1)	K_m (mmHg)	k_Ωl×10^−5 (mL/cm3/mmHg)	(Dk)_Ωl (Barrer)	Ll (μm)
Balafilcon	1.4	2.2	6.7	100.5	100
Thin-Polymacon	0.9	2.2	1.5	90	60
Thick-Polymacon	0.7	2.2	1.5	90	200

The influence of various FDOs, {0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.}, on oxygen pressure distributions is shown in Fig. 2 after 60, 40, and 70 respectively inside the cornea separately coated by Balafilcon, thin- and thick-Polymacon CLs. As can be seen, the FDO value has a significant influence on the pressure distribution within the cornea-CL medium. Moreover, the usual linear distribution of pressure in the standard oxygen diffusion process, =1 within the CL medium, has been converted to a non-linear behavior in the anomalous oxygen diffusion by adjusting different FDOs less than 1. The effect of the same FDOs on the post-lens tear film pressure is illustrated in Fig. 3 for the cornea while wearing different CLs. Again, it can be seen that different values of the FDOs strongly influence the process of the post-lens tear film pressure changes.

Fig. 2

The influence of the FDO on the oxygen pressure distribution in the cornea coated by different CLs

Fig. 3

Effect of the FDO on the PoLTF pressure via various CLs wearing

To select the most appropriate amount of the FDO, the value of the relative error norm (REN) is calculated in terms of experimental, exp, and calculated, cal values based on the following form: where NOE refers to the number of those experimental data reported by Bonanno et al. (2002). Figure 4 shows the value of the REN versus the FDOs for each of the CLs. It can be seen that the REN value for the Balafilcon and thick-Polymacon CLs has been obtained less than 10% at = 0.9. The common denominator of both CLs belongs to their thicknesses which are more than 100 μm. These CLs can therefore reveal that the selection of high thicknesses for CLs can suggest the anomalous oxygen diffusion in the cornea medium. It can also be seen in the same figure that the REN of the thin-Polymacon CL in = 1 is less than 10%. The thickness of this CL is 60 μm. Therefore, it can be claimed that a CL with thickness less than 100 μm can present standard oxygen diffusion in the corneal medium.

Fig. 4

Computation of REN for different FDOs

Based on the appropriate values of the FODs calculated in Fig. 4, the oxygen transfer simulation has been executed again, considering time-fractional modeling in the cornea covered by each of the CLs. To this end, steady-state simulations have been first accomplished using the iterative Newton-Raphson method to compute initial conditions. The convergence diagram of these simulations is shown in Fig. 5. To achieve accurate initial conditions, the stop criterion is set to 10−15. This figure demonstrates that the initial pressure distributions have been obtained via fast iterative processes, finally up to 7 iterations per CL.

Fig. 5

Convergence history of steady-state simulations

The results of the PoLTF pressure distribution are plotted in Fig. 6. The experimental Bonanno et al. (2002) data and the standard diffusion results of Del Castillo et al. (2015) are also given in this figure. It is obvious for the Balafilcon and thick-Polymacon CLs that the results of the time-fractional modeling rather than the standard diffusion results agree better with the experimental data. The improvement of the PoLTF pressure calculation obtained from the time-fractional modeling against the published numerical data of standard diffusion is presented in Table 4 based on the value of REN. It can be deduced that the time-fractional method can improve the modeling of PoLTF pressure distribution by more than 36% for the CLs with thicknesses greater than 100 μm.

Fig. 6

Comparison of the PoLTF pressure obtained from the time-fractional modeling with the published experimental and numerical data

Table 4

Improvement of the PoLTF pressure computation using the time-fractional modeling against the standard modeling proposed by Del Castillo et al. (2015)

CL name	Del Castillo et al.	Time-fractional modeling	Improvement (%)
Balafilcon	4.6945	2.9771	36.58
Thin-Polymacon	4.9749	4.9749	NOCHG*
Thick-Polymacon	15.0190	7.6832	48.84

*NOCHG: No Change

Concluding remarks

In this paper, the process of anomalous oxygen diffusion has been simulated in a certain cornea coated by different CLs. The Balafilcon, thin- and thick-Polymacon CLs have been selected to investigate the results. It is assumed that compressible oxygen penetrates in 1D homogeneous and isotropic porous mediums of cornea and CL without any chemical reaction. The anomalous diffusion process has been modeled by the fractional-time method with various FDOs. To integrate non-integer time derivation, the Caputo approach based on an implicit discretization scheme is selected. The Newton-Raphson algorithm has been employed to solve non-linear discretized equations.

The FDO has a significant effect on pressure distribution in a corneal-CL system. Linear pressure variation in the CL medium can be converted to a non-linear one. Moreover, minimum pressure magnitude inside the cornea changes in the fractional-time modeling. Similarly, PoLTF pressure manipulates based on the time-fractional modeling. Particularly, all significant changes are evident for the CLs with thicknesses greater than 100 μm. The agreement of the results obtained from time-fractional modeling with the experimental data demonstrates that high-thickness CLs lead to form anomalous diffusion process in cornea tissue.