Laser interstitial thermal therapy of lung lesions near large vessels: a numerical study

Objective. Laser interstitial thermal therapy (LITT) is an evolving hyperthermia-based technology that may offer a minimally invasive alternative to inoperable lung cancer. LITT of perivascular targets is challenged by higher risk of disease recurrence due to vascular heat sinks, as well as risk of damage to these vascular structures. The objective of this work is to examine the impact of multiple vessel parameters on the efficacy of the treatment and the integrity of the vessel wall in perivascular LITT. Approach. A finite element model is used to examine the role of vessel proximity, flow rate, and wall thickness on the outcome of the treatment. Main result. The simulated work indicates that vessel proximity is the major factor in driving the magnitude of the heat sink effect. Vessels situated near the target volume may act as a protective measure for reducing healthy tissue damage. Vessels with thicker walls are more at risk of damage during treatment. Interventions to reduce the flow rate may reduce the vessel’s heat sink effect but may also result in increased risk of vascular wall damage. Lastly, even at reduced blood flow rates, the volume of blood reaching the threshold of irreversible damage (>43 °C) is negligible compared to the volume of blood flow throughout the treatment duration. Significance. This investigative simulation yields results that may help guide clinicians on treatment planning near large vessels.


Introduction
Lung cancer remains a significant healthcare challenge in Canada and across the globe [1]. Factors such as education around prevention, improved treatments, and implementation of screening protocols have contributed to a decline in age standardized mortality since peaking in the 1990s [1]. Screening has increased early-stage detection, where treatment is primarily directed towards local control. The preferred standard of care for early-stage disease is surgical resection. Unfortunately, not all patients are eligible for this treatment due to comorbidities. Development and deployment of minimally invasive tools may offer alternative treatment options for such patients.
Stereotactic body radiation (SBRT) is the recommended treatment for medically inoperable patients. SBRT is a highly effective treatment but remains doselimited due to radiation toxicity. Advancements in minimally invasive surgical tools, such as bronchoscopic devices, have introduced the opportunity for direct local treatment using different energy devices, which are not limited by the radiation risks of SBRT.
Laser interstitial thermal therapy (LITT) is a hyperthermia-based treatment that deposits light energy directly into the tissue, which is absorbed and released as vibrational energy to increase local tissue temperature. Tissue temperature increases to ∼55°C result in irreparable damage, such as irreversible protein denaturation, leading to cell death. A concern with LITT, and other hyperthermia-based treatments, is the complexity in treatment planning when vessels are present near the target. Vessels may impact the treatment in two ways. First, the vessels may act as a heat sink, resulting in incomplete treatment and risk for disease recurrence. This concern has been reported by several institutions in retrospective analyses of hepatic radiofrequency ablation (RFA) [2,3], and demonstrated experimentally with RFA of hepatic lesions in pigs [4]. Their studies showed a transition zone at a vessel diameter between 2-4 mm, where smaller vessels are completely occluded, and larger vessels remain patent. These larger vessels act as a heat sink that reduces the size of the final ablation zone. The second concern is potential unintended vascular wall damage during treatment. Lu et al also demonstrated that while larger vessels do remain patent, signs of vascular damage are observed, such as endothelial necrosis and luminal thrombus formation [4]. Under similar treatment conditions, vessels of different types and similar size may exhibit different responses to local hyperthermia, including different expressions of thrombosis-related genes and proteins [5]. Treatments using different types of energy devices are variably impacted by the heat sink effect and in their risk of vascular damage [6,7], indicating that more work is needed to fully characterize the risk in using hyperthermia devices for treatments near large vessels.
To better understand the impact of large vessels on LITT, several investigators have turned to simulations. Modeling vessels embedded in the target volume during LITT showed that vessels placed closer to the source have a larger impact on the asymmetry of temperature distribution [8]. Further, doubling the vessel flow rate had no significant impact on temperature distribution, indicating that the rate of heat removal by the vessel had already reached a threshold at normal physiologic flows. Verhey et al included dynamic changes in tissue properties to their model and evaluated the impact of multiple variables, such as the laser input power, application time, flow rate, and distance to the applicator, on the final coagulation volume [9]. They showed that drastically reducing the flow velocity in the vessel may result in coagulation of blood inside the lumen. This is in line with the notion that vessels under a certain size, which have a smaller capacity to carry away heat, end up being completely occluded during treatment. Several investigators built on this work by developing their own models that consider more complex flow patterns and cases with multiple vessels [10][11][12].
Prior work has shown that the specific tissue properties and anatomy geometry play a significant role in the impact of local vasculature on LITT treatment. The work described herein aims to characterize the influence of vessel distance from target volume, lumen diameter, and wall thickness on the magnitude of the heat sink effect in the context of a treatment planning scenario. Hence, we consider the extra treatment time needed to achieve a successful treatment of the target volume while also assessing the impact on the surrounding normal tissue. We also evaluate the impact of the treatment on the viability of the vessel wall and the blood flowing through the vessel. To our knowledge, this is a novel addition to prior simulation work and has not been thoroughly examined before.

Methods
To evaluate the impact of different variables on LITT treatment near large vessels we developed a finite element model using the COMSOL software [COM-SOL Multiphysics ® v. 5.6. www.comsol.com. COM-SOL AB, Stockholm, Sweden]. The model was similar in structure to the one developed by Verhey et al [9]. The different components to the modelling process are discussed below.

Geometry
The geometry under investigation was comprised of 6 domains: (i) optical fiber (ii) tumor tissue (iii) margin volume (iv) blood in the vessel lumen (v) vessel wall and (vi) normal lung tissue (figure 1).
To represent the emitting surface of a cylindrical diffusing fiber typically used in LITT, the optical fiber was modeled as a cylinder 10 mm in length and 1 mm in diameter. The tumor was modeled as a sphere of 10 mm diameter, representative of early-stage disease that would be a candidate for minimally invasive techniques such as LITT [13]. A margin volume extending 5 mm beyond the edges of the tumor was included to account for the presence of microscopic disease that is not visible on radiological imaging. In radiation oncology, standard practice is to include a 5 mm margin [14]. Importantly, for our model the absolute volume in the margin varied depending on the vessel size and distance (figure 2). The blood vessel was modeled as a cylinder spanning the length of the geometry. The vessel wall thickness depends on the vessel size, type, state, age, and presence of disease. In general, typical wall thicknesses range between 5%-10% of the lumen diameter [15][16][17][18]. We compared both ends of this range in our model. These five domains were then encased in a rectangular volume to represent the surrounding healthy tissue. The volume was 80 mm × 80 mm × 150 mm; the extended length was in the direction of fluid flow to allow for sufficient opportunity for heat dissipation into the surrounding tissue. The size was guided by prior work done on LITT modelling with experimental validation. It was heuristically optimized to minimize computation time and comply with the boundary conditions (BCs) as well as continuity assumptions of the governing equations [19,20].  For example, at the same distance, a case with a smaller vessel will have a larger absolute margin volume due to diminished overlap between the margin and vessel domains; therefore, a larger volume of tissue must reach the damage threshold to achieve the treatment endpoint.

Governing equations & boundary conditions
To simulate LITT we used a dynamic model which coupled the light transport, heat transfer and damage accumulation in the tissue [9,19]. The governing equations and accompanying BCs are described below.

Light distribution
Photon distribution in the domain was calculated using the time-independent light diffusion equation, appropriate for use in turbid media, such as biological tissue [21]: where D is the diffusion coefficient, a m the absorption coefficient, s m¢ the reduced scattering coefficient, Q S is the light source, and f the fluence rate. In our case, the light source is modeled as a spatial BC representing the photon emission from the fiber surface: where P is the laser power and SA is the emission surface area of the fiber. A laser power of 3 watts was used. This was heuristically determined as the maximum permissible laser power that maintained the maximum temperature in the simulation domain below water's phase transition at 100°C. A sample temperature profile for one of the simulations is available in the supplementary documents (figure S1) The second BC accounted for diffuse reflection at the surface of the boundary limits [21], where it was assumed that the relative refractive index is 1 so that any light traversing the surface travelled through it and was not reflected. This assumption was tested by either increasing the volume or setting the relative refractive index to be greater than 1. Neither change affected the results.

Heat transfer
Heat transfer was modeled using a modified version of Pennes' bioheat equation, including an energy transport term that accounts for energy transfer due to flow in the large vessel [22]: where ρ is the density, c is the specific heat capacity, T is the temperature, u is the velocity field in the vessel, k is the thermal conductivity, w is the blood perfusion rate for each tissue type, the subscript B denotes the use of general blood properties, T ¥ is the average tissue temperature at rest, Q M is the heat generated through metabolism, a m is the absorption coefficient, and f the fluence rate. The product of the absorption coefficient and the fluence rate represents the heat generated through light absorption.
An initial condition and two BCs were applied to the thermal model. The initial condition set the starting temperature to match the average body temperature at rest, 37°C. The first BC defines the temperature at the surface of the treatment fibre; we modeled the clinical scenario where the treatment fiber is normally encased in a jacket providing a steady stream of cooling liquid. Like previous investigators, we modeled this surface as a constant temperature of 37°C [20]. This was sufficient to ensure that the temperature in our domain remains below water's phase transition. The second BC set the edges of the geometry at a constant temperature of 37°C to represent that the surrounding tissue is not impacted by the treatment.
Flow velocity was defined using Hagen Poiseuille's law for pressure driven flow of a Newtonian fluid through a cylindrical tube [23]. It was assumed that the flow is fully developed, and therefore there was no net momentum flow. Further, it was assumed that velocity at the walls was zero allowing us to model flow velocity profile in the vessel as: where r is the radial distance from the center of the vessel, u(r) is the velocity as a function of the vessel radius, u avg is the average velocity, and R is the radius of the vessel.

Damage accumulation
Damage accumulation in the tissue was calculated using the Arrhenius equation, which estimates the rate of protein denaturation as a proxy for cell death [24,25]: where ( ) t W is the damage index as a function of time, A is the frequency factor, E D is the activation energy, R is the universal gas constant, T is the temperature of the tissue, t is the total heating time, and c L is the concentration of living tissue as a function of time.

Tissue properties
To improve the clinical relevance of our work, an extensive literature review was conducted to identify the most appropriate tissue properties for the simulation. The complete set of properties are listed in table 1. The optical properties chosen were all at, or near, 1064 nm, a laser wavelength commonly used for LITT due to its deep optical penetration. It was assumed that absorption remains unchanged after coagulation. To estimate the properties of inflated lung, we used the work by Beek et al on the lung's optical properties as a function of respiration to modify in-vivo deflated lung properties collected by Spliethoff et al [26,27]. We assumed that after coagulation the native reduced scattering value is doubled for the lung.
The thermal properties were collected from the IT'IS database for thermal and electromagnetic parameters of biological tissues [35] and assumed to remain constant throughout treatment. The tumor was assumed to have the same properties as the normal tissue. This simplification was set up to reduce the number of dynamic factors affecting the results and maintain the focus on the impact of the vessel's cooling on the LITT response.
To estimate the tissue damage, we chose to use the Arrhenius model for thermal denaturation as opposed to an isotherm contour to allow us to take into consideration the temperature and time dependance of tissue exposure as well as the unique response of the different tissue types. The kinetic parameters used in the Arrhenius model are tissue specific and commonly determined experimentally. They are paired through a linear relationship between the activation energy and the logarithm of the frequency factor [25]. Random errors in experimental measurements can lead to mismatches in the pairings and result in large errors in the denaturation predictions. For this reason, we used the activation energy values from prior literature in conjunction with the empirical correlation established by Qin et al to determine the matching frequency factor pair [25]. To the best of our knowledge, the damage properties for the lung are not available in the literature, therefore we assigned the lung similar values as collagen as it is an abundant protein that is critical for stability and structure of the lung [37,38].
The adaptive mesh refinement tool in COMSOL was used to optimize the model run time for accurate conversion. This tool generates multiple meshes for segments of a time-dependent simulation based on areas with high error estimates. For example, the case with a thick-walled 16 mm vessel diameter, situated 5 mm away from the tumor edge had an initial mesh comprised of 45 885 elements and a final mesh of 187 051 elements after the refinement process. The quadratic Lagrange shape function was used for the light and heat transfer. A quadratic discontinuous Lagrange shape function was used for the thermal damage calculation. The default backwards differentiation formula was used to automatically adjust the time step throughout the solver with a default tolerance set as 0.01. The simulations were performed on an Intel Core i7-8700 ASUS PC with 64 GB of RAM.

Study parameters
Target volume and normal tissue treatment outcome parameters were used to assess the impact of the vessel cooling. Using our previously defined geometry, the target volume (TV) is the tumor volume plus the margin volume (5 mm beyond the tumor perimeter). The treatment end point was set to be the point when 90% of the TV had reached the defined coagulation threshold of an Arrhenius index value of 1; this is defined as t90. The microscopic disease is not expected to be present throughout the full margin volume and hence complete margin coverage is not required. In radiation oncology, acceptable treatment coverage varies between 90% to 98% coverage. While we have chosen 90% for this study, other coverage values may be used depending on clinical experience and need.
The absolute volume of normal tissue (excluding the vessel wall) that reached the defined coagulation threshold of an Arrhenius index value of 1 is noted at t90 and is referred to as NT_Vc (Normal Tissue Volume coagulated). An example treatment calculation is shown in figure 3. The times to reach TV coverage of 90%, t90, are shown as vertical dotted lines that intersect the LITT response curve at the 90% point ( figure 3(a)). These same t90 vertical lines are shown on the normal tissue LITT response curves ( figure 3(b)). The NT_Vc for each treatment is shown as the horizontal dotted lines starting from the intersection of the t90 and response curves back to the vertical axis. An optimal treatment would achieve 90% TV coverage while minimizing NT_Vc. Using these treatment outcomes as metrics, we compared the impact of changes in model as described here.

Vessel distance and diameter
The distance between the edge of the tumor and the outer edge of the vessel wall (as highlighted in figure 1(b)) was varied between 1, 3, and 5 mm, representing vessels sitting near the edge of the tumor and those at the outer edge of the margin volume. These were chosen as they represent vessel distances that may impact the treatment outcome of the clinically relevant target. The vessel lumen diameter was varied between 4, 10, and 16 mm. This represented a wide spectrum of vessels that encompass the upper end of the 2-4 mm transition zone identified by Lu et al and was smaller than vessels expected to exhibit turbulent flow, which was not implemented in our current model [4]. The 16 mm diameter vessel represented an average vessel in the secondary branching of the pulmonary artery, with a Reynolds Number under 2100, which is in the laminar range [23].

Blood flow rate
The flow rate of the vessels was calculated using the square law relationship between blood flow and vessel diameter [39]. The relationship is expected to be valid for the vessel sizes that are used in this work. For example, using a flow rate of 5 l min −1 for the aorta and a diameter of 25 mm, the relationship predicts a flow rate of 2.048 l min −1 for a 16 mm diameter secondary generation pulmonary artery, and 0.128 l min −1 for a 4 mm diameter coronary artery.
These values agree with published literature [39,40]. A reduction in flow rate to 25% of the physiological norm was also modeled to represent potential interventions to reduce the vessel's heat sink effect, such as vessel occlusion or flow reduction through pharmacological means.

Results
Vessel distance from tumor and diameter Figure 4 shows the changes in TV versus treatment time for vessel distances from the edge of the tumor of 1, 3 and 5 mm and vessel diameters of 4, 10 and 16 mm. The case where no vessel is present is shown as a benchmark comparison. The top row shows the relative volume of TV treated versus treatment time, while the lower row shows the absolute volume of normal tissue reaching coagulative thermal dose. The vessel flow rate was maintained at the physiological norm and the wall thickness was set to be 5% of the lumen diameter. Treatment effect was calculated to 2500 s of laser heating.    Figure 5 summarizes the results of t90 and NT_Vc for each vessel diameter versus vessel distance. The distance of the vessel from the tumor has a much greater effect on t90 than vessel diameter. The smallest vessel (4 mm diameter) results in the most significant heat sink effect, from the perspective of longest treatment time. This result may be explained based on our geometry and defined end point. As previously indicated in figure 2, the absolute margin volume is largest for the smallest vessel diameter. The excess tissue in that region is also more difficult to treat because (i) it is shielded from laser light due to the presence of the vessel and (ii) the vessel wall is thinner and so the conduction path for heat to be removed through the blood flow's convective effects is shorter. These effects are reduced as the distance between the vessel and the tumor increases.
The impacts of vessel diameter and the vesseltumor distance on the volume of normal tissue that is coagulated, NT_Vc, is more complex. When the vessel is close to the tumor, the higher t90 for the small vessel diameter case results in more normal tissue being coagulated. As the vessel-tumor separation increases to the perimeter of the margin, NT_Vc is lower than found in the scenario with no vessel. While treatment time is longer with the vessel present, the heat sink effect of the vessel is cooling normal tissue immediately outside the TV, protecting it from coagulation. In general, these results indicate that vessel vicinity is a more critical variable to consider when evaluating the risk of the heat sink effect, as opposed to the size of the vessel itself.
Vessel wall thickness and flow rate Figure 6 shows the effect of vessel wall thickness and blood flow rate on treatment outcome for the case with the most significant cooling effect: tumor-vessel distance of 1 mm. Two alternative scenarios are modeled: an increase in wall thickness from 5% to 10% of the diameter, and a reduction in blood flow to 25% the standard physiological rate. The TV LITT response curves change with increased wall thickness and decreased blood flow, leading to decreases in t90 for all vessel diameters between ∼10%-20%. In these simulations, the distance between the tumor and outer vessel wall remains the same. Hence changes in treatment response due to a thicker wall partly reflect the increase in distance between the heat absorbing fluid and the edge of the TV. The normal tissue LITT response curve is essentially the same for all treatment conditions, but the variations in t90 for the 3 scenarios leads to variations in NT_Vc, which is lower for both increased vessel wall thickness and reduced blood flow by virtue of the accelerated time-to-treatment.
Potential damage to the vessel wall during treatment Vessel wall damage was only observed for the case of 16 mm diameter vessel, with increased wall thickness (10%, i.e., 1.6 mm), and placed 1 mm from the tumor periphery. Results are shown in figure 7 (results for other cases can be found in figure S2). The maximum penetrative depth of vessel wall coagulation volume was 0.35 mm ( figure 7(b)) and confined to the outer wall of the vessel. Decreasing the blood flow rate (figure 7(c)) results in an increased maximum depth of the coagulation volume to 0.7 mm, slightly less than half of the total thickness of the wall. The penetrative extent of the damage seen in this case indicates that extreme care should be taken when treating near thick-walled vessels with artificially reduced flow rates.
The impact of the treatment on the temperature of blood in the vessel Our simulations show that LITT had negligible impact on the blood temperature ( figure S3). Even at reduced blood flow rates, the volume of blood reaching levels of irreversible damage (>43°C) is negligible compared to the volume of blood flowing through the vessel during the duration of treatment [41]. Additional care should still be taken to observe for potential thrombotic complications during and post LITT treatment near large vessels, since a minor volume of elevated temperature may yet act to activate the coagulation cascade.

Discussion & next steps
We assessed the impact of a blood vessel on LITT outcomes, assuming the vessel acts as a heat sink. Prior investigations, such as the work by Verhey et al and Paul et al have shown that the presence of large vessels may result in asymmetric temperature and thermal damage results in response to LITT [9,11,12,20]. In building on this work, we examined the impact of the vessel on a clinically relevant treatment planning scenario. Our work showed that the presence of the vessel affects treatment delivery times needed for adequate coverage of the target volume. Hence, rather than examining the changes in temperature or thermal dose profiles, as commonly done in prior work, the target volume response and the volume of unintended damage to normal tissue were used to assess the impact of the vessel. We assumed that treating 90% of the target volume would be sufficient for local control in most patients. This value may need to be adjusted based on the desired treatment margin to encompass the distribution of microscopic disease within the target volume. A treatment outcome metric of t90 was defined as the LITT time required to achieve treatment of 90% of the target volume. The use of dose-volume histograms (DSV) is standard practice in radiation therapy and is occasionally used in photodynamic and thermal therapies [42]. We chose to use treatment time over DSV as it made it easier to understand the effect that changing treatment conditions had on NT_Vc. Using treatment time was also valid because other treatment parameters, such as laser power and tissue properties, were kept constant. If these properties change, the results can be reformulated with respect to thermal dose.
The separation between the vessel and tumor dominates the heat sink effect, with longer treatment times needed for complete treatment when the vessel is closer to the tumor. The smallest vessel results in the longest treatment time, due to the specific geometry of the vessel location and the defined end point. The increased treatment times also increase the volumes of normal tissue that are coagulated. As the vessel-tumor separation increases, treatment delivery times converge and cooling by the vessel leads to reduced normal tissue coagulation. Although this is a beneficial effect in our defined scenario, it showcases the dangers of treating near large vessels since their cooling effect in proximity to the treatment target risks disease recurrence through undertreatment of nearby cells. Blood flow rate and vessel thickness had less of an effect on treatment outcomes compared to vessel position. In a clinical scenario, this may be overcome by adjusting the location of the treatment fiber or increasing the input power of the laser.
Potential damage to the vessel wall was also estimated in our models. We found the risk of damage to the vessel wall was present only in the largest vessel, is minor under normal physiological flow rates, and decreases with an increasing distance away from the tumor edge. Damage is highest at the outer edge of the vessel. It is unclear whether the magnitudes in these results have clinical significance but are notable since damage to the vascular wall could lead to catastrophic complications, such as aneurysm and/or hemorrhage from structural failure of the wall. This damage may also initiate the coagulation cascade and result in perioperative thrombotic complications. These concerns are exacerbated further with reduced vessel flow rate, as would occur with use of an atraumatic hemostat, catheter balloon, or pharmacological agents that modify flow directed through the vessel. The natural differences in wall thickness between different vessels in the body should be considered when evaluating whether it is safe to treat near vascular structures. The tumor-vessel distance was maintained at 1 mm, the case with the most significant cooling effects. The case where no vessel is present is shown as a benchmark comparison.

Limitations of the study
The current model has some limitations that need to be addressed. The effects of phase changes on the tissue properties and energy transfer are not considered, preventing evaluation of further increased laser power that would result in water evaporation due to heating above 100°C. The model also does not take into consideration the dynamic effects of the tissue's thermal properties with respect to temperature or damage state. This is largely attributed to changes in the tissue's water content, which is not currently considered in the model [43]. Lastly, the model is limited to using laminar flow profiles of blood transport through the tissue. Integrating more complex flow simulation techniques could enable examination of the role of shear stresses, turbulence, and bends in the vascular structure on potential vessel occlusion through overheating. The model is planned to form the basis for a treatment planning platform that includes modifying other treatment delivery parameters, including the use of multiple treatment fibers, as well as variations in light delivery over time or between treatment fibers. For example, in the case with the largest heat sink effect (vessel diameter 4 mm at 1 mm from the tumor edge), varying the number of treatment fibers and fiber location may provide more optimal coverage that minimizes normal tissue damage.
Further validation of the model's results should be carried out. We have confirmed the accuracy of the optical and thermal components of the model through comparison with prior experimental work on coagulating phantoms and ex-vivo tissue [19]. Experiments validating the fluid flow dependent results of the model should be conducted. Tissue mimicking thermochromic or coagulating phantoms may be used to evaluate the extent of thermal damage resulting from the prescribed treatment dose [44,45]. Tubing or exvivo vascular tissue may be used to act as a conduit for fluid flow. Once simulation accuracy has been established, experiments in pre-clinical animal models with large vessels may be carried out for assessment of target coverage, as well as histopathological examination of the extent of the thermal damage of the healthy tissue, including the vessel wall.

Conclusion
The simulated work indicates that large vessels, inside or near a target volume, create a heat sink effect that may increase the thermal dose required to achieve the desired treatment end point. Vessel proximity appears to be the major factor when it comes to meaningful heat sink effects in a treatment. Vessels near the treatment volume may act as a protective measure for reducing healthy tissue damage. Surgical or pharmacological interventions that reduce blood flow may reduce the vessel's heat sink effect but may also increase the risk of vascular wall damage. Vessels with thicker walls appear to be more at risk for damage due to a longer diffusion distance required before excess heat can be expelled through the vessel's outflow. A minor temperature rise in the blood appears to occur throughout treatment. This investigative simulation yields results that may help guide clinicians on treatment planning near large vessels.