Semianalytical thermal model for subablative laser heating of homogeneous nonperfused biological tissue: application to laser thermokeratoplasty

Fabrice Manns; David Borja; Jean-Marie A. Parel; William E. Smiddy; William Culbertson

doi:10.1117/1.1560644

1 April 2003 Semianalytical thermal model for subablative laser heating of homogeneous nonperfused biological tissue: application to laser thermokeratoplasty

Fabrice Manns, David Borja, Jean-Marie A. Parel, William E. Smiddy, William Culbertson

Author Affiliations +

Journal of Biomedical Optics, Vol. 8, Issue 2, (April 2003). https://doi.org/10.1117/1.1560644

1. Introduction

The cornea and sclera of the human eye moderately absorb infrared radiation at wavelengths between 1500 to 2200 nm. Continuous-wave and pulsed infrared lasers emitting in this wavelength range, such as pulsed holmium:YAG and thulium:YAG lasers, and continuous-wave diode lasers, can be used to increase the corneal or scleral temperature above the threshold temperature for thermal shrinkage of collagen (55 to 65 °C). Laser thermal shrinkage of the cornea is used in laser thermokeratoplasty (LTK).¹ ² ³ ⁴ ⁵ ⁶ In LTK, the cornea is heated by applying 8 to 16 laser spots in an annular pattern to induce local shrinkage. Shrinkage induces a change in corneal shape that is used to correct hyperopia. Thermal shrinkage of the sclera has been used experimentally in laser scleral buckling (LSB).⁷ ⁸ In LSB, infrared laser radiation delivered through a fiber optic contact probe induces a local invagination of the sclera due to thermal shrinkage. This procedure is currently being investigated as an alternative to conventional buckle implantation for retinal detachment surgery.

The cornea and sclera consist mainly of water and collagen. Due to their high water content, the absorption curve of both tissues between 1500 and 2200 nm is expected to essentially overlap the absorption curve of water.² ⁹ An extrapolation of the wavelength dependence of the scleral scattering coefficient¹⁰ indicates that the effect of scattering on the light propagation in the cornea and sclera is negligible in this wavelength range. However, there is evidence that scattering may significantly increase during irradiation at wavelengths around 2000 nm, due to dynamic tissue changes induced by denaturation.¹¹ Together, moderate absorption and negligible scattering allow controlled and localized tissue heating, leading to efficient tissue shrinkage without vaporization.

The surgical outcome of LTK suffers from a lack of predictability and uncertain long-term stability, mainly because of a limited understanding of the fundamental thermodynamic and biomechanical processes of laser-induced corneal thermal shrinkage. The amount and kinetics of laser-induced thermal shrinkage or denaturation depend on the biochemical and physical properties of the irradiated tissue and on the irradiation parameters (spot size, laser wavelength, irradiation time, irradiance). In theory, the volume of shrinkage or denaturation can be estimated if the temperature distribution in the tissue and the time-temperature dependence of the denaturation process are known.¹² To provide a reasonable estimate of tissue denaturation, the temperature distribution in the irradiated tissue must be calculated over time with small time increments.

The temperature distribution in tissue is generally calculated by solving the bioheat equation using numerical algorithms (finite-difference or finite-element algorithms).¹² ¹³ ¹⁴ ¹⁵ ¹⁶ Numerical methods are especially advantageous in cases involving complex geometries, inhomogeneities, nonlinearities, complex boundary conditions, evaporation, and dynamic changes in optical or thermal properties. However, when the temperature distribution needs to be calculated over a large volume as a function of time, with small time and space increments, numerical solutions require prohibitive computation times, even with modern computers. Analytical expressions of the temperature distributions in laser-irradiated tissue and other materials were derived using the method of Green’s function.¹⁴ ¹⁷ ¹⁸ ¹⁹ The expression produced with this method is an integral over time and space that can be solved analytically only for a few simple cases.

We present a semianalytical solution of the heat equation for infrared laser irradiation of the cornea and sclera. The semianalytical solution relies on several assumptions that may limit its accuracy and applicability (see discussion), but it provides a rapid solution of the temperature distribution in laser-irradiated tissue at any point in time during and after irradiation. The solution is applicable to nonperfused homogeneous tissue in general, including scattering media, if the fluence rate is known. In Sec. 2, the general solution of the heat equation is presented assuming laser pulses with constant irradiance. In the most general case, the solution is an analytical expression that can be calculated numerically (i.e., semianalytical solution). The solution of the heat equation is also provided for a collimated beam with Gaussian or flat-top intensity distribution perpendicularly incident on a slab of nonscattering tissue. Numerical calculations for a Gaussian beam incident on the cornea are presented in Sec. 4. Implications for LTK and limitations of the model are discussed in Sec. 5.

2. Semianalytical Thermal Model

2.1.

Heat Equation and Geometry

The thermal model solves the linear bioheat equation in Cartesian coordinates under laser irradiation, in a homogeneous isotropic finite tissue slab of dimensions (a,b,c) in the (x,y,z) directions:

Eq. (1)

\frac{1}{α} \frac{\partial T}{\partial t} = (\frac{\partial^{2} T}{\partial x^{2}} + \frac{\partial^{2} T}{\partial y^{2}} + \frac{\partial^{2} T}{\partial z^{2}}) + \frac{g (x, y, z, t)}{k},

where T(x,y,z,t) is the temperature distribution in space and time, g(x,y,z,t) is the heat source term corresponding to absorption of laser radiation, and the other parameters are defined in Table 1. The bioheat equation was solved in a Cartesian coordinate system to provide added flexibility in the choice of irradiation and tissue geometries. For instance, the solution can be directly applied to calculate the temperature during laser heating experiments in rectangular tissue strips, or in tissue irradiated with a rectangular beam. In addition, the final analytical expression is simpler when a Cartesian (rather than cylindrical) coordinate system is used.

Table 1

Notation and parameters. The numerical values of the thermo-optical properties of the cornea were taken from the literature (see, for example, Refs. 9, 16, and 22).
Parameter	Description	Numerical value
a	Tissue x dimension	10 mm
b	Tissue y dimension	10 mm
c	Tissue thickness (z dimensions)	0.55 mm
k	Thermal conduction constant	0.556 W m⁻¹ °C⁻¹
ρ	Density	1 g cm⁻³
c_p	Heat capacity	3.83 J g⁻¹ °C⁻¹
α	Thermal diffusivity	0.001452 cm² s⁻¹
T₀	Initial temperature	35 °C
T_a	Ambient temperature	20 °C
h₁	Anterior convection constant	20 to 500 W m⁻² °C⁻¹
h₂	Posterior convection constant	20 to 2000 W m⁻² °C⁻¹
E₀	Irradiance	50031 W cm⁻²
t_p	Pulse duration	200 μs
w	Beam radius	0.3 mm
-	Energy per pulse per spot	30 mJ
-	Repetition rate	5 Hz
μ_a	Absorption coefficient	20 cm⁻¹
R	Fresnel surface reflectance	2.4

We assume that the beam is incident at normal incidence and that the z direction (depth) is parallel to the axis of the beam (Fig. 1). The initial temperature distribution is assumed to be homogeneous: T(x,y,z)=T₀ at t=0 for all x,y,z.

Figure 1

Irradiation and tissue geometry. The model assumes an isotropic, homogeneous finite tissue slab.

2.2.

Boundary Conditions

In the Cartesian coordinate system of Fig. 1, there are four tissue boundaries at x=0, x=a, y=0, and y=b, and the center of the beam is located at x=a/2, y=b/2. In most practical cases, such as local irradiation of the cornea or sclera, the spot radius (w) is much smaller than the lateral dimensions of the tissue (w≪a,b). Under these conditions, and as long as no external thermal source or forced convection is applied at the boundaries at x=0, x=a, y=0, and y=b, the corresponding boundary conditions will have no significant impact on the temperature distribution in the tissue. We solved the heat equation assuming that these four boundaries were insulated, and that the boundaries at z=0 and z=c are either insulated, at constant temperature, or subject to convection into ambient air for the boundary at z=0 (air-cornea boundary) and into aqueous for the boundary at z=c (cornea-aqueous boundary). We assume that air is at ambient temperature (T_a) and that the temperature of the aqueous is equal to the initial corneal temperature (T₀) (see Table 1).

2.3.

General Solution

A general analytical expression of the heat conduction problem of Eq. (1) can be found using the integral-transform technique or the Green’s function approach.²⁰ If we assume that the heat source term g(x,y,z,t) is constant during laser irradiation (cw irradiation or rectangular laser pulse), and that the irradiation time is t_p, integration over time of the general analytical solution of Eq. (1) yields:

Eq. (2)

T (x, y, z, t) = T_{0} + \frac{1}{k} \sum_{p = 1}^{\infty} {H_{a} (z, η_{p}) + \sum_{n = 0}^{\infty} \sum_{m = 0}^{\infty} H_{t} (t, β_{m}, γ_{n}, η_{p}) \cdot H_{xyz} (x, y, z, β_{m}, γ_{n}, η_{p})},

where the parameters β_m, γ_n, η_p are eigenvalues that are dependent on the boundary conditions, H_t(t,β_m,γ_n,η_p) and H_a(z,η_p) are given in Table 2, and H_xyz(x,y,z,β_m,γ_n,η_p) is given by:

Table 2

Expression of the H functions for collimated beams with Gaussian and flat-top intensity distributions.
	Gaussian beam	Flat-top beam
	$T (x, y, z, t) = T_{0} + \frac{1}{k} \sum_{p = 1}^{\infty} {H_{a} (z, η_{ρ}) + E_{0} \cdot (1 - R) \cdot μ_{a} \cdot H_{z} (z, η_{ρ}) \cdot \sum_{n = 0}^{\infty} \sum_{m = 0}^{\infty} H_{t} (t, β_{m}, γ_{n}, η_{ρ}) \cdot H_{x} (x, β_{m}) \cdot H_{y} (y, γ_{n})}$	$T (x, y, z, t) = T_{0} + \frac{1}{k} \sum_{p = 1}^{\infty} {H_{a} (z, η_{ρ}) + E_{0} \cdot (1 - R) \cdot μ_{a} \cdot H_{z} (z, η_{ρ}) \cdot \sum_{n = 0}^{\infty} \sum_{m = 0}^{\infty} H_{t} (t, β_{m}, γ_{n}, η_{ρ}) \cdot H_{x y} (x, y, β_{m}, γ_{n})}$
Time dependence H_t(t,β_m,γ_n,η_p)	$t \leq t_{p} : \underset{(m, n, p) \neq (0, 0, 0)}{H_{t} (t < t_{p}, β_{m}, γ_{n}, η_{p})} = \frac{1 - e^{- α (β_{m}^{2} + γ_{n}^{2} + η_{p}^{2}) t}}{β_{m}^{2} + γ_{n}^{2} + η_{p}^{2}},$ and H_t(t<t_p,β₀,γ₀,η₀)=α⋅t	$t > t_{p} : \underset{(m, n, p) \neq (0, 0, 0)}{H_{t} (t > t_{p}, β_{m}, γ_{n}, η_{p})} = e^{- α (β_{m}^{2} + γ_{n}^{2} + η_{p}^{2}) t \cdot} \frac{e^{α (β_{m}^{2} + γ_{n}^{2} + η_{p}^{2}) t_{p} - 1}}{β_{m}^{2} + γ_{n}^{2} + η_{p}^{2}},$ and H_t(t>t_p,β₀,γ₀,η₀)=α⋅t_p
Ambient temperature H_a(z,η_p)	$H_{a} (z, η_{ρ}) = h_{1} \cdot (T_{a} - T_{0}) \cdot \frac{2. (η_{p} \cos η_{p} z + \frac{h_{1}}{k} \sin η_{p} z)}{(η_{p}^{2} + \frac{h_{1}^{2}}{k^{2}}) (c + \frac{h_{2} \cdot k}{η_{p}^{2} k^{2} + h_{2}^{2}}) + \frac{h_{1}}{k}} \cdot \frac{1 - e^{- α \cdot η_{p}^{2} \cdot t}}{η_{p}}$
x−y dependence^* H_xy(x,y,β_m,γ_n) H_x(x,β_m),H_y(y,γ_n)	$\begin{array}{l} H_{x} (X, β_{2 m}) = \frac{ε}{a} \cdot \sqrt{2} \cdot w \cdot \cos β_{2 m} X \cdot {(- 1)}^{m} \cdot I (2 m, \frac{a}{\sqrt{2} \cdot w}) \\ H_{x} (X, β_{2 m + 1}) = 0 \\ ε = 1 if m=0, ε =2 if m \neq 0 \\ H_{y} (Y, γ_{2 n}) = \frac{ε}{b} \cdot \sqrt{2} \cdot w \cdot \cos γ_{2 n} Y \cdot {(- 1)}^{n} \cdot I (2 n, \frac{b}{\sqrt{2} \cdot w}) \\ H_{y} (Y, γ_{2 n + 1}) = 0 \\ ε = 1 if n=0, ε =2 if n \neq 0 \end{array}$	$\underset{m \neq 0, n \neq 0}{H_{x y} (X, Y, β_{2 m}, γ_{2 n})} = \frac{8 π w}{a^{2} \cdot {(β_{2 m}^{2} + γ_{2 n}^{2})}^{1 / 2}} \cdot {(- 1)}^{m + n} \cos β_{2 m} X \cdot \cos γ_{2 n} Y \cdot J_{1} [{(β_{2 m}^{2} + γ_{2 n}^{2})}^{1 / 2} \cdot W]$ H_xy(x,y,β_2m+1,0)=0, H_xy(x,y,0,γ_2n+1)=0, H_xy(x,y,β_2m+1,γ_2n+1)=0
x−y dependence^* H_xy(x,y,β_m,γ_n) H_x(x,β_m),H_y(y,γ_n)	H_y(y,γ_2n+1)=0

z dependence^†	z=0: k∂T/∂z=h₁T z=c: k∂T/∂z=−h₂T
$H_{z} (z, η_{ρ}) = \frac{2. (η_{p} \cos η_{p} z + \frac{h_{1}}{k} \sin η_{p} z)}{(η_{p}^{2} + \frac{h_{1}^{2}}{k^{2}}) (c + \frac{h_{2} \cdot k}{η_{p}^{2} k^{2} + h_{2}^{2}}) + \frac{h_{1}}{k}} \cdot \frac{η_{ρ} \cdot (μ_{a} + \frac{h_{1}}{k}) \cdot (1 - e^{- μ_{a} c} \cdot \cos η_{ρ} c) + e^{- μ_{a} c} \cdot \sin η_{ρ} c \cdot (η_{p}^{2} - μ_{a} \frac{h_{1}}{k})}{μ_{a}^{2} + η_{p}^{2}}$
	z=0: ∂T/∂z=0 z=c: ∂T/∂z=0 $if p \neq 0: Η_{z} (z, η_{ρ}) = \frac{2}{c} \cdot \cos η_{ρ} z \cdot \frac{μ_{a}}{μ_{a}^{2} + η_{ρ}^{2}} [1 - e^{- μ_{a^{c}}} \cdot \cos η_{ρ} c]$ $if p=0: Η_{z} (z, η_{0}) = \frac{1}{c \cdot μ_{a}} [1 - e^{- μ_{a^{c}}}]$
^*J₁(x) is the Bessel function of the first kind of order 1 and the function I(q,L), where q is a positive integer and L is a real number, is defined as: $I (q, L) = \int_{u = 0}^{L} \cos (q π \frac{u}{2 L}) \cdot e^{- u^{2}} \cdot d u \cdot$
^† The eigenfunctions, norms, and eigenvalues for all other boundary conditions at z=0 and z=c can be found by setting h_i=0 for an insulated boundary and h_i=∞ for a boundary at fixed temperature (where h_i=h₁ or h₂) .

Eq. (3)

\begin{matrix} H_{xyz} (x, y, z, β_{m}, γ_{n}, η_{p}) = \frac{X (β_{m}, x) Y (γ_{n}, y) Z (η_{p}, z)}{N (β_{m}) N (γ_{n}) N (η_{p})} \cdot \int_{x^{'} = 0}^{a} \int_{y^{'} = 0}^{b} \int_{z^{'} = 0}^{c} \\ X (β_{m}, x^{'}) Y (γ_{n}, y^{'}) Z (η_{p}, z^{'}) \cdot g (x^{'}, y^{'}, z^{'}) \times d x^{'} {dy}^{'} {dz}^{'} \cdot \end{matrix}

The functions X,Y,Z are eigenfunctions and the functions N are norms of the three differential equations obtained when solving Eq. (1) by separation of variables.²⁰ These functions depend on the boundary conditions. The expressions of the eigenfunctions, norms, and eigenvalues corresponding to the boundary conditions of our problem are listed in Table 3. Note that each of the sums in Eq. (2) starts at 0 if both boundaries in the corresponding direction are insulated. In all other cases, the sums start at 1. Equations (2) and (3) are an analytical expression of the temperature distribution, but it can only be calculated numerically, save for a few simple cases.

Table 3

Boundary conditions, eigenvalues, norms and eigenfunctions. For z, the eigenfunctions, norms, and eigenvalues for all other boundary conditions at z=0 and z=c can be found by setting hi=0 for an insulated boundary and hi=∞ for a boundary at a fixed temperature (where hi=h1 or h2) .
Coordinate	Boundary condition	Eigenfunction, norm, eigenvalue
x	x=0: ∂T/∂x=0	X(β_m,x)=cos(β_mx)
	x=a: ∂T/∂x=0	N(β_m)=a/2 if m≠0,N(β_m)=a if m=0
		β_m=mπ/a(m⩾0)
y	y=0: ∂T/∂y=0	Y(γ_n,y)=cos(γ_ny)
	y=b: ∂T/∂y=0	N(γ_n)=b/2 if n≠0,N(γ_n)=b if n=0
		γ_n=nπ/b(n⩾0)
z^*	z=0: k∂T/∂z=h₁(T−T_a)	$Z (η_{ρ}, z) = η_{ρ} \cdot \cos (η_{ρ} z) + \frac{h_{1}}{k} \cdot \sin (η_{ρ} z)$
	z=c: k∂T/∂z=−h₂T	$N (η_{ρ}) = \frac{1}{2} \cdot [\frac{h_{1}}{k} + (η_{ρ}^{2} + \frac{h_{1}^{2}}{k^{2}}) \cdot (c + \frac{k \cdot h_{2}}{η_{ρ}^{2} k^{2} + h_{2}^{2}})]$
		The eigenvalues are all positive roots of the equation:
		$\tan (η_{ρ} c) = \frac{η_{ρ} \cdot k \cdot (h_{1} + h_{2})}{η_{ρ}^{2} k^{2} - h_{1} h_{2}}$

2.4.

Expression of H _xyz for Collimated Irradiation and No Scattering

For irradiation with a collimated beam at normal incidence and with constant irradiance, when the effect of scattering is negligible, the heat source term is given by:

Eq. (4)

g (x, y, z) = (1 - R) \cdot μ_{a} \cdot E_{i} (x, y) \cdot exp (- μ_{a} z),

where E_i(x,y) is the incident irradiance distribution, R is the Fresnel tissue surface reflectance, and μ_a is the absorption coefficient. Under these conditions, the function H_xyz(x,y,z,β_m,γ_n,η_p) can be separated into a product of two functions of separate variables:

Eq. (5)

H_{xyz} (x, y, z, β_{m}, γ_{n}, η_{p}) = (1 - R) \cdot μ_{a} \cdot H_{z} (z, η_{p}) \cdot H_{xy} (x, y, β_{m}, γ_{n}),

where

Eq. (6)

H_{z} (z, η_{p}) = \frac{Z (η_{p}, z)}{N (η_{p})} \cdot \int_{z^{'} = 0}^{c} Z (η_{p}, z^{'}) \cdot exp (- μ_{a} z^{'}) d z^{'},

and

Eq. (7)

H_{xy} (x, y, β_{m}, γ_{n}) = \frac{X (β_{m}, x) \cdot Y (γ_{n}, y)}{N (β_{m}) N (γ_{n})} \cdot \int_{x^{'} = 0}^{a} \int_{y^{'} = 0}^{b} X (β_{m}, x^{'}) \cdot Y (γ_{n}, y^{'}) \cdot E_{i} (x^{'}, y^{'}) d x^{'} {dy}^{'} \cdot

The function H_z(z,η_p) is independent on the beam parameters and depends only on the absorption coefficient and on the boundary conditions at z=0 and z=c. The expressions of this function for selected boundary conditions are listed in Table 2. On the other hand, since we fixed the boundary conditions in the x and y directions, the function H_xy(x,y,β_m,γ_n) depends only on the incident intensity distribution.

2.5.

Expression of H_xy for a Beam with Symmetric Gaussian Intensity Distribution

The incident irradiance of a symmetric Gaussian beam centered on the tissue slab is given by:

Eq. (8)

E_{i} (x, y) = E_{0} \cdot e^{- 2 {(x - a / 2)}^{2} / w^{2}} \cdot e^{- 2 {(y - b / 2)}^{2} / w^{2}},

where E₀ is the peak irradiance and w is the 1/e² spot radius. The peak irradiance, instantaneous power P, and beam radius w, are related by: E₀=2⋅P/π⋅w². Combining Eqs. (8) and (7), we find that the function H_xy can be separated into a product of two functions with separate variables:

Eq. (9)

H_{xy} (x, y, β_{m}, γ_{n}) = E_{0} \cdot H_{x} (x, β_{m}) \cdot H_{y} (y, γ_{n}) \cdot

The expressions of the functions H_x and H_y calculated for the selected boundary conditions (see Table 3) are given in Table 2.

2.6.

Expression of H_xy for a Circular Beam with Uniform Intensity Distribution

The incident irradiance of a circular beam with uniform intensity centered on the tissue slab can be written:

Eq. (10)

E_{i} (x, y) = E_{0} \cdot u [{(w^{2} - y^{2})}^{1 / 2}, \frac{a}{2}, x] \cdot u (w, \frac{b}{2}, y),

where the irradiance (E₀), laser power (P) and beam radius (w) are related by E₀=P/π⋅w² and u(ɛ,d,s) is a one-dimensional rectangular function of half width ɛ centered at s=d. The resulting expression of the function H_xy, calculated for the selected boundary conditions (see Table 3), is given in Table 2. For a circular flat-top beam, the function H_xy cannot be separated into functions of separate variables.

2.7.

Pulsed Irradiation

According to the principle of superposition, the temperature during repetitively pulsed irradiation is the convolution of the temperature distribution produced by one pulse with the Dirac comb function corresponding to the laser repetition rate. For instance, if pulses are applied every second, the temperature at a given point after 2.7 s is the sum of the temperatures produced by a single pulse at times of 2.7 s (contribution of the first pulse), 1.7 s (contribution of the second pulse), and 0.7 s (contribution of the third pulse) after the start of the single pulse.

2.8.

Summary

The calculations show that the temperature distribution produced in a homogeneous, isotropic, purely absorbing tissue slab by a Gaussian collimated beam is the triple sum of the product of four functions of separate variables (H_x,H_y,H_z,H_t), plus a simple sum of an analytical function (H_a) that is present because the boundary condition at z=0 is inhomogeneous (T_a≠T₀). Two of the four functions in the triple sum have analytical expressions. The other two functions are simple integrals with finite boundaries that must be calculated numerically. A similar expression is found when the intensity distribution is uniform, except that the variables x and y can no longer be separated.

3. Numerical Implementation

The temperature during irradiation with a Gaussian beam was calculated with a computer program written in Turbo C/C++. The program calculates the triple sum of the product of the four functions, H_x, H_y, H_z, and H_t in three successive steps, starting by the sum over m, and finishing with the sum over p. Each of the sums is stopped when the value of an asymptote of the summated function falls below a preset convergence limit (the summated functions all converge to zero for large values of m, n, or p). A convergence limit of 0.001 was found to provide sufficient accuracy (error less than 0.1 °C).

When the boundaries at z=0 and/or z=c are subject to convection, the eigenvalues η_p are calculated using a modified Newton-Raphson algorithm. For repetitive calculations with the same values of the absorption coefficient and convection constant, the eigenvalues were first calculated with a separate program and fit with a nonlinear function. This nonlinear function was then used in the program calculating the temperature to calculate the eigenvalues. Curve fitting the eigenvalues significantly reduces the computation time.

The integral function I(q,L) in Table 2 was calculated numerically using the method of trapezoids. An integration increment equal to 0.001L/q (10⁻³ periods of the cosine function in the integral) was used. Since the exponential function in the integral I(q,L) rapidly converges to zero when L increases, I(q,L) was approximated by the integral from 0 to infinity for values of L larger than 5:²¹

Eq. (11)

\int_{x = 0}^{\infty} cos (q π \frac{x}{2 L}) \cdot e^{- x^{2}} dx = \frac{\sqrt{π}}{2} e^{- {(q π / 4 L)}^{2}} \cdot

With these parameters, the relative error in the function I(q,L) was estimated to be less than 0.05 (absolute error less than 10⁻⁵) .

The values of the dimensions a and b must be large enough to avoid that the boundary conditions at x=0, x=a, y=0, and y=b have an influence on the temperature in the volume heated by laser irradiation. Preliminary calculations showed that the temperature increase produced by absorption of laser radiation and lateral diffusion of heat is insignificant at positions that are away from the center by more than three to four times the beam radius. We always used a value of a and b at least five times larger than the beam radius.

4. Application to Pulsed Ho:YAG Laser ThermoKeratoplasty

4.1.

Irradiation Parameters and Effects of the Convection Constants

We used the model to estimate the corneal temperature during Ho:YAG laser thermokeratoplasty with irradiation parameters that are typical of a clinical treatment (Hyperion LTK system, Sunrise Technologies, Fremont, California). The values of all laser and tissue parameters, except the convection constants, are listed in Table 1. The convection constant at the anterior corneal surface is generally accepted to be on the order of 20 W m⁻² K⁻¹.²² To the best of our knowledge, there is no experimental data available on the effect of aqueous flow on the corneal temperature. In their calculations of corneal temperature during radiofrequency heating, Berjano, Saiz, and Ferrero²² modeled the effect of aqueous flow by assuming that the posterior corneal surface is subject to free convection with constants ranging from h₂=20 to 1000 W m⁻² °C⁻¹. They found that variations of the convection constant in this range did not have a significant effect on the temperature distribution (less than 4 °C). On the other hand, because the aqueous flow rate is slow (on the order of μ1/min),²³ earlier thermal models assumed that the convective effect of aqueous flow is negligible.²⁴ Since the thermal and optical properties of the aqueous are essentially equal to those of the cornea, these models assumed that there is no discontinuity at the cornea-aqueous boundary, and that the cornea and aqueous can be treated as a homogeneous continuous structure for the purpose of thermal modeling.²⁴

To determine the effect of the value of the convection constants on the temperature distribution during laser heating, we calculated the corneal temperature at the anterior or posterior corneal surface at the center of the 0.6-mm-diam beam (peak temperature) as a function of time during the treatment, as well as the depth profile of the corneal temperature at the center of the beam at the end of the last pulse (t=1.2002 s) with different values of the convection constants. First, h₁ was fixed (h₁=20 W m ⁻² K ⁻¹) and h₂ was varied from 20 to 2000 W m⁻² K⁻¹. Second, h₂ was fixed (h₂=1000 W m ⁻² K ⁻¹) and h₁ was varied from 20 to 500 W m⁻² K⁻¹. We also calculated the temperature with h₁=20 W m ⁻² K ⁻¹ and the assumption that the cornea and aqueous form a homogeneous and continuous structure, with a thickness c=3 mm (approximately equal to the human anterior chamber depth). Since the penetration depth in the cornea of radiation at 2.1 μm is on the order of 500 μm (μ_a=20 cm ⁻¹), the temperature increase at z=3 mm is expected to be essentially equal to zero. We therefore assumed that the boundary at z=c was insulated (∂T/∂z=0) when the discontinuity at the cornea-aqueous boundary was ignored.

These calculations showed that variations of h₁ do not have a significant effect on the depth and time dependence of the corneal temperature during the treatment, except for h₁=500 W m ⁻² K ⁻¹, which is an extreme value (Fig. 2). On the other hand, variations of h₂ produced significant changes in both the depth and time dependence of the corneal temperature (Fig. 3), especially at the cornea-aqueous boundary. The assumption that the cornea and aqueous form a homogeneous structure produced a temperature distribution that is very close to the temperature distribution produced when h₂=1000 W m ⁻² K ⁻¹. These results are further discussed in Sec. 5.

Figure 2

Effect of the convection constant at the anterior corneal surface (h₁) on the anterior corneal surface temperature during pulsed Ho:YAG laser LTK (top) and on the depth profile at the end of the last laser pulse (t=1.2002 s, bottom).

Figure 3

Effect of the convection constant at the posterior corneal surface (h₂) on the posterior corneal surface temperature during pulsed Ho:YAG laser LTK (top) and on the depth profile at the end of the last laser pulse (t=1.2002 s, bottom). The bold continuous curve in the depth profile is obtained with the assumption that the cornea and aqueous form a continuous homogeneous structure.

4.2.

Effect of Ambient Temperature

According to the expression of H_a in Table 2, the effect on the temperature distribution of the value of the ambient temperature is directly proportional to h₁ and (T_a−T₀). The effect increases with time, and is maximal at the anterior corneal surface (z=0). To determine the effect of ambient temperature on the temperature distribution inside the cornea, the anterior corneal surface temperature (z=0) during Ho:YAG LTK was calculated after 1.4 s for h₂=1000 W m ⁻² K ⁻¹ with values of h₁ ranging from 20 to 500 W m⁻² K⁻¹ and values of T_a ranging from 15 to 35 °C. For T_a=15 °C, the anterior corneal surface temperature after 1.4 s was 61.3, 59.2, and 51.2 °C for h₁=20, 100, and 500 W m⁻² K⁻¹, respectively. For T_a=35 °C, the anterior corneal surface temperature after 1.4 s for the same values of h₁ was 61.6, 60.4 and 55.5 °C. These results show that ambient temperature has a negligible effect (at most on the order of 1 °C) on the calculated temperature distribution during LTK for h₁ ranging from 20 to 100 W m⁻² K⁻¹. The effect becomes significant (>4 °C) only in the extreme case when h₁=500 W m ⁻² K ⁻¹.

4.3.

Corneal Temperature Distribution

The corneal temperature distribution in a plane passing through the center of the beam was calculated at the end of each of the seven pulses, assuming h₁=20 W m ⁻² K ⁻¹ and h₂=1000 W m ⁻² K ⁻¹, and T_a=20 °C (Fig. 4). Figure 4 also shows the temperature at different depths as a function of time during the treatment.

Figure 4

Temperature profiles through the center of the treatment spot at the end of each laser pulse during Ho:YAG LTK and time dependence of the anterior corneal surface temperature (bottom right). The treatment parameters are listed in Table 1. The temperature profiles were calculated for h₁=20 W m ⁻² °C ⁻¹ and h₂=1000 W m ⁻² °C ⁻¹.

5. Discussion

5.1.

Thermal Model

We present an exact semianalytical model that allows a rapid estimation of laser-induced temperature fields in tissue. The model assumes irradiation at constant power (continuous wave or rectangular pulse) of a homogeneous tissue slab subject to free convection, with no perfusion, evaporation, or mass transfer. In its most general form the model can be used to calculate the temperature in both scattering and nonscattering media, as long as the fluence rate in the tissue is known. An analytical expression of the temperature distribution can be found for purely absorbing tissue irradiated with a collimated beam. In its present form, the model is valid as long as the tissue dimensions are much larger than the spot diameter. Otherwise, the assumption that the lateral boundaries are insulated may limit the accuracy of the calculated temperature. However, the solution can be modified to provide the temperature in a finite rectangular tissue strip with all boundaries subject to convections by using the corresponding eigenfunctions, eigenvalues, and norms.²⁰

The main advantage of the analytical technique compared to a numerical technique is that it provides a rapid and accurate mathematical solution of the bioheat equation [Eq. (1)], even during the laser pulse, and even for strong temperature gradients. In the end, however, the accuracy of the calculated temperature will depend mainly on how accurately the tissue temperature can be modeled with the linear bioheat equation using simple convective boundary conditions. The analytical solution cannot easily be adapted to take into account complex geometries, inhomogeneities, dynamic variations of the thermal and optical properties, and phase changes. Some of these factors may affect the accuracy of the calculated corneal temperature during Ho:YAG LTK, especially when the temperature exceeds 60 to 70 °C. During LTK, the cornea shrinks and probably dehydrates. Shrinking and dehydration may induce a significant variation of the optical (absorption, scattering, and refractive index) and thermal properties of the cornea during irradiation. Variations of the absorption coefficient during heating may induce significant errors in the calculated corneal temperature.²⁵ ²⁶ Ith, Frenz, and Weber¹¹ also recently demonstrated that tissue changes during irradiation can induce variations of the tissue scattering properties that can significantly affect the light distribution, even during irradiation with midinfrared lasers. Torres et al.;¹⁵ showed that surface evaporation and vaporization may also significantly affect the thermal response of laser irradiated tissue, even at temperatures as low as 70 °C. Combined, all of these effects probably limit the absolute accuracy of the model when the temperature exceeds 60 to 70 °C.

Other factors that may affect the accuracy of the thermal model for LTK are that the beam does not necessarily have a perfectly Gaussian intensity distribution, and that the laser pulse is not rectangular. The temporal shape of the laser pulse will affect the heating rate during laser irradiation, but it is not expected to have a significant effect on the final temperature distribution, as long as the total energy delivered per pulse remains constant. Clinical LTK systems use fiber optic delivery. The intensity distribution at the corneal surface is generally assumed to be close to Gaussian. As shown by Eqs. (5), (6), and (7), for a collimated beam, the intensity distribution affects only the x and y dependence of the temperature distribution in the cornea, and not the depth dependence. If a more accurate knowledge of the lateral dependence of the temperature distribution is required, the exact intensity profile can be measured and entered in the general expression of the temperature distribution given by Eq. (7).

5.2.

Corneal Temperature During LTK

The calculation of the temperature distribution in the cornea during pulsed Ho:YAG LTK shows that for the duration of the treatment (less than 2 s), free convection at the anterior corneal surface has a negligible effect on the temperature distribution. Except in extreme cases, variations of h₁ and T_a have no significant effect. Therefore, assuming that the anterior corneal surface is insulated would not induce a significant error in the temperature calculation. On the other hand, the values and depth profile of the corneal temperature distribution were found to be strongly dependent on the value of the convection constant at the aqueous-cornea boundary. A better knowledge of the cooling effect of the aqueous is required to more accurately estimate the corneal temperature distribution, and especially to evaluate the risk of damage to corneal endothelium, which does not regenerate. The assumption that the cornea and aqueous form a homogeneous structure produces a temperature distribution similar to the one obtained when the convection constant is h₂=1000 W cm ⁻² K ⁻¹. The temperature distributions in these two cases were similar, most likely because the temperature gradients at the aqueous-cornea boundary produced with these two assumptions are approximately equal. The value of the convection constant h₂ that produces a temperature gradient at the cornea-aqueous boundary equal to the gradient produced when the discontinuity is ignored will be called the “equivalent convection constant.” It is important to note that the equivalent convection constant is dependent on the temperature gradient at the aqueous-cornea boundary. For instance, if a different wavelength with a significantly higher absorption coefficient in the cornea is used, the temperature gradient at the cornea-aqueous boundary will be smaller when the discontinuity is ignored, because the absorbed fluence will decay much more rapidly with depth. In that case, we expect that the equivalent convection constant will be significantly less than h₂=1000 W cm ⁻² K ⁻¹.

Even though the accuracy of the model may be limited for temperatures exceeding 70 °C, (see Sec. 5.1), the calculations indicate that corneal temperature reaches values at the surface and at the endothelium that may be sufficient to induce local superficial vaporization of the epithelium and/or local endothelial cell damage. These findings agree with the results of several experimental studies that showed that there is a risk for local endothelial cell damage²⁷ and local epithelial and anterior stromal scarring²⁸ during pulsed LTK if the energy or number of pulses delivered are too high.

The model also demonstrates that there is significant lateral heat diffusion between the laser pulses. For instance, the position of the isothermal contour at 60 °C at the corneal surface progresses from approximately 0.18 mm from the center of the spot at the end of the first pulse to approximately 0.30 mm from the center of the spot at the end of the last pulse. The temperature distribution progressively reaches a steady state toward the end of the treatment.

Our calculations of the temperature during LTK agree at least qualitatively with experimental temperature measurements using a thermal camera.²⁹ ³⁰ Nevertheless, the model needs to be further validated by directly comparing our predictions with temperature measurements using thermal radi- ometry (for surface temperature measurements) or temperature probes inserted in the cornea during laser irradiation.

6. Conclusions

A semianalytical solution of the bioheat equation during laser irradiation of a finite homogeneous isotropic purely absorbing slab of tissue that allows rapid estimation of the temperature distribution at any point in time is presented. The solution was used to calculate the corneal temperature during laser thermokeratoplasty (LTK). The model demonstrates that the corneal temperature during LTK with a pulsed Ho:YAG laser may be sufficient to induce superficial vaporization of epithelial cells and local thermal damage to the endothelium. A better knowledge of the cooling effect of the aqueous is required to more accurately estimate the corneal temperature distribution during LTK. The model remains to be validated experimentally to determine the effect of shrinkage, evaporation, and other laser-induced tissue changes on the accuracy of the predicted temperatures.

Acknowledgments

This work was supported by the Whitaker Foundation through a Biomedical Engineering Research Grant. This work was also supported by The Henri and Flore Lesieur Foundation and unrestricted grants from The Florida Lion’s Eye Bank and Research to Prevent Blindness New York, New York.

REFERENCES

1.

R. Brinkmann , B. Radt , C. Flamm , J. Kampmeier , N. Koop , and R. Birngruber , “Influence of temperature and time on thermally induced forces in corneal collagen and the effect on laser thermokeratoplasty,” J. Cataract Refractive Surg. , 26 744 –754 (2000). Google Scholar

2.

T. Bende , B. Jean , and T. Oltrup , “Laser thermal keratoplasty using a continuous wave diode laser,” J. Refract. Surg. , 15 154 –158 (1999). Google Scholar

3.

D. D. Koch , T. Kohnen , P. J. McDonnell , R. Menefee , and M. Berry , “Hyperopia correction by noncontact holmium:YAG laser thermal keratoplasty: U.S. phase IIA clinical study with 2-year follow-up,” Ophthalmology , 104 1938 –1947 (1997). Google Scholar

4.

J. M. Parel , Q. Ren , and G. Simon , “Non-contact laser photothermal keratoplasty (LPTK). I. Biophysical principles and laser system,” Refract. Corneal Surg. , 10 511 –518 (1994). Google Scholar

5.

T. Seiler , M. Matallana , and T. Bende , “Laser thermokeratoplasty by means of a pulsed holmium:YAG laser for hyperopic correction,” Refract. Corneal Surg. , 6 335 –339 (1990). Google Scholar

6.

C. Smithpeter , E. Chan , S. Thomsen , H. G. Rylander , and A. J. Welch , “Corneal photocoagulation with continuous wave and pulsed holmium:YAG radiation,” J. Cataract Refractive Surg. , 21 258 –267 (1995). Google Scholar

7.

M. Sasoh , J. M. Parel , F. Manns , I. Nose , J. Comander , and W. E. Smiddy , “Quantification of Ho:YAG and Tm:YAG laser-induced scleral shrinkage for buckling procedures,” Ophthalmic Surg. Lasers , 29 410 –421 (1998). Google Scholar

8.

Q. Ren , G. Simon , J. M. Parel , and W. E. Smiddy , “Laser scleral buckling for retinal detachment,” Am. J. Ophthalmol. , 115 758 –762 (1993). Google Scholar

9.

R. Brinkmann , N. Koop , G. Geerling , J. Kampmeier , S. Borcherding , K. Kamm , and R. Birngruber , “Diode laser thermokeratoplasty: Application strategy and dosimetry,” J. Cataract Refractive Surg. , 24 1195 –1207 (1998). Google Scholar

10.

A. Vogel , C. Dlugos , R. Nuffer , and R. Birngruber , “Optical properties of the sclera and their consequences for transscleral laser applications,” Lasers Surg. Med. , 11 331 –340 (1991). Google Scholar

11.

M. Ith , M. Frenz , and H. P. Weber , “Scattering and thermal lensing of 2.12-μm laser radiation in biological tissue,” Appl. Opt. , 40 2216 –2223 (2001). Google Scholar

12.

F. Manns , D. Borja , and J. M. Parel , “Calculation of corneal temperature and shrinkage during Laser ThermoKeratoplasty (LTK),” Proc. SPIE , 4611 101 –109 (2002). Google Scholar

13.

C. Sturesson and S. Andersson-Engels , “A mathematical model for predicting the temperature distribution in laser-induced hyperthermia. Experimental evaluation and applications,” Phys. Med. Biol. , 40 2037 –2052 (1995). Google Scholar

14.

J. Roider and R. Birngruber, “Solution of the heat conduction equation,” in Optical-Thermal Response of Laser-Irradiated Tissue, van Gemert and Welch, Eds., pp. 385–409, Plenum Press, New York (1995).

15.

J. H. Torres , M. Motamedi , J. A. Pearce , and A. J. Welch , “Experimental evaluation of mathematical models for predicting the thermal response of tissue to laser radiation,” Appl. Opt. , 32 597 –606 (1993). Google Scholar

16.

R. Brinkmann , G. Dro¨ge , N. Koop , A. Wo¨rdemann , G. Schirner , and R. Birngruber , “Investigation on laser thermokeratoplasty,” Lasers Light Ophthal. , 6 259 –270 (1994). Google Scholar

17.

M. K. Loze and C. D. Wright , “Temperature distributions in laser-heated biological tissue with application to birthmark removal,” J. Biomed. Opt. , 6 74 –85 (2001). Google Scholar

18.

M. K. Loze and C. D. Wright , “Temperature distributions in laser-heated semi-infinite and finite-thickness media with convective surface losses,” Appl. Opt. , 37 6822 –6832 (1998). Google Scholar

19.

M. K. Loze and C. D. Wright , “Temperature distributions in semi-infinite and finite-thickness media as a result of absorption of laser light,” Appl. Opt. , 36 494 –507 (1997). Google Scholar

20.

M. N. Ozisik, Heat Conduction, 2nd ed., Wiley and Sons, New York (1993).

21.

I. S. Gradshteyn and I. M. Ryzhik, Table of Integrals, Series and Products, 5th ed., p. 515, Academic Press, New York (1994).

22.

E. J. Berjano , J. Saiz , and J. M. Ferrero , “Radio-frequency heating of the cornea: Theoretical model and in vitro experiments,” IEEE Trans. Biomed. Eng. , 49 196 –205 (2002). Google Scholar

23.

G. O. Sperber and A. Bill , “A method for near-continuous determination of aqueous humor flow; Effects of anaesthetics, temperature and indomethacin,” Exp. Eye Res. , 39 435 –453 (1984). Google Scholar

24.

J. A. Scott , “A finite element model of heat transport in the human eye,” Phys. Med. Biol. , 33 227 –241 (1988). Google Scholar

25.

F. Manns and J. M. Parel , “Scleral heating with pulsed mid-infrared lasers and temperature-dependent absorption coefficient,” Proc. SPIE , 2971 55 –58 (1997). Google Scholar

26.

E. D. Jansen , T. G. van Leeuwen , M. Motamedi , C. Borst , and A. J. Welch , “Temperature-dependence of the absorption coefficient of water for mid-infrared laser radiation,” Lasers Surg. Med. , 14 258 –268 (1994). Google Scholar

27.

C. Wirbelauer , N. Koop , A. Tuengler , G. Geerling , R. Bringruber , H. Laqua , and R. Brinkmann , “Corneal endothelial cell damage after experimental diode laser thermokeratoplasty,” J. Refract. Surg. , 16 323 –329 (2000). Google Scholar

28.

D. D. Koch , T. Kohnen , J. A. Anderson , P. S. Binder , M. N. Moore , R. F. Menefee , G. L. Valderamma , and M. J. Berry , “Histologic changes and wound healing response following 10-pulse noncontact holmium:YAG laser thermal keratoplasty,” J. Refract. Surg. , 12 623 –634 (1996). Google Scholar

29.

T. Papaionnaou , E. Maguen , and W. S. Grundfest , “Spatiotemporal temperature profiling of corneal surface during LTK,” Proc. SPIE , 4611 110 –114 (2002). Google Scholar

30.

J. Rhead , G. Boreman , A. Kar , F. Manns , M. Sasoh , and J. M. Parel , “Infrared measurement of thermal constants in laser-irradiated scleral tissue,” Proc. SPIE , 2673 77 –88 (1996). Google Scholar

Citation Download Citation

Fabrice Manns, David Borja, Jean-Marie A. Parel, William E. Smiddy, and William Culbertson "Semianalytical thermal model for subablative laser heating of homogeneous nonperfused biological tissue: application to laser thermokeratoplasty," Journal of Biomedical Optics 8(2), (1 April 2003). https://doi.org/10.1117/1.1560644

Published: 1 April 2003

Access the abstract

JOURNAL ARTICLE
10 PAGES

DOWNLOAD PAPER SAVE TO MY LIBRARY

GET CITATION

CITATIONS

Cited by 35 scholarly publications.

Explore citations on Lens.org

RIGHTS & PERMISSIONS

Get copyright permission Get copyright permission on Copyright Marketplace

KEYWORDS

Tissues

Convection

Cornea

Laser tissue interaction

Thermal modeling

Pulsed laser operation

Absorption

1.

Introduction

2.