Feasibility and relevance of discrete vasculature modeling in routine hyperthermia treatment planning

Abstract

Purpose: To investigate the effect of patient specific vessel cooling on head and neck hyperthermia treatment planning (HTP).

Methods and materials: Twelve patients undergoing radiotherapy were scanned using computed tomography (CT), magnetic resonance imaging (MRI) and contrast enhanced MR angiography (CEMRA). 3D patient models were constructed using the CT and MRI data. The arterial vessel tree was constructed from the MRA images using the ‘graph-cut’ method, combining information from Frangi vesselness filtering and region growing, and the results were validated against manually placed markers in/outside the vessels. Patient specific HTP was performed and the change in thermal distribution prediction caused by arterial cooling was evaluated by adding discrete vasculature (DIVA) modeling to the Pennes bioheat equation (PBHE).

Results: Inclusion of arterial cooling showed a relevant impact, i.e., DIVA modeling predicts a decreased treatment quality by on average 0.19 °C (T90), 0.32 °C (T50) and 0.35 °C (T20) that is robust against variations in the inflow blood rate (|ΔT| < 0.01 °C). In three cases, where the major vessels transverse target volume, notable drops (|ΔT| > 0.5 °C) were observed.

Conclusion: Addition of patient-specific DIVA into the thermal modeling can significantly change predicted treatment quality. In cases where clinically detectable vessels pass the heated region, we advise to perform DIVA modeling.

Keywords:

• Hyperthermia
• discrete vasculature
• hyperthermia treatment planning
• head and neck
• thermal modelling

Introduction

During hyperthermia cancer treatments, the temperature in the target region is increased to therapeutic levels of 40–44 °C to sensitize tumors cells for the effects of radiotherapy and/or chemotherapy treatments [1]. For head and neck (H&N) cancers, the effectiveness of hyperthermia has been shown in phase III clinical trials [2–5].

Hyperthermia dose–effect relationships show that an increase of temperature will further improve treatment outcome [6–8]. Our H&N hyperthermia applicator (HYPERcollar3D [9]) allows conformal focused microwave heating of tumors located deeply in the entire H&N region. Mandatory in our treatment strategy is patient-specific hyperthermia treatment planning (HTP) to maximize the microwave power absorption inside the target region while minimizing exposure of sensitive healthy tissues. To improve this procedure, we are investigating replacing the power absorption-based optimization by optimization using patient-specific thermal simulations. In this work, we studied if incorporation of vasculature segmentation and discrete vasculature (DIVA) modeling in routine planning would be feasible and if this would result in a relevant change in temperature prediction.

Accurate HTP requires patient specific information: positioning of the patient in the applicator, delineation of the target volume (CTV: clinical target volume as used in radiotherapy planning), tissue segmentation and electrical and thermal properties of the tissues [10,11]. For the H&N region, we have developed an auto segmentation routine [12,13] and showed that it performs within intra-observer variations [14]. Based on this method, we showed that patient specific thermal properties, i.e., blood perfusion and thermal conductivity, can strongly improve 3D temperature simulation accuracy [15]. That study showed the crucial importance of taking patient specific cooling due to perfusion into consideration. Still, the impact of cooling by the large vessels on top of this microvasculature cooling is unknown.

Cooling in tissues is generally modeled using the bioheat equation described by Pennes [16]. Pennes modeled the heat removal by blood using a homogeneous heat sink term, which scales proportional to the tissue temperature increase above the blood temperature. Although the Pennes model (PBHE) takes the cooling by capillaries adequately into account, the effect of DIVA is ignored which may lead to inaccurate temperature predictions [17,18]. Several thermal models to describe the heat exchange between vessels and tissue have been proposed [19–23], but the DIVA implementation described by Kotte et al. [24,25] and validated by Van Leeuwen et al. [26,27] is the only one that connects vessel tree and 3D FDTD thermal modeling.

An integral part of DIVA is an accurate segmentation of the vessels. Recent developments in angiography made clinical imaging of vessels as small as 0.5 mm possible using computed tomography (CT) or magnetic resonance (MR) scanners [28,29]. For routine HTP, contrast enhanced MR angiography (CEMRA) has advantages over CT angiography since it is non-ionizing, i.e., repeatable, and more sensitive.

In this study, we investigated the impact of arteries on the prediction of the temperature distribution using DIVA modeling added to the classical PBHE thermal model. Twelve representative patients underwent CT and MR imaging of the anatomy, as well as CEMRA. The ‘graph-cut’ method, combining information from Frangi vesselness filtering and region growing, was implemented and validated against manual markers (1) to study the feasibility of auto-segmentation as required in routine HTP and (2) to obtain a detailed vessel model. Next, specific 3D anatomy and vessel tree models of 12 patients were constructed. HTP was performed for each patient and the results of PBHE and combined PBHE + DIVA modeling were compared using predicted hyperthermia treatment quality parameters.

Methods

Patient data and scan protocols

Institutional review board approval was obtained and 12 patients eligible for radiotherapy treatment in the H&N region were randomly selected and asked to participate in this study. CT and magnetic resonance imaging (MRI) scan parameters to create the 3D patient models were previously described [30]. For vasculature tree generation, both blanco (pre contrast injection) and arterial sequences were used: TE/TR: 1.93/5.71 ms, FOV 260 mm, slice thickness 0.7 mm, acquisition matrix: 384 × 256, number of slices: 128, flip angle: 25°, contrast injection: 6.5 mL gadolinium contrast agent at 2.5 mL/s, followed by 15 mL saline solution at 2.5 mL/s.

Arterial vessel tree segmentation

Pre-processing of the image data

The images were processed for noise correction using HDCS [31] and bias field correction using the N3 method [32]. To normalize the image intensities, all image histograms were matched to a reference image [33] in the dataset. The default parameters were used for HDCS as available in its ITK implementation (http://www.insight-journal.org/browse/publication/748), and the N3 method parameters were chosen as previously optimized for anatomical sequences of the H&N [34]. For histogram matching, we used 16 landmarks using histograms built with 256 bins.

Graph-cut vessel segmentation

The segmentation method was derived from our previous work [35]. Atlas-based segmentation was combined with intensity modeling using a Graph-cut optimization algorithm for minimization and regularization. Following the terminology of [35], the spatial prior model was based on Frangi vesselness filtering [36] and the intensity model was estimated over the target image using Frangi vesselness combined with region growing based segmentation.

The intensity model was built in three steps:

1. Local maxima were found in the Frangi response and used as seed for region growing.

2. Region growing segmentation was run using the arterial image as input, and seed points from step 1; the lower intensity was determined by getting the 5% level of the intensity histogram of the arterial image in the region of high vessel-probability (defined as having a Frangi response image value higher than 0.05).

3. Using this lower intensity threshold, the region growing algorithm was run using a 26-voxels 3D neighborhood model. The foreground and background intensity models were built by sampling all voxels, within <3 mm from foreground, in or outside the region growing segmentation of the vessels, respectively. The models were constructed using Parzen window estimation with a Gaussian kernel of 50 (intensity value unit).

For the graph-cut optimization, we set the association potential weight λ₁ to 1 and the spatial prior weight λ₂ to 0.5 based on visual inspection.

Vessel centerline and radius reconstruction

The MeVisLab skeletonization algorithm [37] was applied using the binary segmentation of the vessels as input. Vessel centerline points and vessel radii per point were obtained and, using a distance based sorting, we determined the vessels' points connectivity, i.e., which points are starting points, end points, bifurcation and normal vessel points. In Figure 1, we show an example of the skeletonization algorithm (for each point, the vessel radius is available) and the same points with connectivity information. Each vessel tract was defined as the sequence of points from a starting/bifurcation point to an end point/bifurcation point as shown in Figure 1. For each tract, we estimated the radius as the average radius of all vessels centerline point in the tract and we use this information to build the vasculature model for thermal simulation.

Figure 1. Left: vessels skeleton and right: skeleton and connectivity. In this example, green points are starting points, red points are end points and purple points are bifurcations. The color of the squares around the points indicates that the point belongs to a separate tract of the vasculature.

Electromagnetic modeling and optimization

Patient models were generated from automatic segmentation of CT and MR scans [12,13] and positioning inside the HYPERcollar3D was done by an experienced technician. Electromagnetic distributions were computed using Sim4Life (v.4.0.0.2832, Zurich MedTech, Zurich, Switzerland). The resulting 3D electric field distributions were imported into the treatment planning and guidance software VEDO [38] and amplitude and phase settings per antenna were optimized to maximize the target-to-hotspot energy deposition quotient [38]. EM tissue parameters at 434 MHz are given in Table 1 and are from [39].

Table 1. EM and thermal simulation properties of tissues.

Tissue	εr	σ	ρ	c	k	Q	w
Internal air	1.00	0.0	1.2	–	–	–	–
Bone	13.07	0.09	1908	1313	0.32	0.15	10
Brainstem	55.11	1.05	1045.5	3630	0.51	11.4	558.6
Cartilage	45.14	0.60	1099.5	3568	0.49	0.54	35
Cerebellum	55.11	1.05	1045.5	3696	0.55	15.5	763.3
Cerebrum	56.81	0.75	1044.5	3653	0.51	15.7	770
CSF	70.63	2.26	1007	4096	0.57	0.0	763.3
Fat	11.59	0.08	911	2348	0.50	0.51	255
Grey matter	41.66	0.45	1041	3696	0.55	15.54	764
Lens	37.29	0.38	1075.5	3133	0.43	–	–
Lung	23.58	0.38	394	–	–	–	–
Muscle	56.87	0.80	1090.4	3421	0.4	0.96	442.8
Myelum	35.04	0.46	1075	3630	0.51	2.48	160.3
Optical nerve	35.04	0.46	1075	3613	0.49	2.48	160.3
Sclera	57.37	1.01	1032	4200	0.58	5.89	380
Thyroid gland	61.33	0.89	1050	3609	0.52	87.1	5624.3
Tumor	59.00	0.89	1050	3950	1.5	–	848
Vitreous humor	69.00	1.53	1004.5	4047	0.59	–	–
White matter	56.81	0.75	1044.5	3583	0.48	4.32	212

ε_r: relative permittivity; σ: conductivity (S/m); ρ: density (kg/m³); c: specific heat capacity (J/kg/°C); k: thermal conductivity (W/m/°C); Q: metabolic heat generation rate (W/kg); w: perfusion rate (mL/min/kg).

Thermal modeling

Transient 3D temperature distributions were calculated using the PBHE and the combined PBHE + DIVA solvers in Sim4Life. Verification of the Sim4Life DIVA Solver is provided in the Appendix. A 1 mm uniform grid was used in all simulations.

Models were simulated for 900 s to ensure that steady state is reached. The tissue thermal parameters given in Table 1 were used for both models. The listed thermal properties were taken from [40], except for tumor perfusion and thermal conductivity in muscle, fat and tumor, which were taken from [15]. The effect of thermoregulation was not modeled and values for constant perfusion were used. The SAR level, i.e., total input power, was increased until the maximum temperature in the healthy tissue reached 44 °C in the PBHE model. For the first analysis, an identical power level was used in the combined PBHE + DIVA modeling to compare the maximum predicted temperature in healthy tissue. Subsequently, the total input power of the PBHE + DIVA simulation was increased to reach a maximum healthy temperature of 44 °C. The patient initial temperature was set to 37 °C and mixed boundary conditions were applied using the following values for the heat transfer coefficients (h) and outside temperature (T): tissue – background (h = 8 (Wm^–2 °C^–1), T = 20 °C, [41]), tissue – internal air/lungs/(dental) metal implants (h = 50 (Wm^–2 °C^–1), T = 37 °C [41]), tissue – water bolus (h = 292 (Wm^–2 °C^–1), T = 30 °C; derived from in house measurement involving the device’s water bolus).

For DIVA simulations, the inflow temperature of the blood was set to 37 °C. For every vessel, the following settings were used: bucket density = 1000 m^–1, Nusselt number = 3.66, blood initial temperature = 37 °C, blood thermal conductivity = 0.52 Wm^–1 °C^–1, blood heat capacity = 4.05 × 106 Jm^–3 °C^–1. The inflow rate values for the arterial tree were fixed as in [42], i.e., 275 mL min^–1 for the two internal carotid arteries (ICAs) and 90 mL min^–1 for the two vertebral arteries (VAs). To assess the robustness against inflow rate variations, they were modified to represent high (ICA 325 mL min^–1, VA 108 mL min^–1) and low (ICA 225 mL min^–1, VA 72 mL min^–1) flow rates. At bifurcations, the flow rates were distributed over ‘child’ vessels proportional to the cubic ratio of their diameters according to Murray’s law [43].

Each DIVA simulation required 30MCells and took on average 4 h at a standard desktop computer with i5-3550 processor.

Evaluation parameters

The results of the graph-cut vessels segmentation method were benchmarked against manual annotations regarding sensitivity, specificity, accuracy and precision: $$\begin{array}{cc}\text{Sensitivity}=\frac{\text{TP}}{\text{TP}+\text{FN}}\hfill & \text{Specificity}=\frac{\text{TN}}{\text{TN}+\text{FP}}\hfill \\ \text{Accuracy}=\frac{\text{TP}+\text{TN}}{\text{TP}+\text{TN}+\text{FP}+\text{FN}}\hfill & \text{Precision}=\frac{\text{TP}}{\text{TP}+\text{FP}}\hfill \end{array}$$

As ground truth, vessel and background (no vessel) landmarks points were manually placed in pairs by a medical student. In total, 1760 vessel and 5585 background points were placed axially distributed over all patients and over 10 slices per patient. Every vessel landmark is matched with a background landmark near the vessel surface to be sensitive to mis-segmentation.

The impact of the extension of PBHE modeling with DIVA on the temperature simulation was analyzed using the hyperthermia quality parameters T90, T50 and T20, i.e., the iso-temperature levels encompassing 90%, 50% and 20% of the CTV, respectively. A two sample t-test was used to test for statistical significance.

Results

In Figure 2, a visual example of vessel segmentation by the graph-cut method is shown. On the left, the detected voxels (highlighted in purple) are shown overlaid on coronal MRA image. On the right, the resulting full arterial tree ready to be inserted in the thermal simulation, is visualized.

Figure 2. An example result of the graph-cut vessel segmentation algorithm. Pre-processed coronal slice of the MRA image of a patient with vessels segmentation color overlay in purple (left), and the skeletonization of the arterial vessel tree (right).

The sensitivity of the graph-cut segmentation method was 0.85, and specificity was 0.97. These high sensitivity and specificity values are desired to correctly identify the vessel locations. The method also provided high accuracy (0.94) and precision (0.9), which is very desirable since the thermal modeling is severely affected by false positives and false negatives. Hence, the good all-round performance of the method and visual inspection suggest that it is sufficient for thermal modeling.

Figure 3 shows the maximum temperature reached in the healthy tissue with the DIVA model when the same power as in the PBHE is used in the PBHE + DIVA model. In DIVA + PBHE models for five patients (patients 1, 7, 9, 11 and 12), the achieved maximum temperatures in the healthy tissue (using the same power levels as in the corresponding PBHE models) decreased by more than 0.2 °C. The maximum change in the peak temperature was observed in patient 9, where the maximum temperature was 1.4 °C lower than in the PBHE model.

Figure 3. Maximum temperatures in the healthy tissue reached in PBHE + DIVA models with the same input powers as in PBHE models. The maximum temperature in the healthy tissue was 44 °C in the PBHE model.

Figure 4 shows the difference in achieved HTP quality parameters in the CTV. T50 of 5/12 patients decreased notably (|ΔT50| ≥ 0.20 °C) in the presence of arterial cooling. The T50 values in the DIVA simulations were on average 0.30 °C lower, but 1/12 patient showed an insignificant increase in T50 (ΔT50 = 0.05 °C). In two cases, where the heating target-volume was in the vicinity of a blood vessel, cooling was increased by up to 1 °C. There was no clinically relevant difference in the lowest target temperature indicator, i.e., average ΔT90 = −0.16 °C. The maximum temperature surrogate T20 showed the biggest change between the models (average ΔT20 = −0.33 °C). The maximum drop in predicted T50 was 1.00 °C, where the temperature distributions for both models are illustrated in Figure 5. The vessels pass through the target area and create cold tracks in their paths. Still, our statistical analysis did not show a significant reduction in HTP quality parameters (T90: p = .64, T50: p = .45, T20: p = .38) for the total of the 12 patients. The robustness study regarding variations in the flow rates produced changes below 0.01 °C in all cases.

Figure 4. Differences in HTP quality parameters in the CTV (a) ΔT90 (|ΔT90avg| = 0.16 °C), (b) ΔT50 (|ΔT50avg| = 0.30 °C) and (c) ΔT20 (|ΔT20avg| = 0.33 °C). ‘*’ next to the patient number denotes a target volume that includes a part of the vessel tree.

Figure 5. Example temperature distribution maps overlaid on CT images for PBHE (left), and DIVA (right).

Discussion

In this paper, we report our analysis on the effect of DIVA in the H&N region. We found a notable drop in treatment quality parameters with the arterial DIVA model compared to the PBHE model. Although the differences in the values are moderate on average (0.30 °C), for three patients (patients 6, 7 and 8), a large drop in the achieved temperature metrics was observed (average drop 0.92 °C), presumably because the target volumes in these patients are significantly exposed to traversing parts of the arterial trees. Exclusion of these three patients brings the average difference in HTP quality parameters to 0.09 °C. This fact highlights the importance of DIVA modeling in specific cases when the target volume is in the vicinity of the vessel tree.

To our knowledge, these are the first reported results on the use of a clinical workflow to implement patient specific thermal modeling of DIVA in HTP including acquisition of vasculature data. Other authors did study the impact of DIVA modeling for other scenarios or using a non-patient specific vessel network. van Lier et al. [44] evaluated and compared PBHE and PBHE + DIVA models regarding the maximum temperature increase in the head after exposure to a 300 MHz radiofrequency field induced by MRI coils. They reported 0.5 °C difference in the maximum temperature increase predictions. Our approach of focusing on target volume prevents a comparison with their results, which were obtained using a more homogenous B1 field that is desirable in MRI. Our findings confirm that PBHE only predictions can overestimate the temperature increase. For exposure by a mobile phone, Flyckt et al. [45] also reported lower maximum temperature rise with a DIVA model of the eye. Using data from a healthy volunteer, supplemented by a vessel network from a previous study [18], Kok et al. [46] studied the differences in HTP quality parameters in the pelvis region when using DIVA modeling. They saw 0.2 °C, 0.4 °C and 0.6 °C decrease in T90, T50 and T10. The higher difference they found can be explained by the fact that their maximum allowed temperature in the healthy tissue was 45 °C, which is higher than in our clinical protocol.

The graph-cut segmentation is a tradeoff between region growing and Frangi vesselness based segmentation. Frangi has a low sensitivity because the vessel boundaries are not accurately segmented due to the effect of the Gaussian smoothing. Conversely, region growing segmentation tends to 'leak' including high intensity background region as foreground. Graph-cut segmentation provides better boundaries without including spurious high intensity regions in the vessel segmentation.

There are some limitations regarding the accuracy of thermal modeling in this study. Arteries with diameters larger than 0.5 mm were captured in the MRA images. Arteries with diameters in the range of 0.1–0.5 mm may affect the accuracy of PBHE + DIVA models [47]. This resolution is not achievable with current clinically available imaging technologies. Furthermore, the arteries reconstructed in this study (diameter >0.5 mm), nearly always are accompanied by a counter-current vein [48]. When this venous flow originates from a heated volume, counter-current heat exchange will increase the temperature of arterial flow (venous rewarming). A correction coefficient has been suggested to account for this effect [22]. If the venous flow originates from a cold volume; there will be no venous rewarming. As the target volume in the H&N hyperthermia is only about 250 mL, i.e., much smaller than in pelvic hyperthermia, counter-current venous rewarming is not expected. Since venous flow is much slower than arterial flow, very little independent cooling by the veins is expected [47]. Thus, we regard the effect of arteries reported in this study a good estimate of the minimum difference that can be expected in the clinic.

In vivo validation of the vasculature effect on H&N hyperthermia treatment quality has some challenges. The concept of DIVA thermal modeling was validated ex vivo by Raaymakers et al. [49] using a perfused bovine tongue. A specific validation for the H&N region would require in vivo temperature measurements near vessels. During H&N hyperthermia, placement of invasive temperature probes is challenging and risky and therefore often not feasible [50]. In addition, the thermo-sensors are positioned at a safe distance from the vessels to avoid complications like artery puncturing [51]. Hence, noninvasive measurements like magnetic resonance temperature imaging (MRTI) hold the greatest promise for validation of the DIVA modeling. However, the accuracy of MRTI strongly varies depending on the distance to distortions and confounders. Hence, a validation study requires the most advanced MRTI techniques, dedicated MR-hyperthermia equipment and careful planning.

Although the HTP quality parameters are commonly used as surrogates and have been shown to correlate with the treatment outcome [6,52], they are not sensitive to minimum temperatures. Hence, local vessel cooling may create small under-dosed areas where hyperthermia is less effective potentially leading to tumor recurrences. We hypothesize that the correlation between macroscopic quality parameters and treatment outcome is an indicator that ionizing radiation alone is sufficiently effective in those regions, which are highly oxygenized regions. An analysis on the thermal enhancement ratio [53] in these regions might provide different results. However, more research on the biology of hyperthermia on a local scale is needed for definite answers.

Inclusion of DIVA into the HTP scheme is necessary for correct thermal dose estimation. Additionally, to reach maximum allowed temperature limits in the healthy tissue, we showed that power on average needs to be increased by more than 6% when considering DIVA. Failing to deliver the correct power levels can result in a decrease in T50 by on average 0.45 °C. However, the benefits of patient specific vasculature information come with additional burden to the patient. Patients have to stay in the scanner for another 10 minutes. In addition, although fast DIVA simulation tools have been developed [46], these are not commercially available and only predict temperatures in steady state. Transient implementations of the DIVA model involve long simulation times, which makes the integration of DIVA into complaint adaptive online HTP unrealistic [38].

Conclusion

Our analysis, based on 12 patient models, showed no significant decrease of the target temperature isopercentiles when including the effect of patient specific vasculature information into the PBHE model. Nevertheless, in specific cases, where major vasculature traverses the strongly heated target volume, a clinically relevant difference was observed. Based on these findings about the impact of the vasculature above 0.5mm in diameter, we advise to consider DIVA modeling if such vessels are passing through the heated region. For the vessel tree segmentation and extraction from CEMRA, the graph-cut method has good automation potential: it performs well in terms of accuracy, precision and specificity without losing on sensitivity.

Disclosure statement

No potential conflict of interest was reported by the authors.

Additional information

Funding

This study was financially supported by the Technology Foundation STW [Grant 10846], the Dutch Cancer Society [Grant EMCR 2009-4270] and COST Action BM1309 [COST EMF-MED].

Appendix. Verification of Sim4Life DIVA solver

The Sim4Life (and SEMCAD X) DIVA solver is a reimplementation of the approach presented by Kotte et al. [24,25]. The original terminology introduced by Kotte is used throughout this Appendix. 1D vascular network simulations of blood-flow heat transport are coupled with spatio-temporal 4D thermal simulations (based on the Pennes bioheat equation (PBHE) [16]), whereby energy exchange resulting from temperature differences between the vascular blood temperature and its tissue environment is considered. For that purpose, the vascular network is discretized into buckets, while the anatomical domain is voxeled (rectilinear mesh). The heat exchange between the two simulations is realized by combining the temperature in voxels neighboring the vessel bucket with the temperature in the bucket to fit an analytical solution of the radial temperature dependence and to thus derive the energy flux through the bucket surface. The Sim4Life implementation goes beyond the original implementation by allowing for adaptive (inhomogeneous, but still rectilinear) gridding and voxeling. For more information, see the original work of Kotte.

Within this paper, a reduced set of functionalities from the DIVA model is employed: flow rates are defined at in-flows and the division of flow rates at branching points is specified. However, no blood uptake from the tissue along the vasculature is modeled, nor is the blood considered to flow into the tissue at endings of the vascular network, or elsewhere. This simplification is acceptable, as the heating is localized and vessels are starting and ending well outside the volume-of-interest. Therefore, the following aspects need to be verified: transport within the vascular network, energy exchange with the 4D thermal simulation, as well as correct handling of inflow- and branching-points. As the Sim4Life solver is re-implementing the original model from Kotte, the original validation work [17,20,24,25] does not need to be re-performed, provided the verification successfully demonstrates that the implementation was performed correctly. This is the purpose of this Appendix.

The most important benchmark that will be used, is the case of a straight vessel along the axis of a tissue cylinder (Figure 6), for which an analytical formula for the thermal equilibrium length (Leq) has been derived in [24]. Unless specified differently, the simulation parameters provided in Table 2 are used, in line with the original publication.

Figure 6. Principal verification benchmark: a coaxial setup involving a straight vessel inside a tissue cylinder.

Table 2. Basic model parameters (in accordance with Table 1 in [24]).

Tissue conductivity	0.6 W·m^–1·K^–1
Tissue specific heat	4000 J·kg^–1·K^–1
Tissue density	1000 kg·m^–3
Tissue radius	0.003 m
Tissue length	0.01 m
Vessel radius	0.00025 m
Vessel length	0.01 m
Blood velocity	0.025 m·s^–1
Tissue initial temperature	1 °C
Inflow blood temperature	0 °C
Bucket density	5000
Number of grid voxels	31 × 31 × 50

Implementation correctness and convergence with discretization refinement (bucket density): Figure 7(a) shows the theoretical and simulated temperature along the vessel for different bucket densities, while Figure 7(b–d) shows the local equilibration length along the vessel as computed from the simulation results (using Equation (11) from [24]) in comparison to the theoretical equilibration length. In Figure 7(a), blood temperature along the vessel for theoretically calculated and simulated results for different bucket densities are reported. The simulated results show agreement with the theoretical results. As expected, the agreement gets better with increasing bucket density. With increasing bucket density, the computed equilibration lengths approach the theoretical value, with the strongest deviations seen at the in- and out-flow (as expected). More importantly, these results are identical to those shown in Figure 4 from [24], demonstrating the correctness (or equivalence) of the Sim4Life implementation of in- flow, energy exchange, convective transport and 4D thermal simulation.

Figure 7. (a) Theoretical and simulated temperature along the vessel for different bucket densities. Theoretical and simulated L_eq for a bucket density of (b) 5000, (c) 2500 and (d) 1250.

Vessel radius dependence: The vessel radius dependences of the theoretical and numerically obtained equilibration lengths are shown in Figure 8, further verifying the solver. Again, these results are in agreement with those shown in Figure 6 from [24].

Figure 8. L_eq as a function of the vessel radius.

Connected segments and branching: It remains to be demonstrated that connections between vascular segments, as well as branching points are correctly handled. First, connections have been tested by connecting three similarly sized segments back to back. The computed vascular temperature along the connected branches are shown in Figure 9 along with the theoretical values. Small deviations can be seen at the connection points, which are related to how the exchange and estimation voxels are placed and are more apparent when looking at the derived equilibration lengths (sensitivity of derivative). However, the thermal impact is minor and simulation results are close to the analytical solution.

Figure 9. Temperature profile along three connected vessel segments (computed and theoretical; left) and L_eq (right).

To test the implementation of branching, another benchmark was constructed. This benchmark is qualitative, as no analytical solution is known. A parent vessel goes through two branching levels. Both, symmetric and asymmetric branching (radiuses of child vessels) is modeled. Flow is distributed according to Murray’s law. In Figure 10, the resulting tissue temperature distribution is shown. To estimate the four expected temperature at the midpoint of 2nd generation child vessels, the analytical solution for the cylindrical case is used in each segment (with the slab thickness as cylinder radius). The computational results agree within 0.4–6% with the estimated values.

Figure 10. DIVA benchmark illustrating branching with symmetric and asymmetric flow-rate division. Red straight lines are indicating the vessel centerlines.