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).^{1} ^{2} ^{3} ^{4} ^{5} ^{6} 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).^{7} ^{8} 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.^{2} ^{9} An extrapolation of the wavelength dependence of the scleral scattering coefficient^{10} 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.^{11} 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.^{12} 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).^{12} ^{13} ^{14} ^{15} ^{16} 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.^{14} ^{17} ^{18} ^{19} 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:
(1)
$$\frac{1}{\alpha}\text{}\frac{\partial T}{\partial t}={\displaystyle (\frac{{\partial}^{2}\mathit{T}}{\partial {x}^{2}}+\frac{{\partial}^{2}\mathit{T}}{\partial {y}^{2}}+\frac{{\partial}^{2}\mathit{T}}{\partial {z}^{2}})}+\frac{g(x,y,z,\mathit{t})}{k},$$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^{−1} °C^{−1} |
ρ | Density | 1 g cm^{−3} |
c_{p} | Heat capacity | 3.83 J g^{−1} °C^{−1} |
α | Thermal diffusivity | 0.001452 cm^{2} s^{−1} |
T_{0} | Initial temperature | 35 °C |
T_{a} | Ambient temperature | 20 °C |
h_{1} | Anterior convection constant | 20 to 500 W m^{−2} °C^{−1} |
h_{2} | Posterior convection constant | 20 to 2000 W m^{−2} °C^{−1} |
E_{0} | Irradiance | 50031 W cm^{−2} |
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^{−1} |
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_{0} at t=0 for all x,y,z.
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_{0}) (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.^{20} 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:
(2)
$$T(x,y,z,\mathit{t})={\mathit{T}}_{0}+\frac{1}{k}\text{}\sum _{p=1}^{\infty}\left\{{\mathit{H}}_{\mathit{a}}(\mathit{z},{\eta}_{\mathit{p}})+\sum _{n=0}^{\infty}\sum _{m=0}^{\infty}{\mathit{H}}_{\mathit{t}}(\mathit{t},{\beta}_{\mathit{m}},{\gamma}_{\mathit{n}},{\eta}_{\mathit{p}})\cdot {\mathit{H}}_{\mathit{xyz}}(\mathit{x},\mathit{y},z,{\beta}_{\mathit{m}},{\gamma}_{\mathit{n}},{\eta}_{\mathit{p}})\right\},$$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}{\displaystyle \sum _{p=1}^{\infty}\left\{{H}_{a}(z,{\eta}_{\rho})+{E}_{0}\cdot (1-R)\cdot {\mu}_{a}\cdot {H}_{z}(z,{\eta}_{\rho})\cdot {\displaystyle \sum _{n=0}^{\infty}}{\displaystyle \sum _{m=0}^{\infty}{H}_{t}(t,{\beta}_{m},{\gamma}_{n},{\eta}_{\rho})\cdot {H}_{x}(x,{\beta}_{m})\cdot {H}_{y}(y,{\gamma}_{n})}\right\}}$ | $T(x,y,z,t)={T}_{0}+\frac{1}{k}{\displaystyle \sum _{p=1}^{\infty}\left\{{H}_{a}(z,{\eta}_{\rho})+{E}_{0}\cdot (1-R)\cdot {\mu}_{a}\cdot {H}_{z}(z,{\eta}_{\rho})\cdot {\displaystyle \sum _{n=0}^{\infty}}{\displaystyle \sum _{m=0}^{\infty}{H}_{t}(t,{\beta}_{m},{\gamma}_{n},{\eta}_{\rho})\cdot {H}_{xy}(x,y,{\beta}_{m},{\gamma}_{n})}\right\}}$ | |
Time dependence H_{t}(t,β_{m},γ_{n},η_{p}) | $t\le {t}_{p}:\underset{(m,n,p)\ne (0,0,0)}{{H}_{t}(t<{t}_{p},{\beta}_{m},{\gamma}_{n},{\eta}_{p})}=\frac{1-{e}^{-\alpha ({\beta}_{m}^{2}+{\gamma}_{n}^{2}+{\eta}_{p}^{2})t}}{{\beta}_{m}^{2}+{\gamma}_{n}^{2}+{\eta}_{p}^{2}},$ and H_{t}(t<t_{p},β_{0},γ_{0},η_{0})=α⋅t | $$t>{t}_{p}:\underset{(m,n,p)\ne (0,0,0)}{{H}_{t}(t>{t}_{p},{\beta}_{m},{\gamma}_{n},{\eta}_{p})}={e}^{-\alpha ({\beta}_{m}^{2}+{\gamma}_{n}^{2}+{\eta}_{p}^{2})t\cdot}\frac{{e}^{\alpha ({\beta}_{m}^{2}+{\gamma}_{n}^{2}+{\eta}_{p}^{2}){t}_{p}-1}}{{\beta}_{m}^{2}+{\gamma}_{n}^{2}+{\eta}_{p}^{2}},$$ and H_{t}(t>t_{p},β_{0},γ_{0},η_{0})=α⋅t_{p} |
Ambient temperature H_{a}(z,η_{p}) | ${H}_{a}\left(z,{\eta}_{\rho}\right)={h}_{1}\cdot \left({T}_{a}-{T}_{0}\right)\cdot \frac{2.\left({\eta}_{p}\mathrm{cos}{\eta}_{p}z+\frac{{h}_{1}}{k}\mathrm{sin}{\eta}_{p}z\right)}{\left({\eta}_{p}^{2}+\frac{{h}_{1}^{2}}{{k}^{2}}\right)\left(c+\frac{{h}_{2}\cdot k}{{\eta}_{p}^{2}{k}^{2}+{h}_{2}^{2}}\right)+\frac{{h}_{1}}{k}}\cdot \frac{1-{e}^{-\alpha \cdot {\eta}_{p}^{2}\cdot t}}{{\eta}_{p}}$ | |
x−y
dependence^{*} H_{xy}(x,y,β_{m},γ_{n}) H_{x}(x,β_{m}),H_{y}(y,γ_{n}) | $$\begin{array}{l}{H}_{x}(X,{\beta}_{2m})=\frac{\epsilon}{a}\cdot \sqrt{2}\cdot w\cdot \mathrm{cos}{\beta}_{2m}X\cdot {(-1)}^{m}\cdot I\left(2m,\frac{a}{\sqrt{2}\cdot w}\right)\\ {H}_{x}(X,{\beta}_{2m+1})=0\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\epsilon =1\text{\hspace{0.17em}}\text{if}\text{\hspace{0.17em}}\text{m=0,}\text{\hspace{0.17em}}\epsilon \text{=2}\text{\hspace{0.17em}}\text{if}\text{\hspace{0.17em}}\text{m}\ne \text{0}\\ {H}_{y}(Y,{\gamma}_{2n})=\frac{\epsilon}{b}\cdot \sqrt{2}\cdot w\cdot \mathrm{cos}{\gamma}_{2n}Y\cdot {(-1)}^{n}\cdot I\left(2n,\frac{b}{\sqrt{2}\cdot w}\right)\\ {H}_{y}(Y,{\gamma}_{2n+1})=0\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\epsilon =1\text{\hspace{0.17em}}\text{if}\text{\hspace{0.17em}}\text{n=0,}\text{\hspace{0.17em}}\epsilon \text{=2}\text{\hspace{0.17em}}\text{if}\text{\hspace{0.17em}}\text{n}\ne \text{0}\end{array}$$ | $\underset{m\ne 0,n\ne 0}{{H}_{xy}(X,Y,{\beta}_{2m},{\gamma}_{2n})}=\frac{8\pi w}{{a}^{2}\cdot {({\beta}_{2m}^{2}+{\gamma}_{2n}^{2})}^{1/2}}\cdot {(-1)}^{m+n}\mathrm{cos}{\beta}_{2m}X\cdot \mathrm{cos}{\gamma}_{2n}Y\cdot {\text{J}}_{1}[{({\beta}_{2m}^{2}+{\gamma}_{2n}^{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 |
H_{y}(y,γ_{2n+1})=0 | ||
z dependence^{†} | z=0: k∂T/∂z=h_{1}T z=c: k∂T/∂z=−h_{2}T | |
${H}_{z}\left(z,{\eta}_{\rho}\right)=\frac{2.\left({\eta}_{p}\mathrm{cos}{\eta}_{p}z+\frac{{h}_{1}}{k}\mathrm{sin}{\eta}_{p}z\right)}{\left({\eta}_{p}^{2}+\frac{{h}_{1}^{2}}{{k}^{2}}\right)\left(c+\frac{{h}_{2}\cdot k}{{\eta}_{p}^{2}{k}^{2}+{h}_{2}^{2}}\right)+\frac{{h}_{1}}{k}}\cdot \frac{{\eta}_{\rho}\cdot \left({\mu}_{a}+\frac{{h}_{1}}{k}\right)\cdot \left(1-{e}^{-{\mu}_{a}c}\cdot \mathrm{cos}{\eta}_{\rho}c\right)+{e}^{-{\mu}_{a}c}\cdot \mathrm{sin}{\eta}_{\rho}c\cdot \left({\eta}_{p}^{2}-{\mu}_{a}\frac{{h}_{1}}{k}\right)}{{\mu}_{a}^{2}+{\eta}_{p}^{2}}$ | ||
z=0: ∂T/∂z=0 z=c: ∂T/∂z=0$\text{ifp}\ne \text{0:}{{\rm H}}_{z}(z,{\eta}_{\rho})=\frac{2}{c}\cdot \mathrm{cos}{\eta}_{\rho}z\cdot \frac{{\mu}_{a}}{{\mu}_{a}^{2}+{\eta}_{\rho}^{2}}\left[1-{e}^{-{\mu}_{{a}^{c}}}\cdot \mathrm{cos}{\eta}_{\rho}c\right]$$\text{ifp=0:}{{\rm H}}_{z}(z,{\eta}_{0})=\frac{1}{c\cdot {\mu}_{a}}\left[1-{e}^{-{\mu}_{{a}^{c}}}\right]$ | ||
^{*}J_{1}(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)={\displaystyle {\int}_{u=0}^{L}\mathrm{cos}\left(q\pi \frac{u}{2L}\right)\cdot {e}^{-{u}^{2}}\cdot du\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_{1} or h_{2}) . |
(3)
$$\begin{array}{c}{\mathit{H}}_{\mathit{xyz}}(\mathit{x},\mathit{y},z,{\beta}_{\mathit{m}},{\gamma}_{\mathit{n}},{\eta}_{\mathit{p}})=\frac{\mathit{X}({\beta}_{\mathit{m}},\mathit{x})\mathit{Y}({\gamma}_{\mathit{n}},\mathit{y})\mathit{Z}({\eta}_{\mathit{p}},\mathit{z})}{\mathit{N}({\beta}_{\mathit{m}})\mathit{N}({\gamma}_{\mathit{n}})\mathit{N}({\eta}_{\mathit{p}})}\text{}\cdot {\displaystyle {\int}_{{x}^{\prime}=0}^{a}{\displaystyle {\int}_{{y}^{\prime}=0}^{b}{\displaystyle {\int}_{{z}^{\prime}=0}^{c}}}}\hfill \\ \phantom{\rule{1.56em}{0ex}}\mathit{X}({\beta}_{\mathit{m}},{\mathit{x}}^{\prime})\mathit{Y}({\gamma}_{\mathit{n}},{\mathit{y}}^{\prime})\mathit{Z}({\eta}_{\mathit{p}},{\mathit{z}}^{\prime})\cdot \mathit{g}({\mathit{x}}^{\prime},{y}^{\prime},{z}^{\prime})\times d{x}^{\prime}{\mathit{dy}}^{\prime}{\mathit{dz}}^{\prime}\cdot \hfill \end{array}$$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(β_{m}x) |
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(γ_{n}y) |
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_{1}(T−T_{a}) | $Z\left({\eta}_{\rho},z\right)={\eta}_{\rho}\xb7\mathrm{cos}\left({\eta}_{\rho}z\right)+\frac{{h}_{1}}{k}\xb7\mathrm{sin}\left({\eta}_{\rho}z\right)$ |
z=c: k∂T/∂z=−h_{2}T | $N\left({\eta}_{\rho}\right)=\frac{1}{2}\xb7\left[\frac{{h}_{1}}{k}+\left({\eta}_{\rho}^{2}+\frac{{h}_{1}^{2}}{{k}^{2}}\right)\xb7\left(c+\frac{k\xb7{h}_{2}}{{\eta}_{\rho}^{2}{k}^{2}+{h}_{2}^{2}}\right)\right]$ | |
The eigenvalues are all positive roots of the equation: | ||
$\mathrm{tan}({\eta}_{\rho}c)=\frac{{\eta}_{\rho}\xb7k\xb7\left({h}_{1}+{h}_{2}\right)}{{\eta}_{\rho}^{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:
(4)
$$g(x,y,z)=(1-\mathit{R})\cdot {\mu}_{\mathit{a}}\cdot {\mathit{E}}_{\mathit{i}}(\mathit{x},\mathit{y})\cdot \text{exp}(-{\mu}_{\mathit{a}}\mathit{z}),$$(5)
$${\mathit{H}}_{\mathit{xyz}}(\mathit{x},\mathit{y},z,{\beta}_{\mathit{m}},{\gamma}_{\mathit{n}},{\eta}_{\mathit{p}})=(1-\mathit{R})\cdot {\mu}_{\mathit{a}}\cdot {\mathit{H}}_{\mathit{z}}(\mathit{z},{\eta}_{\mathit{p}})\cdot {\mathit{H}}_{\mathit{xy}}(\mathit{x},\mathit{y},{\beta}_{\mathit{m}},{\gamma}_{\mathit{n}}),$$(6)
$${\mathit{H}}_{\mathit{z}}(\mathit{z},{\eta}_{\mathit{p}})=\frac{\mathit{Z}({\eta}_{\mathit{p}},\mathit{z})}{\mathit{N}({\eta}_{\mathit{p}})}\text{}\cdot {\int}_{{z}^{\prime}=0}^{c}\mathit{Z}({\eta}_{\mathit{p}},{\mathit{z}}^{\prime})\cdot \text{exp}(-{\mu}_{\mathit{a}}{\mathit{z}}^{\prime})d{z}^{\prime},$$(7)
$${\mathit{H}}_{\mathit{xy}}(\mathit{x},\mathit{y},{\beta}_{\mathit{m}},{\gamma}_{\mathit{n}})=\frac{\mathit{X}({\beta}_{\mathit{m}},\mathit{x})\cdot \mathit{Y}({\gamma}_{\mathit{n}},\mathit{y})}{\mathit{N}({\beta}_{\mathit{m}})\mathit{N}({\gamma}_{\mathit{n}})}\text{}\cdot {\displaystyle {\int}_{{x}^{\prime}=0}^{a}{\displaystyle {\int}_{{y}^{\prime}=0}^{b}}}\mathit{X}({\beta}_{\mathit{m}},{\mathit{x}}^{\prime})\phantom{\rule{0.56em}{0ex}}\cdot \mathit{Y}({\gamma}_{\mathit{n}},{\mathit{y}}^{\prime})\cdot {\mathit{E}}_{\mathit{i}}({\mathit{x}}^{\prime},{y}^{\prime})d{x}^{\prime}{\mathit{dy}}^{\prime}\cdot $$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:
(8)
$${\mathit{E}}_{\mathit{i}}(\mathit{x},\mathit{y})={\mathit{E}}_{0}\cdot {\mathit{e}}^{-2{(\mathit{x}-\mathit{a}/2)}^{2}/{w}^{2}}\cdot {\mathit{e}}^{-2{(\mathit{y}-\mathit{b}/2)}^{2}/{w}^{2}},$$(9)
$${\mathit{H}}_{\mathit{xy}}(\mathit{x},\mathit{y},{\beta}_{\mathit{m}},{\gamma}_{\mathit{n}})={\mathit{E}}_{0}\cdot {\mathit{H}}_{\mathit{x}}(\mathit{x},{\beta}_{\mathit{m}})\cdot {\mathit{H}}_{\mathit{y}}(\mathit{y},{\gamma}_{\mathit{n}})\cdot $$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:
(10)
$${\mathit{E}}_{\mathit{i}}(\mathit{x},\mathit{y})={\mathit{E}}_{0}\cdot \mathit{u}{\displaystyle [{({w}^{2}-{y}^{2})}^{1/2},\text{}\frac{a}{2},\mathit{x}]}\cdot \mathit{u}{\displaystyle (w,\text{}\frac{b}{2},\mathit{y})},$$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_{0}). 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^{−3} 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:^{21}
(11)
$${\int}_{x=0}^{\infty}\text{}\text{cos}\left(\mathit{q}\pi \text{}\frac{x}{2L}\right)\cdot {\mathit{e}}^{-{x}^{2}}\mathit{dx}=\frac{\sqrt{\pi}}{2}\text{}{\mathit{e}}^{-{(\mathit{q}\pi /4L)}^{2}}\cdot $$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^{−2} K^{−1}.^{22} 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^{22} modeled the effect of aqueous flow by assuming that the posterior corneal surface is subject to free convection with constants ranging from h_{2}=20 to 1000 W m^{−2} °C^{−1}. 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),^{23} earlier thermal models assumed that the convective effect of aqueous flow is negligible.^{24} 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.^{24}
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_{1} was fixed (h_{1}=20 W m ^{−2} K ^{−1}) and h_{2} was varied from 20 to 2000 W m^{−2} K^{−1}. Second, h_{2} was fixed (h_{2}=1000 W m ^{−2} K ^{−1}) and h_{1} was varied from 20 to 500 W m^{−2} K^{−1}. We also calculated the temperature with h_{1}=20 W m ^{−2} K ^{−1} 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 ^{−1}), 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_{1} do not have a significant effect on the depth and time dependence of the corneal temperature during the treatment, except for h_{1}=500 W m ^{−2} K ^{−1}, which is an extreme value (Fig. 2). On the other hand, variations of h_{2} 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_{2}=1000 W m ^{−2} K ^{−1}. These results are further discussed in Sec. 5.
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_{1} and (T_{a}−T_{0}). 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_{2}=1000 W m ^{−2} K ^{−1} with values of h_{1} ranging from 20 to 500 W m^{−2} K^{−1} 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_{1}=20, 100, and 500 W m^{−2} K^{−1}, respectively. For T_{a}=35 °C, the anterior corneal surface temperature after 1.4 s for the same values of h_{1} 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_{1} ranging from 20 to 100 W m^{−2} K^{−1}. The effect becomes significant (>4 °C) only in the extreme case when h_{1}=500 W m ^{−2} K ^{−1}.
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_{1}=20 W m ^{−2} K ^{−1} and h_{2}=1000 W m ^{−2} K ^{−1}, and T_{a}=20 °C (Fig. 4). Figure 4 also shows the temperature at different depths as a function of time during the treatment.
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.^{20}
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.^{25} ^{26} Ith, Frenz, and Weber^{11} 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.;^{15} 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_{1} 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_{2}=1000 W cm ^{−2} K ^{−1}. 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_{2} 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_{2}=1000 W cm ^{−2} K ^{−1}.
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^{27} and local epithelial and anterior stromal scarring^{28} 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.^{29} ^{30} 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.