Temperature-Controlled Hyperthermia with Non-Invasive Temperature Monitoring through Speed of Sound Imaging

Abstract

Hyperthermia therapy (HT) is used to treat diseases through heating of high temperature usually in conjunction with some other medical therapeutics such as chemotherapy and radiotherapy. In this study, we propose a promising temperature-controlled hyperthermia method that uses high-intensity focused ultrasound (HIFU) for clinical tumor treatment combined with diagnostic ultrasound image guidance and non-invasive temperature monitoring through speed of sound (SOS) imaging. HIFU heating is realized by a ring ultrasound transducer array with 256 elements. In this study, tumors in the human thigh were set as heating targets. The inner structure information of thigh tissue is obtained by B-mode ultrasound imaging. Since the relationship between temperature and SOS in different human tissue is available, the temperature detection is converted to the SOS detection obtained by the full-wave inversion (FWI) method. Simulation results show that our model can achieve expected hyperthermia of constant temperature on tumor target with 0.2 °C maximum temperature fluctuation for 5 h. Through simulation, our proposed thermal therapy model achieves accurate temperature control of ±0.2 °C in human thigh tumors, which verifies the feasibility of the proposed temperature-controlled hyperthermia model. Furthermore, the temperature measurement can share the same ring ultrasound transducer array for HIFU heating and B-mode ultrasound imaging, which provides a guiding significance for clinical application.

1. Introduction

HYPERTHERMIA therapy (HT) is used to treat diseases through maintaining high-temperature heating of the whole body or a part for a long period, which is usually used alone or in conjunction with some other medical therapeutics such as chemotherapy and radiotherapy [1]. Many successful clinical trials have applied the therapy to treat different types of cancer tumors, such as recurrent breast cancer, liver cancer, uterine fibroids and bladder cancer [2,3,4,5]. HT is mainly divided into three categories: local HT, regional HT and whole-body HT [6]. The purpose of local HT is to increase the temperature of the tumor without affecting the surrounding normal tissue. Generally, local HT is applied to treat tumors below the skin or in body cavities near the surface of the body through using external or interstitial heating modalities [7]. In contrast, regional HT is more complicated than local HT because its targets are deep-seated tumors in body cavities, especially locally advanced tumors. Non-invasive methods to realize regional HT include heating the blood, irrigating the body cavities and radiofrequency (RF). Whole-body HT is usually applied to treat metastatic cancer by using surface heating, radiation induction or extracorporeal induction. During the entire hyperthermia treatment period, the temperature of the tumor and surrounding tissues needs to be strictly controlled at the expected temperature. Among many techniques such as microwave, infrared and RF, high-intensity focused ultrasound (HIFU) is increasingly being developed as a promising modality to induce hyperthermia for clinical tumor treatment combined with diagnostic ultrasound guidance [8,9].

For many kinds of benign and malignant tumors, HIFU plays an auxiliary role in the current standard treatments such as surgery, gene therapy and radiation. It was first reported by Lynn et al. in the early 1940s that HIFU can pass through tissue without harm and cause tissue ablation in a focal area. In the 1950s and 1960s, significant early work including high-intensity ultrasound and ultrasound visualization was completed. By the late 1990s, commercially available ultrasound-guided HIFU devices had been developed [10]. Extensive advances in phased-array applicator technology have enabled HIFU to perform better in focusing [11]. With the advances in imaging, the accurate focus position makes clinical application possible [12]. Magnetic resonance imaging (MRI) and ultrasound imaging are two guiding modalities with complementary guiding features [13,14,15] which have few side effects in humans. HIFU uses the characteristic that ultrasound can penetrate the human body without damage and focus on tumors in vivo causing a thermal effect (main effect), cavitation effect and mechanical effect. Since the high temperature generated by thermal effect on the focus can destroy the cell structure and damage the protein, it will lead to the coagulation necrosis of the tumor target [16,17,18] instead of burning the surrounding normal tissue. The sensor emits high-intensity ultrasound signals to the tissue target, and then the tissue absorbs ultrasound energy and converts energy into thermal energy during therapy. This heat deposition can cause the rapid temperature rise of tissue. Due to the focusing of the ultrasound beam, thermal energy is mainly concentrated on the tumor target, improving the temperature to the desired range (e.g., 5 °C–6 °C higher than the human body temperature). Therefore, there is almost no obvious heat deposition on the surrounding tissue, which avoids the unpredictable growth of healthy tissue caused by excessive temperature [19].

There are three typical types of HIFU equipment systems that have been put into clinical practice for treatment of tumors: HIFU systems guided by MRI (MRgHIFU), HIFU systems guided by ultrasound (USgHIFU) and HIFU systems guided by both MRI and ultrasound. During the MRgHIFU treatment, MRI equipped with a thermal mapping system can make tumors visible clearly and monitor the thermal ablation effect of target tissue based on the sensitivity to temperature changes in real time [20,21,22,23]. The high spatial resolution of MRgHIFU enables the accuracy of image guidance, which offers abundant anatomic details for tumor detection [24]. However, USgHIFU is a more affordable imaging method and has many successful clinical applications. In addition to its high temporal and spatial resolution, USgHIFU is also insensitive to geometrical distortions, and all these characteristics make ultrasound imaging suitable for guiding HIFU therapy [25]. Considering the advantages of imaging capabilities combining MRI with ultrasound, a new HIFU technology has been proposed, in which ultrasound imaging is used for continuous target tracking and real-time MR-thermometry is performed to guide HIFU heating simultaneously [26]. However, MRgHIFU needs large and expensive systems in comparison with USgHIFU, which is relatively cheap and portable [27]. When USgHIFU systems are equipped with specially designed transducers, it is possible to realize real-time visual temperature monitoring of the target soft tissue [28].

It is hard to realize long-time constant temperature hyperthermia and control the temperature with little fluctuation [29] because of the lack of precisely adjustable heating and temperature measurement techniques [30]. In our study, we aim to simplify the hyperthermia system by using a ring array ultrasound transducer to perform two major functions of HIFU heating and temperature measurement in human tissue.

2. Mathematical Model and Theory

2.1. Acoustic Heating

Acoustic waves are the propagation forms of sound, which are essentially the transmission of energy in the medium. When acoustic waves propagate, the classic case is to pass through a lossless fluid medium with homogeneous density that satisfies the simple conservation of momentum, mass and energy [31]. However, the reality is that the acoustic medium is heterogeneous, and its SOS and ambient density are different. At the same time, the propagation of acoustic waves cannot ignore the loss of acoustic energy due to random thermal motion. In addition, acoustic waves passing through the compressible medium will cause dynamic fluctuations in variables such as density, velocity of acoustic particle, pressure, temperature, etc. The vibration and friction between the acoustic waves and the medium will convert a part of the acoustic energy into thermal energy and generate acoustic absorption physically with dispersion [32].

Since we use a ring ultrasound transducer array in this study, acoustic heating is realized by the dynamic delayed emission of acoustic pulse signals for each ultrasound transducer element according to Equation (1):

$${(x_i-o_{x0})}^2+{(y_i-o_{y0})}^2=R^2,$$

$$l_i=\sqrt{{(x_i-o_x)}^2+{(y_i-o_y)}^2},$$

$$t_0=\frac{l_0}{v_s},$$

$$\mathsf{\Delta }{t}_i=t_0-\frac{l_i}{v_s},$$

(1)

where $(x_i,y_i)$ is the coordinate of a single element on the ring ultrasound transducer array with the radius of R, which has 256 elements ( i ranges from 1 to 256). $(o_{x0},o_{y0})$ is both the center of the ring ultrasound transducer array and the center of the human thigh model whereas $(o_x,o_y)$ is both the focus position and the center of the tumor. In addition, $v_s$ stands for the SOS of tissue and $l_i$ stands for the distance from each element to the focus. We take the time $t_0$ calculated from the farthest distance $l_0$ as a benchmark, which can ensure that the delay value of each element $\mathsf{\Delta }{t}_i$ is greater than or equal to 0 when transmitting signals. When we add the corresponding delay $\mathsf{\Delta }{t}_i$ to each ultrasound transducer element, all the transmitted signals can reach the focus of the target position at the same time.

The ring ultrasound transducer array model applied in this study is illustrated in Figure 1. Furthermore, the simulation model is immersed in the water bath at a constant temperature of 30 °C.

The heat conduction equation is an important partial differential equation in HIFU heating, which describes how the temperature in a region varies over time. Thermal diffusion is the main factor needing consideration during heating. In addition, Pennes’ bioheat equation [33] accounts for many aspects, such as heat deposition due to ultrasound absorption and advective heat loss due to tissue perfusion (blood flowing through the tissue). Therefore, the heat conduction equation including the factor of vascular heat perfusion is given by:

$$\frac{dT}{dt}=\frac{k}{c\rho }·{\nabla }^{2}T-\frac{B}{c\rho }(T-T_a)+\frac{Q}{c\rho },$$

(2)

where k is thermal conductivity in W/m·°C, c is specific heat capacity in J/kg·°C, $\rho$ is medium density in kg/m${}^3$. Q is the volume rate of heat deposition in W/m${}^3$, which can be treated as the initial heat source used for hyperthermia. Furthermore, B is the product of blood density (kg/m${}^3$), blood specific heat capacity (J/kg·°C) and blood perfusion rate ($s^{-1}$) whereas $T_a$ is the blood ambient temperature set to 37 °C by default. ${\nabla }^2$ denotes the Laplace operator and ${\nabla }^{2}T$ needs to consider the specific case of a function $T(x,y,z)$ of three spatial variables $(x,y,z)$ given by:

$${\nabla }^{2}T=\frac{{\partial }^{2}T}{\partial x^2}+\frac{{\partial }^{2}T}{\partial y^2}+\frac{{\partial }^{2}T}{\partial z^2},$$

(3)

where T is the temperature at the point $(x,y,z)$.

We establish a preset relationship equation between the temperature change $\mathsf{\Delta }T$ (°C) and the heat source change $\mathsf{\Delta }Q$ (W/m${}^3$) in Equation (4), which is the key to temperature-controlled hyperthermia model.

$$\mathsf{\Delta }Q=-m·\mathsf{\Delta }T,$$

$$m_i=\left\{\begin{array}{cc}{m}_{i-1}-{\alpha }_1·m_{i-1},\hfill & \mathsf{\Delta }T>0\hfill \\ m_{i-1}+{\alpha }_2·m_{i-1},\hfill & \mathsf{\Delta }T<0\hfill \end{array}\right.$$

(4)

where m is the core operator that allows the constant temperature system to modify the intensity of the heat source as the temperature changes. $m_i$ stands for the value of in the period, which is calculated from $m_{i-1}$ (the value of m in the last period). Furthermore, $\alpha$ is the control factor of m. ${\alpha }_1$ is used after the temperature has risen, and ${\alpha }_2$ is used to prevent temperature dropping. When the temperature rises during the last heating period, $\mathsf{\Delta }T$ is positive and $\mathsf{\Delta }Q$ is negative correspondingly, so it is necessary to reduce the intensity of heat source in the consequent heating period. When the temperature drops, the process is opposite. Using continuous heating by adjusting the intensity of heat source ensures that the temperature fluctuates up and down around the target temperature. To achieve the objective that the temperature fluctuation is within the desired range of the target temperature, it is important to fine-tune m. We first set small upper and lower limits according to the target temperature range. If the temperature exceeds the preset upper limit, m will decrease by $\alpha 1$ times the m of the previous period. The increase by $\alpha 2$ times the m of the previous period occurs when the temperature exceeds the preset lower limit. The increasing thermal diffusion over time requires the control factor $\alpha$ to strengthen the control of m. Therefore, as long as the temperature exceeds the upper or lower limits, $\alpha 1$ and $\alpha 2$ will respectively add 0.01 to their own current values. Through combining the above equations, we can achieve effective and accurate temperature control where the temperature fluctuation is within little range by the frequent iteration of core operator m.

2.2. Speed of Sound Imaging

SOS imaging is a non-invasive imaging modality that can show the acoustic property change in body tissue [34]. Previous studies have shown that the SOS in body tissue is sensitive to the temperature change [35]. In our study, the temperature control is necessary for the constant temperature hyperthermia system. Due to the available relationship between temperature and SOS (temperature and SOS vary in direct proportion) [36], the temperature monitoring in the medium is realized by the SOS imaging.

There are basically three categories of SOS imaging methods. The first one is filtered back projection (FBP) reconstruction technique, which is also one of the main techniques in image reconstruction of clinical CT applications [37]. By employing the Radon transform, FBP technique can complete the image reconstruction in frequency domain, which avoids the direct solution of complex inverse problem including finding the mapping function that matches the measured projection data and the SOS distribution. In this way, the FBP reconstruction technique can provide fast computing speed with a reduction of computation complexity by using the Radon transform, which is the foundation of the FBP algorithm. A high-quality image can be produced when enough transmission data have been acquired, which is difficult to realize in practical application. Furthermore, FBP technique assumes that the propagation of acoustic waves travels along straight lines in body tissue [38]. The approximation causes the limitation of the precision in SOS reconstruction. The second one is bent-ray tracing method (BRTM) based on ray theory, which employs time of flight (TOF) of acoustic wave propagation paths to achieve the SOS distribution reconstruction [39]. In this method, the acoustic wave propagation paths can be estimated by using the principles of geometrical optics and these paths are not straight when traveling through different medium interfaces [40]. The precision of the TOF measurement is crucial to the accuracy of the SOS distribution estimation. In our application, the average length of acoustic wave propagation paths is around 120 and the average SOS is about 1580. Assuming that the sampling frequency of the ring ultrasound transducer array is 20, the average error of SOS is more than 1, which means that the deviation of temperature monitoring is not satisfied with the temperature fluctuation limitation. The third one is full-wave inversion (FWI) based on the theory of acoustic wave propagation, which has many applications in medical imaging and geophysics [41,42,43]. When applying the FWI method, the SOS image reconstruction is realized by estimating the SOS distribution using the acoustic wave equation and the measured data. The measured acoustic wave data can be acquired by the receiver ultrasound sensor elements and the estimated acoustic wave data can be acquired through pseudospectral k-space method. Given the gradient of a sum of squared norms between the measured data and the estimated data, a gradient descent-based optimization algorithm can be employed to update the SOS estimate at every iteration [44].

Considering high order diffraction and scattering, FWI method can provide the better image quality in robustness, resolution and accuracy compared with two previous methods that have an approximate hypothesis [45]. In our work, temperature monitoring in the soft tissue is realized by the FWI method, which depends on the optimization iteration algorithm to achieve the SOS distribution reconstruction. At the beginning of iterations, the initial assumed SOS values of the medium play an important role in calculation speed and accuracy of SOS image reconstruction. Through a correspondent B-mode ultrasound image, the boundary information and position information of human thigh tissues can be acquired. Combined with cross-sectional anatomy images of human thigh and the acoustic properties of the human body, the precise initial values of SOS can be employed as a priori information to estimate the actual SOS values, which is helpful to improve the calculation efficiency.

3. Experiment and Analysis

3.1. Simulation Tools

Our simulation mainly relies on K-Wave to build the model and perform HIFU heating. K-Wave as an open-source toolbox can provide users with a MATLAB script interface to simulate the time-domain propagation of acoustic waves in one-dimensional, two-dimensional or three-dimensional space. In addition to explaining linear and nonlinear acoustic propagation, K-Wave toolbox also has various functions including the calculation of the distribution of heterogeneous material parameters and power law acoustic absorption. More details about the K-Wave model including its main governing equations and the boundary conditions used to our simulation can be found in the K-Wave user manual [46,47,48].

3.2. Simulation Setup

Our simulation is performed in a three-dimensional space supported by K-Wave. 3D simulation space is defined as $N_x\times N_y\times N_z$ ($490\times 490\times 40$) in grid size, in which the spacing of each grid point is about 0.449 mm. The actual physical size of the simulation model is 220 × 220 × 18 mm. By observing Figure 2a, we can find that the internal structures of the human thigh including: skin, subcutaneous fat, muscle, thigh bone, bone marrow, tumor and blood vessels. All the tissues are defined respectively with reference to the properties such as size, SOS, density, attenuation coefficient, thermal conductivity and heat capacity in

Table 1. Properties of medium.Size describes the internal structure in the simulation model we used. SOS, Density and Attenuation Coefficients were used to calculate the acoustic field, and Thermal Conductivity and Specific Heat Capacity were used to calculate the biological heat transfer.

Medium	Size (mm)	SOS (m/s)	Density (kg/m${}^3$)	Attenuation Coefficient (NP/Mhz·m)	Thermal Conductivity (W/m · °C)	Specific Heat Capacity (J/kg · °C)
Water	220 × 220 × 18	1482	1000	0.025	0.60	4178
Skin	1.70 (Thickness)	1595	1109	21.158	0.37	3391
Fat	10.00 (Radius)	1430	911	4.358	0.21	2348
Muscle	60.00 (Radius)	1580	1090	7.109	0.49	3421
Bone	5.80 (Radius)	2198	1178	47.000	0.32	1313
Bone Marrow	6.70 (Radius)	1372	980	4.358	0.20	2065
Tumor	6.50 (Radius)	1450	1050	5.650	0.51	3540
Blood Vessel Wall	0.50 (Thickness)	1570	1100	7.020	0.46	3306
Big Vessel (Blood)	2.90 (Radius)	1550	1060	2.368	0.52	3617
small Vessel (Blood)	0.65 (Radius)	1550	1060	2.368	0.52	3617

[49]. In addition, the blood perfusion rate is set to 0.01 s${}^{-1}$. In our study, we assume the human thigh is placed in the ring ultrasound transducer array with 256 elements (with 100 mm radius and 9 mm height) as shown in Figure 1. The whole thigh is immersed in water at a constant temperature of 30 °C. All the ultrasound elements transmit pulse signals with five cycles. In addition, the center frequency of the transmitted signals is 1.5 Mhz whereas the sampling frequency is 12 MHz. Therefore, the wavelength is approximately 1mm calculated from the center frequency and SOS. Considering that the height of the ring ultrasound transducer array is 9 mm, the focal shape should be a linear acoustic source composed of nine elements instead of one single acoustic source. Figure 3a illustrates the focus of ultrasound signals emitted by three of the 256 elements. The gray level represents the intensity of sound pressure and the focus spot is in black. It is illustrated that the profiles of ultrasound signals generated by 2 of 256 elements in Figure 3b. Figure 3c shows the pressure field in Pa generated by the ring ultrasound transducer array with 256 elements.

To confirm the universality of the simulation, we set three human thigh models with different muscle radius including 55.00 mm, 60.00 mm and 65.00 mm. Furthermore, we change the position of the tumor according to the growth characteristic that tumors may grow close to blood vessels because they rely on vessels to provide a lot of oxygen and nutrients [50]. The multiple models of human thigh referring to Figure 2a are given in Figure 4, which displays the density maps of different human thigh models.

3.3. Imaging

Figure 2a as a cross-sectional anatomy image provides a preliminary reference [51] for us to understand the human thigh structure. In order to build a simulation model with more details, the inner structure information of human thigh needs to be obtained by B-mode ultrasound imaging, especially the details of the tumor region.

The B-mode ultrasound image of human thigh tissue has been constructed by the Delay and Sum (DAS) [52] algorithm, which is supported by the K-Wave simulation tool with ring ultrasound transducer array, as shown in Figure 2b. B-mode ultrasound image helps us to find the focus position for HIFU heating and provides us with the a priori information of the internal medium in SOS imaging. In the process of SOS image reconstruction, we first split up the human thigh model into seven parts including skin, subcutaneous fat, muscle, thigh bone, bone marrow, tumor and blood vessels with the boundary information from the B-mode ultrasound image. Due to the existence of thermal diffusion, the temperature gradually decreases along the direction of thermal diffusion with the focusing position as the center. Therefore, the SOS in each part is not consistent. The simplification model of SOS distribution reconstruction is established using the average value of each part in the actual SOS image as the value of the target SOS in the corresponding part. In the iterative estimation of SOS distribution with FWI method, the SOS in each part is unrelated with other parts and the initial SOS in each part refers to the acoustic properties of the corresponding human tissue. Furthermore, we set the SOS in thigh bone, bone marrow and blood vessels to the standard values (SOS values at room temperature) and keep the SOS in these regions constant. The reasons for the constant SOS in the above regions are as follows. Firstly, since the regions of bone and bone marrow are outside the effective heating area, HIFU heating almost has little influence in the two regions. Secondly, considering the size of vessel regions, the SOS change in the small region hardly contributes to the accuracy of the whole SOS iteration. Meanwhile, the SOS iteration in the small region is sensitive to noise which may lead to the wrong iteration direction of SOS, so it is necessary that SOS in vessel regions remains constant. These assumptions help us to improve the accuracy in estimating the SOS values and the reconstruction image of SOS is shown in Figure 2c. Since the accurate SOS in the tumor region is essential to the temperature-controlled hyperthermia, the SOS in the internal region of the tumor is averaged to improve the robustness and accuracy of SOS estimation. During the SOS imaging iteration, the averaged SOS value in tumor region is illustrated in Figure 5. From Figure 5b (a zooming curve of the black rectangle region in Figure 5a); we can discern that the SOS fluctuation in the tumor region can be limited to less than 0.2 m/s which in return can yield the temperature monitoring accuracy to 0.2 °C.

3.4. Temperature-Controlled Hyperthermia

Through the focusing process of HIFU, the ultrasound beams with high intensity are focused in the center of the tumor to generate acoustic pressure field that demonstrates a good focusing effect. Long-term constant temperature hyperthermia can be divided into two heating stages. The first heating stage aims to quickly heat the human thigh tissue until the temperature of the tumor center is close to the target temperature. The second heating stage is constant temperature control realized by the combination of the heat conduction equation (as shown in Equation (2)) and the equation of temperature-controlled hyperthermia model (as shown in Equation (4)).

The key of the second heating stage is dynamical modification of the heat source according to Equation (4) every 60 s. The core operator is set to multiple initial values suitable for different $\mathsf{\Delta }T$ based on the relationshipthat and are inversely proportional (the product of m and $\mathsf{\Delta }T$ ranges from 8 × 10${}^4$ W/m${}^3$ to 2 × 10${}^5$ W/m${}^3$, which is also 0.005 to 0.01 times the amplitude of the heat source of the previous period). According to Equation (4), when the temperature of the tumor center exceeds the upper and lower limits of the target temperature, the operator will self-adjust by $\alpha$ times that of the previous period. It will self-decrease when the temperature rises and exceeds the upper limit, and it will self-increase to prevent the area of tumor cooling down when temperature drops. The initial values of ${\alpha }_1$ and ${\alpha }_2$ in our study are both set to 0.01. Since the relationship between temperature and SOS in the different human tissue is available (temperature and SOS vary in direct proportion), the temperature detection in the last heating period is converted from the SOS value which is updated accordingly.

We set target temperatures at 43 °C for each model and the temperature control curve of long-term constant temperature hyperthermia is illustrated in Figure 6. It can be observed that the temperature can be kept within ±0.2 °C range of the target temperature with slight fluctuations for 5 h. Furthermore, we are also concerned with the temperature distribution in three-dimensional perspective during the whole heating process. We select model B in Figure 4b with the target temperature of 43 °C as an observation object, using COMSOL Multiphysics${}^{\circledR }$, for simulation. As shown in Figure 7, the heat is mainly concentrated in the center of the tumor and there is almost no obvious heat deposition in other areas, which proves that our HIFU heating can achieve precise focusing on the tumor.

4. Conclusions

This paper proposes a temperature-controlled hyperthermia model for tumor therapy with non-invasive temperature measurement through SOS imaging. We use the ring ultrasound transducer array to achieve a precise and variable focus position with good focusing effect at the tumor center. In addition, the results of the temperature control process prove that we have solved the critical problem of long-time temperature-controlled hyperthermia with the temperature being kept constant. Furthermore, the temperature monitoring and the constant-temperature hyperthermia sharing the same system are designed to guide the production of corresponding medical equipment, which realizes the integration of functions and reduces the difficulty of clinical implementation. When performing the simulation, the program divides the simulation space into several grids. Due to the small size of the individual array elements of the toroidal array transducer proposed in this paper, the deviation of this discretization is small. However, in practical applications, larger curved transducers are often used to increase the power of the transducer, which may lead to greater deviations in the sound field simulation. In this case, it is necessary to reduce the size of individual grids to improve the accuracy of the simulation, but this will lead to an increase in the resources occupied by the simulation. Balancing simulation accuracy and resource occupation is an issue worthy of continued discussion.

Next, we will produce the ring array transducer as designed and gradually test it on phantom, experimental animals to verify the reliability of our proposed thermal therapy model. Since the model is cost-controllable and non-invasive, it can extend the clinical application of hyperthermia. Furthermore, the thermotherapy model has good compatibility with several human sites, and with appropriate improvements, it may be possible to extend the treatment of tumors from the thigh to other parts of the body, such as the breast, providing promising clinical applications.