Modeling the efficacy profiles of UV-light activated corneal collagen crosslinking

Objective Analysis of the crosslink time, depth and efficacy profiles of UV-light-activated corneal collagen crosslinking (CXL). Methods A modeling system described by a coupled dynamic equations are numerically solved and analytic formulas are derived for the crosslinking time (T*) and crosslinking depth (z*). The z-dependence of the CXL efficacy is numerically produced to show the factors characterizing the profiles. Results Optimal crosslink depth (z*) and maximal CXL efficacy (Ceff) have opposite trend with respective to the UV light intensity and RF concentration, where z* is a decreasing function of the riboflavin concentration (C0). In comparison, Ceff is an increasing function of C0 and the UV exposure time (for a fixed UV dose), but it is a decreasing function of the UV light intensity. CXL efficacy is a nonlinear increasing function of [C0/I0]-0.5 and more accurate than that of the linear theory of Bunsen Roscoe law. Depending on the UV exposure time and depth, the optimal intensity ranges from 3 to 30 mW/cm2 for maximal CXL efficacy. For steady state (with long exposure time), low intensity always achieves high efficacy than that of high intensity, when same dose is applied on the cornea. Conclusions The crosslinking depth (z*) and the crosslinking time (T*) have nonlinear dependence on the UV light dose and the efficacy of corneal collagen crosslinking should be characterized by both z* and the efficacy profiles. A nonlinear scaling law is needed for more accurate protocol.


Introduction
Photochemical kinetics of corneal collagen crosslinking (CXL) and the biomechanical properties of corneal tissue after CXL have been extensively explored in both animal and human models [1][2][3][4][5]. The dynamics of CXL and the safety and efficacy issues of CXL have been explored theoretically [6][7][8][9][10][11][12][13] and clinically [14,15]. To increase the speed of CXL procedures, accelerated CXL using high UV power (9 to 45 mW/cm 2 ) are also reported [16]. However, the accelerated CXL based on the linear theory of Bunsen Roscoe law (BRL) [17] to shorten the irradiation time is still a debating issue. More recent clinical studies have indicated that the efficacy based on BRL is actually lower than the non-BRL [18,19]. Prior modeling work of Schumacher et al. The study in [6] assumes a constant riboflavin concentration during the UV light exposure time. This assumption ignoring the dynamic of the riboflavin concentration leads to major errors in the calculations of the UV light intensity profile (both spatial and temporal), the rate of photopolymerization and the increase of stiffness. The CXL efficacy is characterized by multiple factors including the concentration and diffusion profile of the RF in the stroma, the absorption constant of the photolysis products and the quantum yield of the CXL process.
The critical formulas to be developed in this study include the crosslinking time, the crosslinking rate and the efficacy defined by the conversion of monomer-polymer characterized by the following parameters: the three extinction coefficients, concentration and diffusion depth of the riboflavin solution, the UV light intensity and irradiation duration, and the corneal thickness.
This study will present, for the first time, the nonlinear law for the CXL efficacy, in contrast to the conventional Bunsen Roscoe law [17].

The modeling system
The corneal model system [11,12] consists of its epithelial layer and the underlying stroma collagen as shown in Fig 1A. The UV light is normally incident to the corneal surface. The theory developed in this study can apply to both epithelium removed (epi-off) and epithelium intact (epi-on) case by slight revisions. However, we will focus on the more efficient epi-off case as shown in Fig 1B, where the riboflavin (RF) solution has an initial surface concentration (at z = 0) C 0 and distribution in the stroma defined by C 0 (z) = C 0 F(z, D), with F(z) = 1-0.5z/D having a diffusion depth D defined by half-width of the half-maximum, or C 0 (z = D) = 0.5 C 0 . The distribution function, F(z), is chosen to obtain analytic formulas for the roles of the diffusion depth D. One may also use a more realistic distribution based on measured data which, however, will require numerical simulations and many of the important features to be found analytically will not be available. We also assume that the UV light has a uniform intensity distribution across the stroma surface.

The dynamic intensity
In the modeling system shown in Fig 1B, the UV light intensity I(z, t) and RF concentration C (z,t) in the corneal stroma (for the epi-off case) may be described by a set of coupled firstorder differential equations, or by the integral equations [8][9][10][11] as follows, Where a = 83.6λϕε 1 , with ϕ being the quantum yield and λ being the UV light wavelength; ε 1 and ε 2 are the extinction coefficients of RF and the photolysis product, respectively. I 0 is the initial UV light surface intensity, or I(z = t = 0) = I 0 .; and C 0 is the initial RF surface concentration assuming a distribution profile given by F(z) = 1-0.5z/D, or C(z,t = 0) = C 0 F(z), with a diffusion depth D in the stroma. We note that the above described epi-off model may be easily applied to the epi-on case as follows. Given the epithelial thickness d, the intensity in Eq (1.b) is revised by redefining the corneal thickness as z' = z + d, that is, extra absorption of the epithelium reduces the intensity on the stroma surface (at z = 0) by a factor of exp(-Ad), where A is the absorption coefficient of the epithelium which may be slightly different from that of stroma.
The prior work of Schumacher et al [6] based on a non-depleted RF concentration, i.e., the assumption of aE = 0 in Eq (1.b), significantly overestimates the RF concentration for t>0. This assumption is also based on the continuing instillation of the RF drop during the CXL process such that the RF is always in its saturated state without depletion. That is, after the pretreatment time, they treat the CXL as a steady process without solving the dynamics of CXL. They also use an oversimplified model to assume no absorption from the photolysis products, i.e.,ε 2 = 0. Therefore, their calculated profiles, Figs 2 and 3, significantly deviate from our exact numerical profiles to be shown later.
The reported measurements [20,21] provide the parameters ε 1 = 204 (%Ácm) -1 and Q = 13.9 (cm -1 ), whereas ε 2 is not yet available in human, but is estimated to be in the range from 80 to 120 (%Ácm) -1 by the RF depletion test [11]. The quantum yield ϕ will be treated as a free parameter in our calculation and have a range from 0.3 to 0.5. The following units are used: C(z,t) in weight percent (%), I(z,t) in (mW/cm 2 ), λ in cm, Q in (cm) -1 and ε j (for j = 1,2) in (%Ácm) -1 . As shown by Eq (1), there are three major UV absorption components in the CXL The initial (at t = 0) RF concentration distribution inside the stroma given by a distribution function F(z,D), having a diffusion depth (D) [11,12].
process: the absorption of the stroma tissue (Q), which is independent to the RF concentration; the absorption of the unreacted RF solution (ε 1 C 0 ), and the photolysis product (ε 2 C 0 ), both are proportional to the initial RF concentration C 0. The initial UV light intensity (at t = 0) is obtained by the integration of Eq (1.a) with C(z,0) = C 0 F(z) as follows:  where (for j = 1, 2) where A 1 (for j = 1) is the initial state (at t = 0) absorption coefficient, independent to ε 2 ; A 2 (for j = 2) is the steady-state absorption coefficient, which is independent to ε 1 because of the complete concentration depletion, C(z, t = 1) = 0. Notably, G(z) in Eq (2.d) is the integration of F(z) over z.

The effective dose
The effective dose (or fluence) applied to the cornea collagen (at a depth z) for a UV exposure time (t) is defined by The UV light intensity is a dynamic function of time and z and requires a full numerical simulation of Eq (2). For comprehensive analytic formulas, various numerically fitting techniques are presented earlier [11]. In this study we will use the linear approximation of the light intensity I(z, t) given by the mean value of the initial (I 1 ) and steady intensity (I 2 ) defined by I (z, t) = 0.5[I 1 (z) + I 2 (z)] for t<T and I(z, t) = I 2 (z) + 0.5[I 1 (z)−I 2 (z)] for steady-state, or t>T, where T is the steady-state cross-linking time given by T = m/(aI 0 ), with m = 12 fit numerically with the exact solution of Eq (2) [12] Time integration of the light intensity, the effective dose (at z) Eq (3) becomes [12] EðzÞ where E 0 = I 0 t is the surface dose (at z = 0), and g is a correction factor for the transient state given by g = 0.5(I 2 -I 1 )(T/E 0 ), T = f/(aI 0 ), with the fit parameter f = 12. With ε 1 = 204 (%Ácm) -1 , ε 2 = 102 (%Ácm) -1 , Q = 13.9 (cm -1 ) and C 0 = 0.1%, we obtain an approximated value of the correction factor , which is about 10% to 25% correction to the steady-state when E 0 = 3.0 to 4.0 J/cm 2 and at z = 400 um.

The crosslinking time
Crosslinking time (T Ã ) may be defined in various ways. Basically, it is defined when the crosslinking procedure is mostly completed, or when the RF initial concentration is mostly depleted by the UV light and the procedure reaches a steady state having a very low polymerization rate. We choose to define T Ã by one of our most preferred definition based on the level of RF concentration depletion as follows.
Substituting Eqs (4) to (1.b), we obtain the analytic formula for the RF concentration as follows Eq (5) provides us the formula for the crosslinking time defined by C(z, t = T Ã )/C 0 F(z)(e −N ), which leads to the formula for the cross linking time given by T Ã ðzÞ ¼ T 0 ð1 þ ag=MÞexpðA 2 zÞ: ð6:aÞ Using the numerically fit expression of g, we obtain T Ã ðzÞ ¼ T 0 ½expðA 2 zÞ þ 1:5ð1 À expðÀ A 0 zÞÞ; ð6:bÞ where A' = (A 1 -A 2 ) and T 0 is the surface crosslinking time given by The efficacy profile where K is the ratio of the growth and termination rate constant of the polymer chain. Solving for Eq (8) we obtain the degree of the monomer-polymer conversion, which is proportional to the CXL efficacy (Ceff) as follows.
Rðz; t 0 Þdt 0 ð9:bÞ Analytic formulas have been discussed for the simplified case that the RF depletion is ignored, and for the case that absorption of UV light in the stroma (without RF) is ignored (or Q = 0) [22]. Analytic formulas including both Q and the RF depletion can be obtain by using the mean UV light intensity such that I(z,t) = 0.5 (I 1 + I 2 ) with I j defined by Eq (2), as follows.  (8) and (10) shows that photoinitiation rate (R) is a product of two competing factors, the RF concentration and the UV light intensity. Taking dS/dz = 0, we may find the optiomal z Ã for maximum S and Ceff, and z Ã is proportional to the dose, or exposure time for a given UV light intensity. Eq (11) also shows that the steady state efficacy (when atHI 0 >>1), is inverse proportional to the square root of [aI 0 ] which will be further demonstrated by our numerical results later. Defining a steady state by when P(z,t) < 0.133 and almost independent to time, we obtain the condition for a steady-state dose as follows: E 0 = tI o > 0.026 exp [H(z)], for E 0 in (mJ/cm 2 ).

Rðz
We may solve for Eqs (9) and (11) to obtain the exposure time formula for a given CXL efficacy as follows. Taylor expending of the ln term in Eq (12.a), we obtain a nonlinear scaling law of I 0 -0.5 Bunsen Roscoe law (BRL) [17]. Eq (12) also provide the formula for the crosslink time (T Ã ) which is nonlinearly proportional to X -0.5 , which may be compared to Eq (6) (13) is given by when E<0.133, or when the dose E 0 = tI o > 2.5 (mJ/cm 2 ). Which is similar to the condition defined by Eq (11.b). We shall note that there is no exact solution for Eqs (1) or (9) when Q or p is included. Therefore, we solve Eq (1) numerically by finite element method which is justified by comparing the numerical results wit the exact solution, Eq (13) for the special case of Q = p = 0 and D>>1. The roles of Q, p and D in Eq (13) will be fit numerically and presented elsewhere.

Results and discussions
Unless specified, the following calculated results are based on the measured parameters [17,20] of ε 1 = 204(%Ácm) -1 and Q = 13.9 cm -1 and the assumed quantum yield ϕ = 0.5 and ε 2 = 50 (%Ácm) -1 .  (11), that for the same UV dose (or energy), higher intensity achieves lower efficacy, and the efficacy is proportional to the square root of [C 0 /I o ]. The efficacy steady state value are analytically available from Eqs (9) and (11) with P = 1. Fig 5 shows the optimal intensity of CXL efficacy at depth of 250 and 400 um, where the saturation and optimal features, with optimal intensity (I Ã ) ranging 3 to 12 mW/cm 2 (for z = 250 um), 3 to 12 mW/cm 2 (for z = 400 um). These features are also predicted by Eq (11), and they are inverse proportional to the exposure time.

The nonlinear law
As predicted by Eqs (9) and (11), and our numerical data, the CXL efficacy at transient state (for small dose) is proportional to [t 2 I 0 ] 0.5 and at steady-state, it is a nonlinear increasing function of [C 0 /I 0 ] 0.5 or [t/E 0 ] 0.5 . This nonlinear scaling law is clinically more accurate than that of the linear theory of Bunsen Roscoe law (BRL) [17] and may be used as the new guidance for the protocol of accelerated CXL in relating the exposure time and the UV light intensity. The conventional accelerated CXL protocol based on BRL, therefore, has undervalued the exposure time (t) for higher intensity using the linear scaling of t = [E 0 /I 0 ], rather than t = [E 0 /I 0 ] 0.5 , if based on our nonlinear law. To achieve the same CXL efficacy higher RF concentration requires higher UV light intensity; and for the same dose, higher UV light intensity requires a longer exposure time.
The recent clinical data [23] showed that for the same dose, 9 W/cm 2 and 10 minutes is more efficient than 30 W/cm 2 and 3 minutes (based on BRL). This clinical new finding is consistent with our theory (as shown in Figs 2 to 4). However, the theory presented by Kling and Hafezi [23] based on a simplified model same as that of Schumacher [18]  shows thin cornea has higher efficacy than thick cornea, which is opposite to our predicted feature shown in our Figs 2 to 4 that CXL efficacy always lower at small z (or corneal depth), except at the transient state which shows an optimal z Ã . The simplified model [18,23] ignoring the RF depletion, or assuming a constant C(z,t) in our Eqs (1.a) and (1.b), totally eliminate the important competing feature between the UV light intensity and the RF concentration. Therefore their efficacy profiles are incorrect, although Kling and Hafezi [23] claimed that their clinical data (shown by their Fig 4) are consistent with their modeling. We believe that a wrong modeling can still fit to correct clinical data when there is free adjusting parameters in the modeling.

The demarcation line and crosslink depth
Gatzioufas et al. [24] and Moramarco et al. [25] reported the CXL efficacy and its relation to the stromal demarcation line depth (SDD). The mean SDD were reported at 149 um and 213 um, respectively, in continuous wave (CW) and pulsed light-assisted CXL. Therefore pulsed light-assisted CXL could be more clinically effective than standard CW light. These findings also support the importance of oxygenation for effective CXL, as suggested by Richoz et al. [26].
However, Gatzioufas et al. [24] stated the open debating issues that whether the SDD was a true indicator of CXL efficacy; "the deeper, the better" principle may apply to CXL or other factors interfered with the clinical interpretation of the SDD after CXL. It was reported that the high UV light intensity achieved more superficial the stromal demarcation line [18]. However, there was no proportional decrease in CXL efficacy with increasing UVA light intensity [19].
The above controversial issues may be interpreted by our theory as follow. As shown by Figs 2 and 3, the optimal crosslinking depth (z Ã ) is an increasing function of E 0 , D and ϕ, but it is a decreasing function of C 0 . In comparison, the CXL efficacy, Eqs (9) and (11), is an increasing function of C 0 and the UV exposure time (for a fixed UV intensity, as shown by Fig 3), but it is a decreasing function of UV light intensity (at a fixed dose as shown by Fig 4). As expected, higher concentration results in a stronger absorption of the UV light and reduces the crosslinking depth (z Ã ). Therefore, the CXL efficacy should be characterized by both z Ã and Ceff, which have opposite trend with respective to the UV light intensity and RF concentration.
Our nonlinear scaling laws demonstrate that, for the same dose, higher UVA light intensity achieves lower Ceff than that of lower intensity. It is known that the elasticity of the anterior corneal stroma is significantly higher than that of the posterior corneal stroma. Therefore, thicker corneas require a larger z Ã (or longer exposure time) than that of thin corneas to achieve the same CXL efficacy. As suggested by Gatzioufas et al. [24], the minimal (or optimal) amount of dose (energy) and shorter treatment time (using higher intensity) might further improve the CXL safety and reduce its potential complications.

Corneal thinning effects
It was reported that corneal thinning occurs during intra-operatively during the UV exposure. The thinning effect may be analyzed by the role of the diffusion depth (D) defined in Eq (2.d) as follows. The thinner cornea results a larger effective riboflavin diffusion depth (D) in the stroma and therefore increases the crosslink depth (z Ã ) which can be estimated by Eqs (2) and (10). For example, a typical thinning of 100 um and D increases from 300 um to 400 um, our calculations show a larger z Ã . However, the currently available clinical data are not sufficient for our modeling to include the actual influence of the thinning effects on the CXL efficacy. Further investigations are required.
Our one-dimensional modeling assumes a uniform distribution of riboflavin across the corneal surface (horizontal direction on the x-y plane), but focuses on the distribution along the depth (z) is included. A two-dimensional model is in progress to include the horizontal features which may be compared to topography data with variation on the x-y plan.
Maximal efficacy CXL efficacy may be described by either the measured demarcation line depth or by the increase in stiffness of the crosslinked collagen. Maximal CXL efficacy and crosslink depth would require optimization of the UV light intensity (or dose) for a given riboflavin concentration. As shown in Fig 5, the optimal intensity ranges from 10 to 30 mW/cm 2 for the maximal CXL efficacy. To improve the CXL efficacy, various techniques have been proposed. These include: pulsed mode operation of the UV light, extra oxygen supply to the corneal surface; and enhancement of the riboflavin diffusion such as diffusion in the de-epithelialized stroma (standard method); diffusion through the epithelium into the stroma (transepithelial method); or direct introduction of riboflavin into the stroma (pocket technique, ring technique, needle technique); and enrichment of riboflavin in the stroma by iontophoresis. In addition to UVA light activated riboflavin, other photosensitizers using blue light (at 430 nm) and green light (at 532 nm) were also proposed.
Detailed photochemical kinetics of CXL was reported showing the roles of oxygen in both type-I and type-II reactions [3,22]. In addition, combination of supplemental oxygen and pulsing UV exposure providing substantial improvement of CXL efficacy was also reported by Muller et al. (Avedro Inc internal Report). Epi-off technique has been known to be much more efficient than the epi-on which has much less oxygen in the stroma.
We should note that our nonlinear scaling law and the related stromal demarcation line depth require further clinical studies to customize its value under various CXL conditions. These factors include the concentration and diffusion profile of the RF in the stroma, the absorption constant of the photolysis products and the quantum yield of the CXL process, besides the operation modes (CW or pulsed) of the UV light which may influence the oxygen environment in a type-II process.