Fractional Beer-Lambert law in laser heating of biological tissue

1.   Introduction

Laser heating is a common thermal therapy in clinical practice. It is known that the use of a laser as a heating source is employed in the cochlea stimulation ^[32], the treatment of different cancer tumours ^[7,11,17,19], laser-assisted drug permeation ^[9], endovenous ablation ^[30], skin remodeling ^[25] and other applications ^[2]. This paper focuses on laser heating of the cornea and one of its applications: laser thermal keratoplasty (LTK). LTK enables the alteration of the intrastromal collagen in a precise and controlled manner, and, consequently, it can be employed to reduce hyperopia by several diopters.

In addition to experimental procedures, theoretical modeling is a usual methodology to assess the suitability of laser heating as a thermal therapy (see ^[6,12,13]). Accordingly, theoretical models have also been developed to achieve a deeper understanding and a broader knowledge of laser heating of the cornea, and, moreover, to make this technique more effective and confidently applicable ^[1,20,21,26]. However, to consider theoretical modeling as a reliable tool, the results from theoretical models should be in agreement with experimental ones. As the authors stated in ^[20], some factors may affect the accuracy of theoretical models. Therefore, in the last decades, many efforts have been made to improve mathematical models from different viewpoints. One of them corresponds to the change of the physics formulation.

The parabolic heat transfer equation (PHTE) regularly used to model laser heating is

where $T(\mathbf{x}, t)$ is the temperature at point $\mathbf{x}$ at time $t$, $k$ is the thermal conductivity, $\rho$ is the density, $c$ is the specific heat, and $S(\mathbf{x}, t)$ is the heat source in the material. In ^[18] we proposed an alternative formulation using a time-fractional derivative in the first term of equation (1.1), which we called the time-fractional derivative formulation of the heat transfer equation (FDHTE). Our hypothesis was that the FDHTE, which can assume a non-ideal heat conduction behavior, could provide more accurate theoretical results with a better agreement with experiment. The results in ^[18] corroborated the hypothesis for FDHTE in cases in which a heat source was not applied. However, we did not find significant differences between the solutions of PHTE and FDHTE when the source was applied. The reason relies on the fact that the value of the term $S(\mathbf{x}, t)$ in (1.1) is several orders higher than the rest of the terms during the heating phase (when the laser source is applied). This occurs for the laser application exposed in ^[18], but also for LTK and a lot of applications of thermal laser heating. Thus, in cases in which the source term is significantly greater than the others, the effect caused by changes in terms of time or spatial partial temperature derivatives is almost negligible.

In the present work, our interest is also to explore alternative formulations of the PHTE for laser heating in thermal therapies. Due to the important role of the heat source $S(\mathbf{x}, t)$ in equation (1.1), we assume that by introducing an alternative and contrasted formulation in this term, we can obtain a different temperature distribution that can provide more accurate results compared to those obtained experimentally.

The heat source for laser heating is usually obtained from the Beer-Lambert Law (BLL). According to ^[23], the BLL predicts a laser intensity assuming almost ideal conditions, such as homogeneous sample without scattering, or no changes in the absorption despite the interactions between particles. These assumptions are unrealistic in biological tissue. In the case of the cornea, ^[31] stated that its mechanical properties are not isotropic in a plane perpendicular to the surface, and a resulting mechanical anisotropy is a defining characteristic of the cornea. For this reason, some extensions of the BLL have been considered. In ^[15] authors also stated that BLL was not appropriate to model anomalous natural behaviours (for instance, anomalous transport that can take place in turbulent plasmas, random lasers and biological cells), and a more realistic approach to this phenomenon can often be provided by a Fractional Beer-Lambert Law (FBLL). In fact, this idea was also underlined in ^[4], in which the authors marked the use of FBLL as an instrument to transmit randomness of microscopic dynamics to the macroscopic scale. In ^[4] authors gave a generalization of the Beer-Lambert Law derived from a fractional Poisson process, and it was later extended in ^[28]. These results were employed in ^[10] to introduce the FBLL by considering an extension of the Srivastava-Owa operators. The FBLL is useful for providing a generalized diffusion model based on the fractional radiative transport equation, as shown in ^[15,16]. Regarding nonbiological tissues, some recent works have shown the FBLL to be effective for exhibiting sub-exponential effects in aluminum foil half-value thickness ^[27] and, in ^[5], for approximating light attenuation inside photosynthetic cultures in photobioreactors.

This study aims to assess the temperatures obtained from the PHTE using a heat source term based on the FBLL. As far as we know, no studies have evaluated the influence of the FBLL inside the PHTE. Moreover, we have applied the theoretical formulation to the laser heating of the cornea in LTK.

2.   Mathematical modeling of the laser heating

We consider a 1-dimensional problem in which the laser incidence is produced at $x = 0$, and the heating process is produced along the $x$-axis. Assuming constant thermal conductivity, density and specific heat in Eq (1.1), the heat conduction problem from the viewpoint of the PHTE is given by

where $\alpha$ is the thermal diffusivity ($\alpha = \frac{k}{\rho c}$), $T_0$ is the initial temperature, $T_a$ is the ambient temperature, and $h$ is the thermal convection constant at the interface of the tissue with the ambient temperature. The diffusivity of the material $\alpha$ measures its ability to conduct thermal energy relative to storing it. The parameter $\alpha$ represents the speed at which the heat is diffused through a material, and its unit is m$^2$/s. Therefore, if the material has a high value of diffusivity, it will respond quickly to thermal changes. Meanwhile in a material with a small $\alpha$ value, the heat is mainly absorbed, and it will provide a slower response to thermal changes. At the surface $x = 0$, the tissue is affected by free convection of the surrounding air, as can be observed in the last condition of (2.1). This condition corresponds to Newton's cooling law, so the thermal convection constant $h$ is the rate of heat transfer between a solid surface (in this case, the cornea) and a fluid (the surrounding air) per unit area and per unit of temperature difference (W/m$^2\cdot$ K).

To solve problem (2.1), we need to introduce the $S(x, t)$ expression in the governing equation. Specifically, the relationship between the heat source and the laser intensity distribution $I(x, t)$ is:

where $b$ is the absorption coefficient in m$^{-1}$ and $H(t)$ is the Heaviside function, defined by

The parameter $b$, which is also called the attenuation coefficient, quantifies the absorption of laser light per unit length. In this sense, $b$ represents the depth at which the laser light of a specific wavelength can penetrate in a material before it is absorbed.

At this step we considered two cases for obtaining $I(x, t)$: 1) using the BLL and 2) using the FBLL.

2.1.   Modeling the laser source using the Beer-Lambert Law

The intensity of the laser $I(x, t)$ is obtained using the BLL from the problem:

where $I_0$ is the laser intensity applied at point $x = 0$, and $R$ is the reflection coefficient (%), which accounts for the fraction of laser light that is reflected by the material.

It is easy to verify that the analytical solution of (2.4) is given by $I(x, t) = I_0(1-R)e^{-bx}.$ Therefore, according to Eq (2.2),

Now, substituting expression (2.5) in the governing expression of (2.1), this problem can be solved as in the reference [29, Section 3.2].

2.2.   Modeling the laser source using the fractional Beer-Lambert Law

The intensity of the laser $I(x, t)$ is obtained using the FBLL from the problem:

where we replace the ordinary derivative with the Caputo space-fractional derivative of real order $\beta$ defined by

and

where $\Gamma$ denotes the standard Gamma function.

It should be noted that $b$ represents in (2.6) a constant coefficient with unit m$^{-\beta}$ to make the Eq (2.6) dimensionally consistent. We observe that, according to the existing literature ^[3,32], this dimensional property could also be explicitly highlighted in the coefficient $b$ as $b^{\beta}$ or $b_{\beta}.$

In order to solve Eq (2.6), we will use the Laplace transform as the main tool. Given a function $f:\mathbb{R}_+ \to \mathbb{R}$, the Laplace transform of $f$ is defined by

for every $s\in \mathbb{C}$ such that $Re(s) > \omega$ for some $\omega \in \mathbb{R}.$

The Laplace transform of the fractional derivative of a function $f$ is given by

and

We also recall the definition of the Mittag-Leffler function ^[14]:

where $\Gamma$ denotes the standard Gamma function, and the following Laplace transform which involves the Mittag-Leffler function:

We proceed to solve the FBLL equation (2.6) by applying the Laplace transform, obtaining:

Using (2.7), we get:

and then

Consequently, the solution of the FBLL is given by

where we have used formula (2.9) with $\mu = 1$, $\rho = \beta$ and $\lambda = b$.

According to (2.2), the heat source is

Observe that the case $\beta = 1$ corresponds to the exponential solution of the classical BLL given by expression (2.5).

We point out that problem (2.6) for $1 < \beta \leq 2$ can be stated as follows:

where an extra initial condition is needed. Proceeding as before, and applying formula (2.8), the solution of (2.12) is given by:

We will now proceed to solve the problem (2.1) when considering the FBLL in the heat source, that is, $S(x, t)$ given by (2.11). In order to facilitate the resolution of the problem, we will make the following change of variables:

Thus, the resulting problem is:

In what follows, let $M = bI_0(1-R)$. Taking the Laplace transform with respect to $t$ and denoting $D(x, s): = \mathcal{L}_t[V(x, t)]$, we get

Applying now the Laplace transform with respect to $x$ to (2.16) and denoting $F(r, s): = \mathcal{L}_x[D(x, s)]$, we obtain:

Substituting the initial conditions (2.17) and (2.18) in (2.19) and isolating $F(r, s)$, we get

Using formulas $15$ and $16$ from ^[24] and the convolution theorem of Laplace transforms, we get that the inverse spatial Laplace transform of (2.20) is given by:

Taking into account that $\lim_{x\longrightarrow\infty} D(x, s) = 0$, we obtain the reduced expression

Applying now the inverse temporal Laplace transform in (2.22) using formulas $181,183$ and $184$ from ^[24] and the convolution theorem of Laplace transforms, we obtain:

where $Erfc(z)$ is the complementary error function. According to the changes made in (2.14), the solution of (2.1) for $0 < \beta\leq1$ is given by

The case $\beta = 1$ corresponds to the solution of (2.1) when using the classical BLL. Following a similar procedure, if $1 < \beta\leq2$, the function $V_{\beta}(x, t)$ is given by:

2.3.   Laser thermal keratoplasty

As we have mentioned in the introduction, we apply the expression (2.24) to a specific treatment of laser heating called laser thermal keratoplasty (LTK). LTK is a minimally invasive refractive surgery employed to reshape the cornea and correct hyperopia. A holmium laser is usually employed to focus the energy on multiple points around the cornea, forming a ring. The delivery of energy is transformed into heat in the tissue, and the heat modifies the curvature of the tissue, and then reshapes the cornea, as shown in Figure 1.

We consider the typical parameters of LTK using a holmium:YAG laser as in the previous study ^[29]: i.e., the initial laser intensity $I_0 = 5\times 10^{8}$ W/m$^2$, the reflectance coefficient $R = 2.4\%$ and the absorption coefficient of the cornea $b = 2000$ m$^{-1}$. As thermal characteristics of the cornea, we use those reported in the same study ^[29]: $\rho = 1060$ kg/m$^3$, $c = 3830$ J/kg$\cdot$K and $k = 0.556$ W/m$\cdot$K. The initial temperature of the cornea ($T_0$) is assumed to be $35$ $^\circ$C, the ambient temperature ($T_a$) is assumed to be $20^\circ$C, and the thermal convection constant is $h = 20$ W/m$^2\cdot$K. As we are interested in the heating phase, we consider a time for the laser application of $200$ $\mu$s. The cornea is approximately $540\; \mu$m, and for this reason we have made the analysis in points from $0$ to $500\; \mu$m. Where the corneal tissue finishes, the aqueous humor begins, and it is followed by the crystallinelens, the vitreous humour, etc. Most of these layers are much thicker than the cornea. Due to the difference of thickness between the cornea and the other tissues, we can assume that the boundary condition of $T_0$ (body temperature) can be modeled as a limit when $x$ tends to infinity. From a physical point of view, this can be also modeled as a finite bar setting as a boundary condition $T(x = x_{max}, t) = T_0$, but in this case $x_{max}$ should have a very high value.

From the analytical solutions obtained for the FBLL, we use Comsol Multiphysics to represent the solutions.

3.   Results and discussion

We are interested in the temperatures of the cornea during all the heating laser process to assess the differences between using BLL and using FBLL. In this sense, Figure 2 represents the temperature in a portion of the spatial domain $x\in[0,500]$ $\mu$m at $t = 50\; \mu$s and at the end of the heating ($t = 200\; \mu$s); and Figure 3 represents the temperature evolution at the points $x = 0$, $x = 200\; \mu$m and $x = 500\; \mu$m. Both figures are made using the FBLL to compute the heating source produced by the laser application, $\beta$ is the fractional derivative parameter, using expression (2.24). We choose $\beta$ values next to the unit since the solution for $\beta = 1$ corresponds to the case of the BLL source.

Figure 2 shows that when $\beta < 1$ in the FBLL ($\beta = 0.95$ and $\beta = 0.90$), the heat conduction is predicted to be slower than using the BLL; and when $\beta > 1$ in the FBLL ($\beta = 1.10$ and $\beta = 1.05$), the heat conduction is faster. In terms of the heat source, with $\beta$ values smaller than $1$, a smaller quantitatively laser heat source is obtained, which provokes slower heat conduction. On the contrary, with $\beta$ values greater than $1$, a higher value of the heat source is obtained, and faster heat conduction is produced. In Figure 2A), we see moderate differences in temperatures for $t = 50\; \mu$s. In fact, we find that the maximum differences ($3.3\; ^{\circ}$C, which represents approximately an $8\%$) are between $\beta = 1$ and $\beta = 1.10$ at $x = 500\; \mu$m. However, we observe in Figure 2B) that the differences increase for longer times. At $t = 200\; \mu$s differences between $\beta = 1$ and $\beta = 1.10$ at $x = 500\; \mu$m are $12.5\; ^{\circ}$C which represents approximately a $21\%$ difference. So, Figure 2 shows that slight variations in the $\beta$ value can produce great changes in temperatures.

Figure 3 shows that the temperature behavior is very similar at $x = 0$ for all $\beta$ values considered. This fact was also observed in Figure 4 of ^[18] using the FDHTE and in Figure 2 of ^[29] using the hyperbolic heat transfer equation (HHTE). However, Figures 3B) and 3C) show that differences between different $\beta$ values' curves grow when we move away from the point $x = 0$. Note that in ^[18] and ^[29], these differences were not observed during the heating phase at $x = 0$, nor in the rest of the points of the spatial domain. Specifically, Figure 4 of ^[29] showed identical temperature results for PHTE and HHTE during the heating phase and in all the spatial domain. In the previous studies ^[18] and ^[29], we changed the heat conduction models using the PHTE, the FDHTE, and the HHTE formulations, but all of them considered the laser heat source provided by the BLL. Due to the predominant effect of the laser heat source over the rest of the terms, the changes in the temporal derivatives assumed in HHTE and FDHTE formulations, provided negligible differences during the heating phase compared to the PHTE formulation. Nevertheless, we observed that using the FBLL, slight changes in the values of $\beta$ translate into significant temperature variations even during the heating phase. Figure 3C) shows that, assuming the FBLL instead of the BLL, we obtained differences of approximately $22\%$ between $\beta = 1$ and $\beta = 1.10$ at $t = 200\; \mu$m (Figure 3C).

Figures 2 and 3 corroborate our hypothesis, since changes in the formulation of the laser heat source produce significant variations in temperatures. In this sense, the FBLL offers a new approach for modeling the heat conduction of laser heating which is more flexible and could better fit the data in cases where the BLL approach does not fit the data well. This study assumes a new insight to improve the formulation of heat conduction problems, since it is focused on the changes in the heat source, instead of changes in the spatial and temporal derivatives involved in the PHTE formulation.

Consequently, our goal in this article was to explore a more flexible model that could fit data better than the classical approach using the BLL.

The expression (2.24) that we have applied to LTK, can also be applied to other laser heating problems mentioned in the introduction (^{[9,11,17,19,32]}, for example). In the case of laser heating of biological tissues involving perfused organs, only the blood perfusion term should be added in the PHTE according to the Pennes Bioheat Equation used in ^[22], for instance.

$\beta$ and $b$ are two independent parameters which arise in the model in different ways. $\beta$ is the fractional parameter considered in the FBLL, and $b$ corresponds to the absorption coefficient, which appears in the BLL problem formulation (2.4), in the FBLL problem formulation (2.6) and in the heat source (2.2). The parameter $b$ is known in most of the cases, but when it is unknown, it needs to be adjusted. The processes to adjust $\beta$ and $b$ values have to be independent. In the case of LTK presented in the study, first we would apply an optimization method to adjust the value of $b$. Specifically, we would use the experimental values obtained in $x = 0$, since at this point there is no influence of the $\beta$ value (as we can see in Figure 3A), but $b$ has an important effect. Secondly, once the value of $b$ is adjusted, we would choose the $\beta$ value which better matches with the experimental results in points different from $x = 0$. It is important to note that the FBLL provides an alternative expression to the solution obtained with the BLL, differing from the resulting solution by adjusting the value of the absorption coefficient parameter $b$. Actually, adjusting the value of $b$ does not solve the mismatch between experimental and theoretical results when dealing with non-ideal laser heating conditions. As an example, we have modified the $b$ value in our model to assess its influence on the results.

Figure 4A) shows the temperature evolutions at $x = 200\; \mu$m for $\beta = 1$ and $\beta = 1.10$ considering $b = 2000$ m$^{-\beta}$ and $b = 3300$ m$^{-\beta}$, respectively. The case $\beta = 1$ corresponds to the BLL solution, and $b = 2000$ m$^{-\beta}$ is the value we have used for the LTK modeling. So, the case $b = 2000$ m$^{-\beta}$ with $\beta = 1$ is our reference curve. The curves corresponding to $\beta = 1.10$ with $b = 2000$ m$^{-\beta}$ and the reference curve are the same that appear in Figure 3B) for those $\beta$ values. What is new in Figure 4A) is the curve for $b = 3300$ m$^{-\beta}$ and $\beta = 1$. As can be seen in Figure 4A) it is possible to obtain the same behavior by changing the value of $\beta$ or changing the value of $b$, since the curves for $b = 3300$ m$^{-\beta}$ and $\beta = 1$ and $b = 2000$ m$^{-\beta}$ and $\beta = 1.10$ are coincident, and both differ in the same way from the reference curve. However, in Figures 4B) and 4C) we show that this statement is false. If we consider the temperature evolution at the point $x = 200\; \mu$m (Figure 4B)) and the same pairs of $\beta$ and $b$ values as those considered in Figure 4A), we observe that in this case the curve for $b = 3300$ m$^{-\beta}$ and $\beta = 1$ does not coincide with the curve for $\beta = 1.10$ and $b = 2000$ m$^{-\beta}$. Indeed, the curve for $b = 3300$ m$^{-\beta}$ and $\beta = 1$ is above the other two curves, assuming much faster heat transfer than the others. Therefore, the result of changing the $\beta$ and $b$ values does not produce the same differences compared to the reference curve. Moreover, changing the value of $b$, we cannot achieve that theoretical results agree better with the experimental in cases where the BLL overestimates or underestimates the solution. The explanation is clearer if we observe Figure 4C), in which the temperature evolution at the point $x = 500\; \mu$m is shown considering the same value pairs for $\beta$ and $b$ as in Figures 4A) and 4B). As we observe in Figure 4C) again the curve for $b = 3300$ m$^{-\beta}$ and $\beta = 1$ and the curve for $\beta = 1.10$ and $b = 2000$ m$^{-\beta}$ do not coincide. However, the most important fact is that at this point the curve for $b = 3300$ m$^{-\beta}$ and $\beta = 1$ is under the other curves, assuming a slower heat transfer than in the other cases. Therefore, we cannot establish a correlation between $\beta$ and $b$ values. Furthermore, changing the value of $b$, we cannot obtain at all the points of the spatial domain the same behavior of the temperature above or below the reference curve ($b = 2000$ m$^{-\beta}$ and $\beta = 1$).

In the introduction, we mentioned some interpretations of the FBLL, although the physical meaning of the fractional-order $\beta$ is still an open problem ^[8]. However, $\beta$ in our study is only a parameter that can be adjusted to obtain a better agreement between the theoretical and experimental results.

4.   Conclusions

In this work, a new approach to the modeling of laser heating has been carried out. The study focuses on an alternative formulation of the laser heat source based on the FBLL. Introducing the solution of the FBLL into the heat transfer equation, we obtain different temperature profiles than those obtained with the classical BLL. These differences are also observed in all the points of the spatial domain, even in the heating phase, something that other alternative heat equation formulations (such as HHTE or FDHTE) did not show. Therefore, the FBLL is suitable for modeling laser heating, and it could be applied in cases where the BLL formulation does not fit the data well. The use of the FBLL is especially important in laser heating cases where the laser heat source is a term quantitatively larger than the rest of the terms of the heat transfer equation, since the possible non-ideal effects of the laser heating are more noticeable.

Although the FBLL solution has been obtained in previous studies, the novelty of the present work is that this solution is incorpored into a heat transfer problem, and the resulting equation is solved analytically. In fact, it is applied to obtain the temperature profiles produced in the laser heating of the cornea called LTK.

The results obtained can be extended to other techniques involving laser heating of biological tissue and, in general, other processes of laser heating that suggest the use of FBLL.

Acknowledgments

C. Lizama is partially supported by ANID Project FONDECYT 1220036. M. Murillo-Arcila is supported by MCIN/AEI/10.13039/501100011033, Project PID2019-105011GBI00, and by Generalitat Valenciana, Project PROMETEU/2021/070. M. Trujillo is supported by Grant RTI2018-094357-B-C21 funded by MCIN/AEI/10.13039/501100011033 (Spanish Ministerio de Ciencia, Innovación y Universidades / Agencia Estatal de Investigación)

Conflict of interest

All authors declare no conflicts of interest in this paper.