Simulation study of the cooling effect of blood vessels and blood coagulation in hepatic radio-frequency ablation

Abstract

Purpose

Computer simulations of hepatic radio-frequency ablation (RFA) were performed to: (i) determine the dependence of the vessel wall heat transfer coefficient on geometrical parameters; (ii) study the conditions required for the occurrence of the directional effect of blood; and (iii) classify blood vessels according to their effect on the thermal lesion while considering blood coagulation. The information thus obtained supports the development of a multi-scale bio-heat model tailored for more accurate prediction of hepatic RFA outcomes in the vicinity of blood vessels.

Materials and methods

The simulation geometry consisted of healthy tissue, tumor tissue, a mono-polar RF-needle, and a single cylindrical blood vessel. The geometrical parameters of interest were the RF-needle active length and those describing blood vessel configuration. A simple, novel method to incorporate the effects of blood coagulation into the simulation was developed and tested.

Results

A closed form expression giving the dependence of the vessel wall heat transfer coefficient on geometrical parameters was obtained. Directional effects on the thermal lesion were found to occur for blood vessel radii between 0.4 mm and 0.5 mm. Below 0.4 mm blood coagulation blocked the flow.

Conclusions

The closed form expression for the heat transfer coefficient can be used in models of RFA to speed up computation. The conditions on vessel radii required for the occurrence of directional effects on the thermal lesion were determined. These conditions allow the classification of blood vessels. Different approximations to the thermal equation can thus be used for these vessel classes.

Keywords:

• Radio-frequency ablation
• modeling
• directional effect
• blood coagulation
• Nusselt number

Introduction

In recent years, minimally invasive techniques such as radio-frequency ablation (RFA), microwave ablation, and cryotherapy have become popular in clinical practice for the treatment of liver cancer. In RFA, in particular, an alternating current flowing between the active region of an RF-needle inserted in the tumor and a ground pad is responsible for localized tissue ablation through heating.

During RFA, the challenge is to ablate the complete tumor without harming the adjacent healthy or vital structures [1]. The final lesion is affected by fine-scale perfusion in the liver and by large vessels which have a strong localized cooling effect. The micro-vasculature in the ablated tissue is destroyed and blood coagulation occurs in the vessels neighboring the targeted area with diameter smaller than 3 mm [2]. Moreover, imaging techniques currently used in clinical practice (e.g., CT and MRI) cannot detect small vessels in the liver (dia < 1 mm).

Mathematical modeling of the vessel effect could predict the temperature distribution in the liver and hence aid treatment planning. Such models must be fast yet accurate to be useful in practice. Chen and Saidel [3] used the transient Pennes bio-heat equation with a 1D advection-diffusion approximation for the large blood vessel to study the effect of the vessel radius, distance between the RF-needle and vessel, and the angle between the RF-needle and vessel on the vessel cooling effect. They assumed the blood vessel wall heat transfer coefficient to be a constant independent of vessel radius. Huang used the steady Pennes bio-heat equation with a steady 1D convection-diffusion approximation for the large blood vessel for a similar study, using a constant blood vessel wall heat transfer coefficient. Shao et al. [4] studied the cooling effect of complicated vessel trees. They examined the impact of varying the average distance of the needle from vascular structures, the total vascular volume, and the number of branches in the vascular structure on the cooling effect. They found that the cooling effect increases with decreasing distance, increasing vascular volume, and increasing number of blood vessel branches. Existing works on multi-scale modeling have not considered the possibility of directional effects of isolated vessels [1,5,6]. While low blood flow-rates were found to cause a directional effect on the thermal lesion [7], no detailed study has been performed yet on the influence of the vessel radius. A study by dos Santos et al. [8] found the vessel wall heat transfer coefficient to be independent of temperature, but it did not comment on the dependence of this coefficient on geometrical parameters. To the best of our knowledge, simulation studies on RFA have not considered the possibility and effect of blood coagulation in vessels of diameter larger than 1 mm.

This work aims to address these shortcomings in the state of the art of RFA modeling, particularly focusing on (i) understanding the dependence of the heat transfer coefficient on geometrical parameters, (ii) determining the conditions at which directional effects and blood coagulation occur, and finally (iii) understanding the dependence of the thermal lesion on the blood vessel diameter.

Methods

Geometry

The computational geometry (Figure 1(a)) consisted of a cylindrical healthy tissue region of both diameter and height 12 cm containing a tumor (radius 0.5 cm) positioned in the center, a simplified UniBlate^® RF-needle (from Angiodynamics, New York, NY, USA), and a blood vessel of variable radius r. Both height and diameter of the cylindrical tissue were chosen to be equal to 12 cm, similarly to previous works found in literature [9–11] and considered enough to avoid boundary condition effects. To reduce computational complexity, the internal structure of the RF-needle was assumed to be homogeneous with averaged material properties (density, specific heat capacity, thermal conductivity, and electrical conductivity). The exact outer dimensions of the RF-needle (diameter 1.5 mm) were used. The blood vessel position with respect to the needle was varied according to three parameters: d₁, d₂, and α (Figure 1(b)). This parametrization took advantage of the cylindrical symmetry of the domain. During the RFA procedure, the current passed through the active surface of the needle with variable length λ (Figure 1(c)), and toward the ground pad. Simulations of RFA were performed using different values of the variable parameters described above, and were grouped into two sets (SET A and SET B) based on the way the parameter values varied.

Figure 1. (a) Simulation geometry. The bottom surface of the tissue cylinder was electrically grounded. (b) Views of the simulation geometry showing the parameters α, d₁, and d₂. Note that the distance d₁ is negative if the blood vessel axis is located below the center of the RF-needle active region in the side view shown and vice versa. (c) Enlarged view of the RF-needle geometry showing the definition of parameter λ. The axial length from the top surface to the active surface is 50 mm.

Governing equations

Bio-heat models for hepatic RFA require separate treatment of fine-scale perfusion and large blood vessels. A number of models describing the fine-scale perfusion effects have been proposed, but their experimental validation is currently difficult [12–14]. It was decided to omit the fine-scale perfusion from the governing equations, since the focus of the present work is the effect of large blood vessels. The governing equations for solid tissue, blood, and RF-needle were given by (1) $$\begin{array}{ll}{\rho }_t{C}_t\frac{\partial T_t}{\partial t}\hfill & =\nabla ·k_t\nabla T_t+{\sigma }_t|\nabla {\varphi }_t{|}^2,\hfill \\ {\rho }_b{C}_b\left(\frac{\partial T_b}{\partial t}+{\overrightarrow{u}}_b·\nabla T_b\right)\hfill & =\nabla ·k_b\nabla T_b+{\sigma }_b|\nabla {\varphi }_b{|}^2,\hfill \\ {\rho }_n{C}_n\frac{\partial T_n}{\partial t}\hfill & =\nabla ·k_n\nabla T_n+{\sigma }_n|\nabla {\varphi }_n{|}^2,\hfill \end{array}$$(1) where the subscripts t, b, and n stand for tissue, blood, and RF-needle respectively. The constants ρ_i, C_i, T_i, k_i, σ_i, and ϕ_i are respectively the density, specific heat, temperature, thermal conductivity, electrical conductivity, and electric potential for $i\in \{t,b,n\}.$ Their values are listed in Table 1. The thermal and electrical conductivities of tissue and blood were considered temperature dependent. All specific heats and densities were assumed constant. A convection boundary condition, $-k_t\nabla ·T_t=6(T_t-T_{\text{ref}})$ W m^–2, was applied to the tissue boundary, where $T_{\text{ref}}=310.15$ K. The blood inlet and outlet boundaries were kept at 310.15 K. The temperature and heat-flux were assumed continuous across all interfaces. Blood was assumed to be a Newtonian fluid and shear effects were neglected. Blood flow was assumed to be steady and laminar [15].

Table 1. Table showing material property values used.

Region	Material property	Value
Tissue	ρ kg m^–3	1079 [22]
	C J kg^–1 K^–1	3540 [22]
	$k(T{}^{°}\mathrm{C})$ W m^–1 K^–1	$4.69e-1+1.16e-3\times T$ [23]
	Healthy Tissue
	$\sigma (T{}^{°}\mathrm{C})$ S m^–1	$-1.26e-7\times T^4+2.57e-5\times T^3$
	Healthy tissue	$-1.84e-3\times T^2+6.35e-2\times T$
		$-4.79e-1$ [24]
	$k(T{}^{°}\mathrm{C},r_{\text{tum}}\text{mm})$	$0.561:r_{\text{tum}}<3$ [25]
	W m^–1 K^–1	$0.419+3.91e-3\times T:r_{\text{tum}}\ge 3$ [23]
	Tumor tissue
	$\sigma (T{}^{°}\mathrm{C})$ S m^–1	$-1.26e-7\times T^4+2.57e-5\times T^3$
	Tumor tissue	$-1.84e-3\times T^2+6.35e-2\times T$
		–0.378 [26]
	A s^–1	$3.1e98$
	$\Delta E_a$ J	$6.285e5$
	R J mol^–1 K^–1	8.314
Blood	ρ kg m^–3	1050 [22]
	C J kg^–1 K^–1	3617 [22]
	$k(T{}^{°}\mathrm{C})$ W m^–1 K^–1	$0.552-2.26e-3\times (37-T)$ [27]
	$\sigma (T{}^{°}\mathrm{C})$ S m^–1	$100\times {(326.0-5.63\times T)}^{-1}$ [28]
	Healthy blood
	σ S m^–1 Coagulated blood	0.3114 [29]
	A s^–1	$7.66e66$
	$\Delta E_a$ J	$4.225e5$
	R J mol^–1 K^–1	8.314
RF	ρ kg m^–	1
	C J kg^–1 K^–1	1000
Needle	k W m^–1 K^–1	0.02
	σ S m^–1	1e7
Phantom	ρ kg m^–3	1033 [30]
	C J kg^–1 K^–1	3939 [30]
	k W m${}^{}$ K^–1	0.59 [30]
	σ S m^–1	1

Liver tumors are known to have different thermal conductivities in their core and periphery. The tumor core was assumed to be the region $r_{\text{tum}}<3$ mm.

In Equation (1)(1) $$\begin{array}{ll}{\rho }_t{C}_t\frac{\partial T_t}{\partial t}\hfill & =\nabla ·k_t\nabla T_t+{\sigma }_t|\nabla {\varphi }_t{|}^2,\hfill \\ {\rho }_b{C}_b\left(\frac{\partial T_b}{\partial t}+{\overrightarrow{u}}_b·\nabla T_b\right)\hfill & =\nabla ·k_b\nabla T_b+{\sigma }_b|\nabla {\varphi }_b{|}^2,\hfill \\ {\rho }_n{C}_n\frac{\partial T_n}{\partial t}\hfill & =\nabla ·k_n\nabla T_n+{\sigma }_n|\nabla {\varphi }_n{|}^2,\hfill \end{array}$$(1) the term ${\rho }_b{C}_b{\overrightarrow{u}}_b·\nabla T_b$ is responsible for heat advection, where ${\overrightarrow{u}}_b$ is the blood velocity which depends on the blood flow rate Q according to ${\overrightarrow{u}}_b(\eta )=\frac{2Q}{\pi r^4}(r^2-{\eta }^2)$ (where η is the radial coordinate in the blood vessel cross-section) [15]. The flow rate is assumed to be dependent on blood vessel radius r according to Murray’s law. This is an heuristic expression describing the relationship between the pipe diameter and the flow rate within a flow network [16], given by: (2) $$Q={\left(\frac{1024\mu }{d_i^6{\pi }^2\xi [{\rho }_{\text{tube}}(c^2+2c)+{\rho }_{\text{fluid}}]}\right)}^{-1/2}.$$(2)

Here d_i, Q, μ, ξ, and c are the vessel diameter, volume flow-rate, dynamic viscosity of the fluid, power-to-mass ratio of the pumping system, and a constant material property of the pipes respectively. The constants ${\rho }_{\text{tube}}$ and ${\rho }_{\text{fluid}}$ are the pipe material density and fluid density respectively. The constant material property $\xi [{\rho }_{\text{tube}}(c^2+2c)+{\rho }_{\text{fluid}}]$ is independent of the vessel radius and therefore its value has been derived back from Equation (2)(2) $$Q={\left(\frac{1024\mu }{d_i^6{\pi }^2\xi [{\rho }_{\text{tube}}(c^2+2c)+{\rho }_{\text{fluid}}]}\right)}^{-1/2}.$$(2) using known values for the portal vein for the other parameters. In particular, the values used are: $Q=1.57e-5{\mathrm{\hspace{0.17em}m}}^3{\mathrm{\hspace{0.17em}s}}^{-1}$ and $\mu =0.00365$ Pa s for hepatic vessels with d_i = 10 mm [17], thus resulting in a value of $\xi [{\rho }_{\text{tube}}(c^2+2c)+{\rho }_{\text{fluid}}]=93.44{\text{\hspace{0.17em}kg s}}^{-3}{\mathrm{\hspace{0.17em}m}}^{-1}.$

The expression ${\sigma }_i|\nabla {\varphi }_i{|}^2$ represents the RF heat source densities for $i\in \{t,b,n\}.$ The electric potential is governed by the quasi-static approximation (3) $$\nabla ·{\sigma }_i\nabla {\varphi }_i=0$$(3) where $i\in \{t,b,n\}$ [1]. The power input to the domain was chosen to be dependent on the active length lambda and equal to $P_{\text{in}}=301.2\lambda$ W where λ is expressed in meters. The chosen value is similar to the one used for the model validation with an active length of 3 cm (see Section ‘Model validation’) and enough to lead to typical temperatures in RFA procedures (below 100$°$ C).

The top and curved surfaces of the tissue boundary were electrically insulated, while the bottom surface had ${\varphi }_t=0$ V. A time-dependent voltage (4) $${\varphi }_{\text{top}}(t)=\frac{P_{\text{in}}}{{\int }_{\text{Active Surface}}\overrightarrow{J}(t)·\widehat{n}dS}$$(4) was applied at the top boundary of the RF-needle, where $\overrightarrow{J}(t)$ is the current density. All external boundaries of the blood vessel domain were electrically insulated. All tissue-vessel interfaces were assumed to be electrically conducting. All tissue-needle interfaces except the active surface were electrically insulated.

The damage to the tissue and blood due to the RF heat source was modeled using the Arrhenius model [18,19]. The Arrhenius damage index Ω is related to the survival fraction ω according to (5) $$\omega =e^{-\Omega }.$$(5)

The tissue/blood survival fraction $\omega (\overrightarrow{x},t)$ evolves in time according to (6) $$\frac{\partial \omega (\overrightarrow{x},t)}{\partial t}=-A·\hspace{0.17em}\text{exp}\hspace{0.17em}(\frac{-\Delta E_a}{RT(\overrightarrow{x},t)})\omega (\overrightarrow{x},t).$$(6)

For both tissue and blood, the initial condition $\omega (\overrightarrow{x},0)=1$ (i.e., Ω = 0) was used. The values for the activation energy $\Delta E_a,$ ideal gas constant R, and frequency factor A for the tissue and blood Arrhenius models are given in Table 1. Tissue and blood were considered to be thermally damaged for survival fraction $\omega (\overrightarrow{x},t)<$ e^–1 (0.37) [19].

A simplified model was used to describe the effect of thermal coagulation of blood on blood flow. At any time step, the latest temperature field inside the blood vessel was used to calculate the survival fraction $\omega (\overrightarrow{x},t)$ of the blood. For this calculation, a number of blood ‘particles’ distributed over the vessel cross-section were considered, as shown in Figure 2. The local temperature $T_{\text{bl}}(z)$ seen by each blood particle as it traversed the length of the blood vessel was recorded. Here, z is the coordinate along the path of the blood particle, obtained using the relation (7) $$z=\left|{\overrightarrow{u}}_{\text{bl}}\right|t$$(7) (where t is the time and ${\overrightarrow{u}}_{\text{bl}}$ the velocity of the blood particle). This temperature-time data was used to evaluate the Arrhenius damage integral to solve Equation (6)(6) $$\frac{\partial \omega (\overrightarrow{x},t)}{\partial t}=-A·\hspace{0.17em}\text{exp}\hspace{0.17em}(\frac{-\Delta E_a}{RT(\overrightarrow{x},t)})\omega (\overrightarrow{x},t).$$(6) for each blood particle. If numerous such particles are used, an approximate representation of the coagulation state of the blood vessel cross-section is obtained. When this calculation showed all the blood particles in the cross-section of the vessel as coagulated, the blood flow was switched off and the simulation was resumed using the electrical conductivity of coagulated blood (given in Table 1). All other properties of coagulated blood were assumed to be the same as those for healthy blood. To the best of our knowledge, this is the first time a weak coupling between blood coagulation and RF ablation has been considered.

Figure 2. Paths followed by blood particles traversing the blood vessel. The definition of the z coordinate used in the calculation is also shown.

The transient simulation was run for a total physical time of 540 s. This is the typical RFA treatment time used in clinical practice. The tissue and blood were kept at an initial temperature of $37°$ C, i.e., the normal body temperature. The finite volume method based simulation package STAR-CCM + was used to solve the equations.

Model validation

An RF ablation experiment was performed in a cube (side length 10 cm) of thermochromic gel phantom using a RITA^® 1500X RF generator and the UniBlate^® RF-needle (both by Angiodynamics, New York, NY, USA) with active length set to 3 cm. Temperature set-point, ablation time, and maximum power of the RF generator were 103$°$ C, 10 min, and 35 W, respectively.

For the phantom, 671 g of distilled water was mixed with 175 g of a 40% (Bis-)Acrylamide (National Diagnostics, Atlanta, GA, USA) 19:1 solution. Fifty grams of irreversible change ink concentrate (Kromagen Magenta MB60-NH, TMC Hallcrest, Flintshire, UK) and 7 g of sodium chloride (Sigma Aldrich, St. Louis, MO, USA) were added to the solution and mixed under continuous stirring. The ink is known to discolor proportional to the maximum achieved temperature [20]. The catalyst N,N,N0,N0-tetramethylethylenediamine (1.4 g) was added before 1.4 g of ammonium persulphate (both from Sigma Aldrich, St. Louis, MO, USA) to initiate polymerization. The mixture was poured into a cube-shaped container and allowed to polymerize at room temperature.

The thermochromic phantom was prepared as a cube of side length 10 cm, as seen in Figure 3, which includes a 3D schematic of the experimental geometry. The RF needle was inserted perpendicularly at the center of the phantom cube face (as seen in Figure 3(a)), with insertion depth as shown in Figure 3(b). During the ablation, the phantom was placed on a table such that the inserted RF-needle was parallel to the table surface. The RF generator ground-pad was fixed to the surface of the phantom resting on the table.

Figure 3. (a) 3D schematic of the geometry showing the phantom cube and RF-needle orientation. (b) Schematic of the top view of the phantom geometry. Cube dimensions and RF-needle dimensions are shown.

In the simulation, the geometry of the experimental setup was replicated. Equations (1)(1) $$\begin{array}{ll}{\rho }_t{C}_t\frac{\partial T_t}{\partial t}\hfill & =\nabla ·k_t\nabla T_t+{\sigma }_t|\nabla {\varphi }_t{|}^2,\hfill \\ {\rho }_b{C}_b\left(\frac{\partial T_b}{\partial t}+{\overrightarrow{u}}_b·\nabla T_b\right)\hfill & =\nabla ·k_b\nabla T_b+{\sigma }_b|\nabla {\varphi }_b{|}^2,\hfill \\ {\rho }_n{C}_n\frac{\partial T_n}{\partial t}\hfill & =\nabla ·k_n\nabla T_n+{\sigma }_n|\nabla {\varphi }_n{|}^2,\hfill \end{array}$$(1) and (3)(3) $$\nabla ·{\sigma }_i\nabla {\varphi }_i=0$$(3) were used with the phantom material properties given in Table 1. There was no blood subdomain. The top surface of the RF-needle was kept at room temperature. Convection boundary conditions of the form (8) $$-k\nabla T·\widehat{n}=h(T-T_{\text{ref}}),$$(8) were applied to all other external boundaries. Symbols k, T, $\widehat{n},$ h, and $T_{\text{ref}}$ are the phantom thermal conductivity, temperature, unit surface vector, boundary convection coefficient, and reference temperature respectively. The values of the convection coefficient h used are shown in Table 2. For all these surfaces $T_{\text{ref}}$ was set as the room temperature. The room temperature (and initial temperature) was ${20}^{°}$ C. The bottom surface of phantom cube was electrically grounded. The voltage corresponding to the input power obtained from the RF generator logging data was applied to the top surface of the RF-needle.

Table 2. Convection coefficient values for external boundaries. [21].

Boundary	Convection Coefficient h (W m^–2K^–1)
Top Phantom Face (phantom-air interface)	12
Bottom Phantom Face (face resting on table)	9
Phantom Side Faces (parallel to ablation needle and flow-channel)	12
Phantom Side Faces (perpendicular to RF-needle)	10
RF-needle – air Interface	12

The average power supplied by the RF generator in the experiment was 10.22 W, which is similar to the power used in this parametric study. The thermochromic phantom turns magenta in regions where the temperature exceeds 70${}^{°}$ C. The comparison of the experimental and simulation results is shown in Figure 4. The magenta coloration of the thermochromic phantom marks the region where temperature exceeded 70${}^{°}$ C during the ablation. The blue line marks the boundary of this zone as predicted by the simulation. The experimental $T>{70}^{°}$ C zone area was 3.71 cm², while the simulated one was 3.17 cm². These results gave confidence on the method and on the simplification of the needle geometry.

Figure 4. Axial cross-section of post-ablation thermochromic gel phantom. During the experiment, the temperature in the phantom exceeded 70$°$ C in the magenta zone. The blue outline marks the same zone predicted by the model.

Simulations performed

Blood vessel configurations – SET A

A full central composite Design of Experiments (DoE) has been considered for the generation of SET A of the blood vessel configurations. Five different values for each parameter were considered which are given in Table 3. Ten (2 × 5) simulations were performed with one parameter taking either a ‘Max’ or ‘Min’ value and the other parameters taking ‘Center’ values. Thirty-two (2⁵) simulations were performed with each parameter taking either a ‘High’ or a ‘Low’ value. One simulation was performed with all parameters taking center values. This led to a total of 43 simulation runs. One reference simulation, which had no blood vessel region, was performed for every value of the RF-needle active length λ.

Table 3. Parameter values used for simulations of SET A.

Parameter	Min. Value	Low Value	Center Value	High Value	Max. Value
d1 (mm)	−7	−5	0	5	7
2 d2 (mm)	0.5	1	3	5	5.5
α (${}^{°}$)	0	30	45	60	90
r (mm)	0.5	1	3	5	5.5
λ (mm)	5	10	15	20	25

Blood vessel configurations – SET B

In these simulations, the blood vessel was parallel to the RF-needle. The distance d₁ between the RF-needle and blood vessel was kept 0. The perpendicular distance between the RF-needle and the vessel, d₂, was 1 mm. The RF-needle active length was λ = 20 mm. Seventeen different values were used for the vessel radius (Table 4).

Table 4. Vessel radii used for simulations of SET B.

No.	r (mm)	No.	r (mm)
1	0.25	10	0.70
2	0.30	11	0.75
3	0.35	12	1.00
4	0.40	13	1.25
5	0.45	14	1.50
6	0.50	15	1.75
7	0.55	16	2.00
8	0.60	17	5.00
9	0.65

Metrics and method used to quantify blood vessel effects

Possible effects of the blood vessel on the thermal lesion are:

• reduction in size of the lesion due to cooling (as shown in Figure 5)

  Figure 5. Schematic of the cooling effect (black outline) and the directional effect (green line) of a blood vessel on a thermal lesion. The cooling effect reduces the lesion size, while the directional effect causes stretching in the flow direction. Definition of ${\Delta }_{\text{up}},{\Delta }_{\text{down}},$ and $l_{\text{ref}}$ used to define the Δ_A directional effect metric shown. The cross-section is the one containing the vessel axis and the shortest line segment connecting the vessel axis and the RF-needle axis.

  

• stretching of the lesion in the direction of the blood flow (as shown in Figure 5, referred to as ‘the directional effect’ in this paper)

• combined effect resulting in both size reduction and lesion stretching.

To capture these effects and to quantify the impact of parameter variation on the thermal lesion, three metrics have been defined:

1. Average Lesion Boundary Displacement Δ_A (%): (9) $${\Delta }_A=\frac{{\Delta }_{\text{up}}+{\Delta }_{\text{down}}}{2l_{\text{ref}}}$$(9) where, ${\Delta }_{\text{up}},{\Delta }_{\text{down}},$ and $l_{\text{ref}}$ are defined as shown in Figure 5. ${\Delta }_{\text{up}}$ and ${\Delta }_{\text{down}}$ are positive in the blood flow direction and negative otherwise. The metric Δ_A represents the directional stretching of the thermal lesion as a percentage of reference lesion length. A large positive value of Δ_A corresponds to a large directional effect and vice versa.

2. Maximum Nusselt number ${\text{Nu}}_{\text{max}}:$ The maximum Nusselt number on the blood vessel surface at the end of the ablation was used as a proxy to study the effect of parameter variation on the heat transfer coefficient [21]. It was calculated using the expression on the extreme right: (10) $${\text{Nu}}_{\text{max}}=\underset{{}_{{\Gamma }_v}}{\max }\left(\frac{hL}{k}\right)=\underset{{}_{{\Gamma }_v}}{\max }\left(\frac{\dot{q}L}{k(T-T_{\text{ref}})}\right).$$(10) Here, Γ_v is the blood vessel surface, $h(\overrightarrow{x})$ is the heat transfer coefficient on the tissue-vessel interface, $\dot{q}(\overrightarrow{x})$ is the radial heat flux through Γ_v at the end of the ablation, L is the characteristic length (assumed to be the vessel diameter), k is the blood thermal conductivity (0.52 W m^–1 K^–1), T is the temperature on the vessel surface at the end of the ablation, and $T_{\text{ref}}$ is the reference temperature (310 K).

3. Relative ablated volume $V_{\text{rel}}:$ (11) $$V_{\text{rel}}=\frac{V}{V_{\text{ref}}}$$(11) This is the ratio of the ablated volume for a particular vessel and RF-needle configuration to the appropriate reference ablated volume.

Minitab^® Statistical Software¹ was used to fit second-order response surfaces to the three metrics (obtained from the simulations of SET A) as functions of d₁, d₂, r, α, and λ. Backward elimination was used to remove the insignificant terms in the response surface expressions. The terms with a p value lower than 0.1 (i.e., a confidence level of 90%) were removed during the backward elimination procedure.

Results

Quantitative data about the thermal lesions (i.e., lesion volume and maximum diameter) obtained in each simulation are given in the online supplemental material.

Results for SET A

The most important contribution to the directional effect metric Δ_A is of parameter d₁, as visible in the main effect plots obtained from the fitted response surface in Figure 6(a). A non linear relationship with the radius r was found. For the simulations of set A, no directional effect was observed on the thermal lesion (thermal lesion plots not shown), and the Δ_A values were within the range [−3%, 3%].

Figure 6. Results obtained from simulations of SET A.

Figure 6(b) shows the main effects plot for the maximum Nusselt number (Nu${}_{\text{max}}$) on the vein surface obtained from the fitted response surface. The important parameters affecting this metric were d₂, r, and λ. The most important effect was of variation in r. The following closed form expression, which approximates the relationship between Nu${}_{\text{max}}$ and the parameters d₂ (m), r (m), and λ (m) was obtained (R-sq = 96.75%): (12) $${\text{Nu}}_{\text{max}}(d_2,r,\lambda )=4.80-763d_2+9340r-444\lambda -1262379r^2+329993r\lambda$$(12)

(valid in the parameter ranges analysed).

The relative ablated volume $V_{\text{rel}}$ decreased with decreasing distance between the RF-needle and the blood vessel. The cooling effect was the largest when the blood vessel was parallel to the RF-needle and decreased with increasing angle till the vessel became perpendicular to the RF-needle. The relative ablated volume increased with increasing RF-needle active length.

Results for SET B

With decreasing blood vessel radius, two opposing effects become stronger:

1. Lesion stretching due to the directional effect

2. Blood coagulation

The objective of the simulations of SET B was to investigate which of these two effects dominates at which vessel radius and its implications for the ablated zone.

Figure 7 shows the estimated blood coagulation plots for all the simulations of SET B at the end of the ablation (Figure 7(b–h)) and the blood vessel cross-section at which thermal damage to blood is computed (Figure 7(a); red line). The complete vessel cross-section coagulated for 0.25 mm $\le r\le$ 0.40 mm (Figure 7(b)), while the blood remained unaffected for 0.7 mm $\le r\le$ 5 mm (Figure 7(h)). For 0.45 mm $\le r\le$ 0.65 mm (Figure 7(c–g)) the blood cross-section suffered partial thermal damage which was not sufficient to cause flow blockage. For 0.25 mm $\le r\le$ 0.40 mm the time at which the blood cross-section completely coagulated is given in Table 5.

Figure 7. (a) Location of the vessel cross-section at which thermal damage is calculated. (b–h) Arrhenius survival fraction ω for the vessel cross-sections for SET B computed at the blood outlet. A value of $\omega <0.37$ corresponds to blood coagulation. The RF-needle was located below the cross-section.

Table 5. Time required for complete coagulation of the blood vessel cross-section.

Radius (mm)	Time for Complete Blood Coagulation (s)
0.25	65
0.30	76.5
0.35	126
0.40	238

Figures 8 and 9 show the thermal lesions obtained at the end of the simulation. For r = 0.25 mm, r = 0.3 mm, and r = 0.35 mm, the thermal lesion was ellipsoidal as a result of early flow-blockage due to blood coagulation. For r = 0.4 mm, the blood coagulation occurred late resulting in a tail-like structure, which was not fully enveloped by thermal lesion growth after flow blockage. For r = 0.45 mm, the thermal lesion had a long tail-like structure. For r = 0.5 mm a small tail was visible in the tissue thermal lesion, whereas for r > 0.5 mm no tail like structure was present. Based on these findings, a classification depending on the likelihood of the directional effect can be derived as shown in Figure 10. For vessels of ‘region 1′ there is no directional effect of blood flow on the thermal lesion, so an advection term is not required in the blood subdomain. A perfusion term with a higher perfusion coefficient (e.g., the model from Kröger et al. [6]) would suffice. For vessels of ‘region 2′ directional effects do occur, and an advection-diffusion equation would be necessary to model heat transfer. The blood vessels of ‘region 3′ have neither a directional effect nor a cooling effect, and can be modeled similar to solid tissue, but with different material properties.

Figure 8. Thermal lesion at the end of the simulation for 0.25 mm $\le r\le$ 0.7 mm. The tissue section containing the vessel axis and the RF-needle axis is shown. The top horizontal band in each image is the blood vessel section. Blood flow was from right to left. The bottom horizontal ‘half-band’ is the RF-needle. Blue indicates healthy tissue and red indicates ablated tissue. The white circle is the tumor boundary.

Figure 9. Tissue thermal damage at the end of the simulations for 0.75 mm $\le r\le$ 5 mm. Sub-figure (h) shows the reference lesion. The tissue section containing the vessel axis and the RF-needle axis is shown. The top horizontal band in each image is the blood vessel section. Blood flow was from right to left. The bottom horizontal ‘half-band’ is the RF-needle. Blue indicates healthy tissue and red indicates ablated tissue. The white circle is the tumor boundary.

Figure 10. Classification of blood vessels into regions based on the occurrence of a directional effect.

Figure 11 shows the effect of radius variation on the three scalar metrics in SET B. The directional effect metric Δ_A first increased up to r = 0.45 mm and then decreased as r increased (Figure 11(a)). For the cases where a tail-like structure was visible in the thermal lesion the value of Δ_A was larger than 5%. The maximum Nusselt number increased with increasing radius (Figure 11(b)). For the radii 0.3 mm, 0.35 mm, and 0.4 mm the lesion volume was larger than the reference lesion (Figure 11(c)). This could be a combined effect of the high electrical conductivity of healthy blood, the low blood flow-rate, and the blockage of blood flow during the simulation. From vessel radius 0.45 mm onward $V_{\text{rel}}$ decreased with r.

Figure 11. Effect of vessel radius on (a) directional effect Δ_A, (b) Nusselt number ${\text{Nu}}_{\text{max}},$ and (c) relative ablated volume $V_{\text{rel}}$ metrics for simulations of SET B.

Discussion

Our work extends the results of dos Santos et al. [8] through an expression for the heat transfer coefficient as a function of d₂, r, and λ. In SET A for the variation of cooling effect with parameters r and α, our findings agree with those of Chen and Saidel [3] and Huang [7]. For SET B, where blood coagulation came into play, our results showed that the vessel cooling effect first decreased with increasing r and then increased. We observed that, in cases where blood coagulation occurred during the ablation, the ablated volume was larger than the reference. Our study found that, as in the one by Chen and Saidel, the vessel cooling effect decreases continuously as d₂ increases. The study by Huang found that the cooling effect of a blood vessel first increases as d₂ is increased and reaches a maximum value, after which it decreases again. To the best of our knowledge the present study is the first one to investigate how varying the parameter d₁ influences the blood-vessel effects.

Huang found that, ignoring blood coagulation the directional effect occurs for small blood vessels (irrespective of flow rate) or large vessels with low flow rate. We observed that directional effects occurred for vessel radii between 0.4 mm and 0.5 mm for flow rates dictated by Murray’s law. To the best of our knowledge, the present study is the first one that has taken into account blood coagulation in a systematic investigation of the directional effect of blood vessels in RFA.

The present work uses temperature dependent expressions for the thermal and electrical conductivities of tissue and blood. This is in contrast with related works in literature, see e.g., [3,7], which use constant values. The strong dependence of Nu${}_{\text{max}}$ on the vessel radius suggests that a constant Nusselt number would either over- or under-estimate the cooling effect depending on the vessel radius. Hence, the closed-form expression for Nu${}_{\text{max}}$ could potentially improve the accuracy of the simulation results. The findings concerning directional effects highlight the need for using an appropriate, advection-based, mathematical model for 0.4 mm $\le r\le 0.5$ mm. Furthermore, the findings suggest that the bio-heat model used by Kröger et al. [6] is valid for 0.5 mm $\le r\le$ 5 mm.

The results we obtained depend on the RF-needle shape and are specific to the Angiodynamics UniBlate^® RF-needle. The effect of a single blood vessel on the thermal lesion is considered. The number of blood vessels surrounding the RF-needle is known to have an impact on the thermal lesion during RFA. Another limitation of the model used is the weak coupling of the blood coagulation, blood flow, and energy models. The classification of blood vessels does not include the effect of parameters other than the vessel radius. Future work can focus on: (i) addressing limitations discussed above and (ii) experimental validation of the conclusions obtained.

Conclusions

A mathematical model has been used to determine the impact of blood vessel parameters and RF-needle active length on: (i) the ablated volume, (ii) the directional effect, and (iii) the heat transfer coefficient at the vessel wall.

Firstly, the vessel radius r, distance d₂, and the RF-needle active length λ are the most important parameters affecting the maximum Nusselt number on the blood vessel surface. An expression for the maximum Nusselt number as a function of these parameters has been obtained, which can be used to determine the heat transfer coefficient in approximate models.

Secondly, directional effects of blood are important for 0.4 mm $\le r\le 0.5$ mm. As a first approximation, mathematical models should allow for directional effects for blood vessels of radii between 0.4 mm and 0.5 mm.

Thirdly, a vessel classification scheme based on the occurrence of a directional effect has been proposed, which can be used to develop a mathematical model tailored for fast estimation of vessel cooling effects.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Additional information

Funding

This work was supported by the European Commission under the Marie Skłodowska Curie Actions – Innovative Training Networks (MSCA-ITN) grant of the type European Industrial Doctorate (EID) [project number 642445].

Notes

1 MINITAB^® and all other trademarks and logos for the Company’s products and services are the exclusive property of Minitab Inc. All other marks referenced remain the property of their respective owners. See minitab.com for more information.