Introduction

Since its inception, controversy has surrounded the exact etiology of thermal damage to the saphenous vein [1-5]. There are proponents of both conduction energy (laser tip contact) and convective energy (steam bubbles) as being the etiology of resultant tissue damage. This controversy has been debated not only from the scientific aspect, but in the legal arena as well.

As reported in previous publications [1, 5], the hallmark of acute damage at the time of thermal ablation is complete loss of endothelium. This finding is present when laser firing at 12-15 watts is done at intervals of 3-4 mm apart. If it is assumed that the laser point contacts the surface of the vein, which occurs when a conventional bare tip fiber is used, then the question arises of whether the contact and the resultant conduction by tissue is responsible for the histologic findings. The absence of endothelium immediately following laser thermal ablation is a constant finding regardless of wavelength and pulse duration. Based on our histologic observations it is our hypothesis that heat conduction would not be the primary etiology of these observed luminal changes. The purpose of this article is to present the results of testing our hypothesis by simulations based on the construction of a relevant mathematical model which will accurately calculate both temperature and distance traveled of energy from the laser application site.

Construction of the model

The key parameter of the model is the initial temperature to which the laser probe elevates the affected region of the vein wall. In the model we will keep this temperature as a free parameter, but when we present results of simulations based on the model, we will assume that at the time of laser firing the initial temperature at the laser tip is 100 C. This assumption is reasonable since (1) steam bubbles occur at this point, and (2) we have measured the temperature of these bubbles impinging on the vein wall. We have run simulations at other initial temperatures, such as 700 C and the associated results do not significantly change our final conclusions.

The first step in building a valid mathematical model is the enunciation of the assumptions on which the model rests. We base our model on five postulates derived from both experience and theory.

Postulate 1: We assume that the probe tip contacts the partially collapsed vein wall in a rectangular region. In this region the transfer of a short burst of energy raises the temperature to an initial value T0.

Discussion: The rectangular shape is chosen so as to model the region where the probe tip contacts the partially collapsed vein wall. Our analysis will indicate that the time evolution of the heat conduction smooths out the boundary of the initial region so that the assumption of a rectangular shape is not of particular importance in determining the region of cell death. As indicated above, we will present the results of simulations for a value of T₀ of 100 C.

Postulate 2: We assume that no transport of mass occurs along the vein wall during the process, so that the only transport involved is the transport of energy in the form of heat.

Postulate 3: We assume that the vein wall is a right circular cylinder of infinite length.

Discussion: That the saphenous vein is a uniform right circular cylinder is a simplifying assumption that models the actual vein as it appears in ultrasounds of the procedure and in clinical experience. The tumescent solution outside of the vein certainly exerts no tensile forces on the vein wall, and any compression it might provide has minimal effect on the surface area. The assumption of infinite length merely indicates that the regions of laser application are far enough from the ends of the vein that the longitudinal boundaries play no role in the local heat conduction process.

Postulate 4: We assume that after the initial thermal contact the vein wall remains a closed system.

Discussion: This postulate, along with Postulate 2, restricts the phenomenon at hand to one of heat conduction along the vein wall. The assumption that heat conduction is the primary mechanism causing endothelial necrosis thus requires Postulates 2 and 4 as its underlying support. From this point of view the effect of energy transfer through the blood and into the vein wall is a separate problem. It is true that dissipation of energy into the ambient anesthetic and tissue occurs. Experimental results [6] indicate that such heat loss is minimal. In any case assumption of Postulate 4 ensures that whatever region of cell death we obtain in this model will be an upper bound for the actual region. Any model incorporating other modes of energy loss will certainly produce smaller regions of cell death.

Postulate 5: We assume Fourier's law of heat conduction: Heat flows from regions of higher temperature to lower temperature along the thermal gradient -T, where

T = T(x, y, z, t) is the temperature at the point with coordinates x = (x, y, z) at time t.

Discussion: Postulates 4 and 5 imply that heat flows along the vein wall in such a way that energy is conserved. Standard energy conservation arguments indicate that the local behavior of the heat flow is described by the heat conduction equation [7].

Figure 1

Here ΔT is the Laplacian of T and the constant k is the thermal diffusivity of the vein wall. We take for the value of thermal diffusivity that of human myocardium tissue: k = 1.289 x 10-7 m2/sec.

The heat conduction problem: the planar approximation

Consider a right-handed coordinate system with the z-axis along the central longitudinal axis of the vein and the x-axis through the center of the rectangle initially raised to a temperature T₀ on the vein wall. In cylindrical coordinates equation (1) becomes

Figure 2

Here φ is the azimuthal angle measured counterclockwise from the x-axis and ρ is the radial distance out from the origin. Along the vein wall ρ is a constant, R, so that on the cylindrical surface the equation becomes

Figure 3

We take the value of the radius of the vein as ρ = R = 4 mm.

We may shift to a local coordinate system on the cylinder and centered on the rectangle initially at temperature T₀. The local coordinates then become the circumferential distance

s = Rφ and the longitudinal distance z as shown in Figure 1. Also shown is the initial condition rectangle in red. The length 2a is the diameter of the probe tip, while 2b is the length of the probe tip. Typically a = .25 mm and b = 2 mm.

Figure 4

Figure 1 Local Coordinates on the Vein Wall

In these local coordinates equation (3) assumes the form:

Figure 5

The initial condition problem associated with equation (4) has circumferential boundary conditions which are periodic with period 2pR. Since equation (4) is a diffusion-type equation, we expect the temperature distribution to spread out from its initial rectangular region. When it spreads as far as ±pR ~ 12.57 mm there will therefore be an overlapping of the left and right traveling “tails” of the solution. The linearity of equation (4) ensures that the overlap solution is the superposition of the tails. On the other hand, for the time periods of importance overlap effects are very small and, as we shall see, the spreading heat distribution is essentially contained within the region - pR < s < + pR. Provided we keep this domain and time periods in view, we may replace this domain with -∞ < s < +∞; consequently, the relevant planar approximation problem is:

Figure 6

Here Θ(x) is the Heaviside function.

Solution of the planar approximation problem

Our method of solution will be to construct the two-dimensional Green’s function for equation (4) and then to integrate the product of the Green’s function and the initial condition.

[See for example 8, 9]

The relevant Green's function problem is to solve equation (4) subject to the initial condition, T(s, 0) = (s – s ₀), where s = (s, z) is the local coordinate position vector and

(s – s ₀) is the Dirac delta function centered at s ₀. This problem may be solved by using the two-dimensional Fourier transform of T(s, t):

Figure 7

Taking the Fourier transform of (4) yields an ordinary differential equation:

Figure 8

Here a² = a ∙ a. The accompanying initial condition becomes:

Figure 9

Solving for τ gives

Figure 10

Taking the inverse transform yields the Green's function:

Figure 11

The initial condition problem has solution:

Figure 12

Upon substitution of the initial condition this becomes:

Figure 13

The solution we have obtained is the product of two independent Gaussian distributions spreading along the vein wall. One Gaussian spreads longitudinally while the other spreads azimuthally. The formula obtained for T(s, z, t) permits the simulation of the conduction of heat away from the laser tip along the vein wall. Figures 2 and 3 show the results of such a simulation for an initial temperature of T₀ = 100 C and at times t = .63 sec and t = 5.00 sec. In each figure the region of possible cell death is indicated in white. We estimate that, in each 15 mm longitudinal segment of the saphenous vein, not more than 14 per cent of the vein wall experiences temperatures above the cytotoxic threshhold. In the figures one can also observe the smoothing out of the boundary of the initial condition rectangle. Such smoothing of boundaries under heat conduction is well-known and indicates that in our Postulate 1 other assumptions on the shape of the initial temperature distribution would give results similar to ours.

Figure 14

Figure 2 Temperature distribution as seen looking down the T-axis at t = .63 sec. s denotes circumferential distance around the vein wall, while z denotes longitudinal distance along the vein wall. The region above the cytotoxic threshold, taken as 45 C, is indicated in white.

Figure 15

Figure 3 Temperature distribution as seen looking down the T-axis at t = 5.00 sec. s denotes circumferential distance around the vein wall, while z denotes longitudinal distance along the vein wall. The region above the cytotoxic threshold, taken as 45 C, is indicated in white.

Figure 4 displays the azimuthal Gaussian distribution for the temperature ratio, R = T/T₀, shortly after the burst of energy from the laser probe at t = 1.25 sec. Figure 5 shows the same distribution at a later time t = 10.0 sec. Since overlapping of the tails of the temperature distribution occurs for s = ± .000785 m, these graphs indicate that for small times there is no appreciable overlap effect. At large times, however, the overlap effect may be as much as twenty per cent of the maximum temperature at the corresponding time. In the next section we shall construct a solution that will be valid at these longer times.

Figure 16

Figure 4 Azimuthal spreading of the distribution of the temperature ratio R as a function of circumferential distance s at t = 1.25 sec. The horizontal scale gives the distance s along the cell wall in meters.

Figure 17

Figure 5 Azimuthal spreading of the distribution of the temperature ratio R at t = 10.0 sec. The horizontal scale gives the distance s along the cell wall in meters.

The initial condition heat conduction problem with periodic boundary condition

The problem at hand is to incorporate the circumferential boundary conditions into the initial condition problem; consequently, we seek a solution which is periodic in s with period 2πR. In this case the problem may be formulated as:

Figure 18

The sequence of functions

Figure 19

for n belonging to , forms an orthonormal set on the interval [-πR, + πR]. We expand T(t, s, z) in terms of this basis:

Figure 20

This Fourier series expansion for T(t, s, z) manifestly satisfies the periodic boundary condition.

Substitution of the expansion (7) for T(t, s, z) into the partial differential equation (6) yields

Figure 21

The linear independence of the functions {e_n (s)} gives for each n:

Figure 22

The initial condition dependent on z is specialized from f(s, z) and may be written:

Figure 23

where is a constant to be determined. Let σ_n (t, ω) be the Fourier transform of τ_n (t, z) with respect to z. The Fourier transform of equation (8) yields an ordinary differential equation for σ_n which has solution:

Figure 24

Here σ_0n is a constant determined by the initial condition for τ_n (t, z).

The Green’s function solution for τ_n (t, z) may be obtained by assuming that for each n τ_n (0,z) = δ(z – ζ). σ_n (t, ω) accordingly becomes

Figure 25

Taking the inverse Fourier transform yields the Green’s function solution:

Figure 26

The solution for τ_n (t, z) may be obtained by integrating Г_n (z|ζ, t) against the initial condition (9). Substituting the resulting expression for τ_n (t, z) back into equation (7) produces the general solution for T(t, s, z) with t > 0:

Figure 27

The coefficients τ_0n may be determined by use of the azimuthal initial condition which requires

Figure 28

The orthogonality of the basis functions {e_n (s)} permits calculation of the τ_0n :

Figure 29

and for n ≠ 0:

Figure 30

Substitution into equation (10) provides the desired solution for the temperature function:

{image:31}

Figures 6 and 7 show the results of simulations of this solution for an initial temperature of T₀ = 100 C and at times t = .63 sec and t = 5.00 sec. In both simulations the series in s are expanded to include the first 1500 terms. As before the region of possible cell death is indicated in white. The regions above the cytotoxic limit are smaller than the regions at the same times in the simulation results given in the previous section. The similarity of the results thus indicate the validity of the planar approximation approach.

{image:32}

{image:33}

Discussion

Ultrasounds taken during endovenous ablation consistently show the formation of steam bubbles as demonstrated by Proebstle et al. [10] and Yilmaz et al. [3]. Since the erythrocyte is 80% water and considering the affinity of these wavelengths for water, it is reasonable to assume that these bubbles represent the formation of steam.

There are three types of energy that can account for the observed tissue damage. These are conductive (heat transfer into the target tissue), convective (heat transfer into the liquid), and radiant (light absorption). Much of the radiant energy is rapidly absorbed by the heme protein and water. The model we have examined in this paper is based on the assumption that all of the energy from the laser probe results in conductive energy occurring at the point of laser contact. This assumption has provided a way of testing the most extreme case of conductive damage to the endothelium. In view of the results of our simulations it is reasonable to conclude that large-scale endothelial necrosis is a result of energy transfer in the convective and radiant channels.

Based on our mathematical model, tissue conduction by laser application only accounts for about 15% of cellular damage when temperatures of 45 C are used as the initial point of irreversible cell damage. This point is well below the 50 C needed for cell death. When 50 C is considered, the area of destruction would be even smaller.

Conclusion

We conclude from the preceding discussion that massive cell death of the endothelium cannot be due solely to conduction along the inner surface of the vein wall. Laser tip contact is not necessary for thermal ablation to be effective. Experimental evidence [5] indicates that the fluid in the vein flows in a highly turbulent manner and that portions of the fluid undergo a liquid-gas phase transition. The ordinary diffusion equation is no longer valid in such a regime. Modelling heat transfer in this regime remains an open question.