Microwave thermometry with potential application in non-invasive monitoring of hyperthermia

1 Introduction

Microwave hyperthermia [17] has shown to be a good candidate to help increase therapeutic gains of traditional cancer therapies like radio- and chemo-therapy [11, 13]. It aims to increase the tumor temperature to therapeutic levels of 40–44 °C while keeping healthy tissue at the normothermic temperatures. To control and evaluate the quality of the treatment in terms of achieved temperatures in tumor as well as healthy tissues, thermal dose monitoring is a critical instance of the treatment. To this end, a robust and preferably real-time feedback control that helps to adjust the focal point in the target is warranted. Hence, temperature monitoring needs to be considered as an integral part of every reliable thermotherapy system [31]. In past decades, several methods have been proposed in the literature among which the invasive measurements in tumor-related reference points [7, 21, 27] and MR-guided thermometry [14, 10, 32] are currently in clinical use.

While use of cheap and readily accessible multipoint thermal probes with thermocouples, thermostors or fiberoptics sensors is still considered as a gold standard, the invasive nature of the procedure limits the information about temperature distribution to few monitoring positions. MR-guided thermometry provides a 3D non-invasive temperature monitoring; nevertheless, the acquisition is not fast enough to be integrated into real-time schemes [32]. Furthermore, the hybrid HT-MR systems have other limitations such as high expense, access and compatibility issues. MR is also not suitable for intraoperative monitoring in patients with MR incompatible pacemakers or artificial heart valves [9]. Microwave thermometry, is an alternative and appealing method as it could provide information about 3D thermal profiles while using the same antennas that are used for heating. Such a theranostic system would resolve the complexity and costs accompanied with MR-thermometry while providing the clinicians with sufficient information about thermal profiles during the treatment [16].

There is a tangible change in dielectric properties of water within the temperature range of interest [22]. The fact that dielectric properties of soft-tissues are mainly dominated by water fraction makes it possible to cast this problem into the framework of differential microwave imaging proposed in several publications [12, 25, 16, 8]. Following a similar methodology, we consider in this study an annular phased array system of 16 monopole antennas designed to work at the ISM band of 915 MHz. At this frequency, a 1 % dielectric change between 28 °C and 53 °C according to in vivo measurements in animal liver tissue has been reported [24]. Although the range of dielectric changes in the relevant scenarios may seem not to be very large, it is still considerably higher than other non-invasive thermal monitoring options. The challenge to reconstruct a relatively small changes is at least partly compensated by possibility to use Born or Distorted Born approximations to linearize the inherently nonlinear volume integral equation at hand.

As far as the inversion technique is concerned, the linear least square method (LLSM) is used, which has a smooth kernel here, meaning that the solution is highly sensitive to small perturbations in the data and is therefore ill-posed. To deal with the ill-posedness, a wide variety of regularization schemes have been proposed in the literature among which Tikhonov is one of the popular and efficient methods [28, 30, 29]. The regularization parameter can be chosen either by an a priori rule like Tikhonov’s regularization strategy [1, 2, 20] or by a posteriori rules like the balancing principle [18] and Morozov’s discrepancy principle [26]. For image quality restoration, an adaptive finite element method, developed in [23], is combined in this study with a balancing principle for optimal choice of the regularization parameter. The main idea of an adaptive finite element method is to obtain better image quality via local mesh refinements through minimization of Tikhonov’s functional. Our numerical investigations show that an adaptive finite element method improves image quality of reconstructed images compared with the standard least square method.

The paper is outlined as follows. Simulation setup and the standard least square method are described in Section 2, followed by a formulation of the adaptive finite element method described in Section 3. Balancing principle and fixed point algorithm are also formulated in this section. Section 4 is devoted to discussion about the results obtained. Finally, Section 5 concludes the main findings of this note.

2 Materials and methods

2.1 Forward problem and data acquisition

A cylindrical scenario has been considered with 16 monopole antennas immersed in a coupling fluid. The matching liquid is a 20 : 7 (alcohol : water) mixture chosen based on [16] to balance conductive loss, relative permittivity and antenna matching. We used CST Microwave Studio [33] to simulate the scenario whose CAD model is shown in Figure 1, along with the characteristic of the antenna which is designed to have its best performance at 915 MHz. The selection of this frequency is to achieve a good compromise between both penetration depth and resolution of the final image.

The original scenario, which has an off-centered cylindrical inclusion as target, is considered as the baseline (i.e. the initial state at t=0 to which the dynamic cooling process of the heated target will be compared at each time step). Then, to imitate the dynamic process of the cooling over a 10-minute window of time, the dielectric parameters of the target have been changed in a parametric sweep of CST according to Table 1. These parameters are also taken from [16] and adjusted accordingly to fit the scenario at hand.

Figure 1

(a) The cylindrical scenario consisting of 16 monopole antennas immersed in a matching liquid tank. The ROI is also highlighted as an orange box. (b) Return loss (S11) of the designed antenna for the frequency band of 915 MHz.

(a)

(b)

Table 1

Changes in permittivity and conductivity of the target as it starts to cool down from 55 °C to 29 °C over a ten-minute window of time.

Timeline (min)	ϵr⁢(t)	σ⁢(t) (S/m)
t=0 (baseline)	26	0.12
t=2	28	0.15
t=4	29	0.19
t=6	29.8	0.21
t=8	30.5	0.23
t=10	31	0.24
Coupling fluid 20 : 7 isopropyl : water (ϵr=24.5,σ=0.46)

3 Adaptive finite element method with balancing principle

For improvement of the solution 𝑚 obtained via minimization of functional (2.4) using formula (2.7), we will apply an adaptive finite element method (AFEM) developed in [23] combined with a balancing principle for choosing the regularization parameter 𝜆, for minimizing the functional

where F⁢(u) is the convolution operator defined similarly to [23], 𝑚 is the solution obtained via (2.7) and ∥⋅∥2 is the L2 norm. AFEM was tested on experimental data for reconstruction of a two-dimensional signal measured by microtomograph and has shown significant improvement of the image restoration quality [23].

The finite element method for minimization of Tikhonov functional (2.4) was formulated in [4]. Here, we briefly present an AFEM algorithm for minimization of functional (3.1). To find function 𝑢 which minimizes functional (3.1), we look for a stationary point of this functional with respect to 𝑢 such that Jλ′⁢(u)⁢(b)=0, where Jλ′⁢(u) is the Fréchet derivative of functional (3.1) acting on 𝑏; see details of derivation in [4].

To proceed further, we set-up the dimensionless computational domain Ω={(x,y,z)} and decompose it into two-dimensional slices Ωz for any desired value of 𝑧. Then every slice is discretized by the mesh Kh which consists of non-overlapping triangles 𝐾 satisfying the mesh regularity assumption [19] with the mesh function ℎ which is defined as a piecewise-constant function such that h=hK, where hK is the diameter of 𝐾 which we define as a longest side of 𝐾.

Further, we introduce the finite element space Vh consisting of piecewise linear functions in space. The finite element method for (3.1) reads: find uh∈Vh such that, for all v∈Vh,

where Jλ′⁢(uh) is the Fréchet derivative of functional (3.1) derived in [4].

Let us briefly present the balancing principle for choosing the regularization parameter in (3.1). According to [28, 18], for the Tikhonov functional

the value function G⁢(λ) is defined as

If there exists derivative Gλ′⁢(λ) at λ>0, then from (3.3) and (3.4) follows that

The balancing principle for finding the optimal regularization parameter finds λ>0 such that the following condition is satisfied: φ¯⁢(λ)=C⁢λ⁢ψ¯⁢(λ), where 𝐶 is determined using a priori information about the solution [18]. We will consider the case when C=1. In our computations, we have used the fixed point algorithm, summarized in Algorithm 1, whose convergence was studied in [18].

This method is based on the next theorem which follows from [23, Theorem 5.1].

According to the above theorem, the finite element mesh Kl:=(Kh)l, with l=0,1,…,M denoting number of mesh refinements, is locally refined in the neighborhood of those points of slice Ωz where the computed residual vector R⁢(uh):=|h⁢uh| attains its maximal values. To get the function ml:=(mh)l on the new refined mesh Kl, l=1,…,M, we use interpolation of ml-1 from the mesh Kl-1 to the new refined mesh Kl using the procedure described in [3]. For ease of access and more clarity, the whole process of AFEM is tabulated in Algorithm 2.

4 Results and discussion

Reconstruction results using a regularized least squares method for the scenario described in Section 2.1 are shown in Figure 2. From the plots, it can be seen that the shape and location of the heated target have been retrieved successfully and that the contrast grows over time as the scenario keeps deviating more and more from the baseline. As explained earlier, this is a qualitative indicator of changes in constitutive parameters of the target, which indirectly refers back to its temperature changes over time. However, the changes in the dielectric property of the target, as reported in Table 1, are somewhat arbitrary in this study and more importantly are uncorrelated to any thermal changes in the target because they are just simply taken from [16] and then rescaled. Therefore, in a future study, a proper correlation between these properties is necessary so as to map the changes in |Δ⁢O⁢(r)| back to temperature.

Figure 2

Time series of |Δ⁢O⁢(r)| every 2 minutes from t=2 up to 10 minutes. First and second rows are 2D cuts in the transverse and lateral planes, respectively. The third row is a 3D visualization of the same scenario as the contrast changes over time. The border of the target is also shown with a white dashed line in the plots.

(a)

(b)

(c)

Figure 3

Left figures: convergence of fixed point algorithm for computation of the optimal regularization parameter 𝜆. Here, 𝑙 is the number of mesh refinement, nno is the number of nodes and nel is the number of elements in the mesh Kl. Right figures: convergence of AFEM on adaptive locally refined meshes. All results are presented for reconstructions shown in Figure 4.

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

(i)

(j)

To test feasibility of AFEM, we set-up the dimensionless computational domain

and decompose it into two-dimensional slices Ωz for any desired value of z∈[0,25]. The optimal value of the regularization parameter 𝜆 was computed by the fixed point algorithm for an initial guess for the regularization parameter λ0=0.1. Convergence of fixed point algorithm is presented in the left column of Figure 3. These figures show that the optimal 𝜆 was found on the interval λ∈[0.46,0.48] on all refined meshes. Then the AFEM computations were performed in the software package WavES [34] with β=0.5 in (3.11) which allows to refine meshes locally at places where the inclusion was located; see Figures 4, 5.

In the transverse plane, Figure 4 compares the original reconstructions on the unrefined initial FEM mesh with that of AFEM on locally adapted meshes at each time step. One can see from these plots how the density of the mesh changes around the border of the target and how this change becomes more and more pronounced as the contrast grows bigger and bigger due to thermal changes of the target as time goes by. Figure 5 does the same but in the lateral plane. All in all, one can argue that these figures clearly show improvement in the shape detection of the target. Contrast, which means better resolvability of the target to be imaged from its background, has also been improved significantly. Finally, convergence of AFEM is presented in the right column of Figure3. All convergence results are given for reconstructions shown in Figures 4. Analyzing results of Figures 3–5, we can conclude that reconstructions obtained by AFEM improve the initial reconstructions of the target.

We note that an adaptive finite element method in three dimensions with mesh refinement in the whole three-dimensional domain Ω is an ongoing work, and it will be presented in a forthcoming publication.

Figure 4

Reconstructions in the transverse plane. The first two figures on the left are the reconstruction obtained by formula (2.7) and the same reconstruction on the coarse finite element mesh. The next two figures are the reconstruction obtained by AFEM using Algorithm 2 with β=0.5 and the same solution on the locally refined mesh.

(a)

t=2 min

(b)

t=4 min

(c)

t=6 min

(d)

t=8 min

(e)

t=10 min

Figure 5

Reconstructions in the lateral plane. The first two figures on the left are the reconstruction obtained by formula (2.7) and the same reconstruction on the coarse finite element mesh. The next two figures are the reconstruction obtained by AFEM using Algorithm 2 with β=0.5 and the same solution on the locally refined mesh.

(a)

t=2 min

(b)

t=4 min

(c)

t=6 min

(d)

t=8 min

(e)

t=10 min

5 Conclusions

A feasibility study on the numerical characteristic integration of an FDTD-based forward solver with an inversion scheme has been presented. The standard scheme has then been supplemented by AFEM to gain higher quality reconstructions during the course of the cooling process of a heated object. More specifically, regularized least squares, combined with an adaptive finite element method, is used for image restoration. A balancing principle is also applied for optimal choice of the regularization parameter. Numerical tests show that an adaptive algorithm improves image quality compared with that of the standard techniques.

Current and future work considers development of fully 3D adaptive refinement algorithm and its application on real experimental scenarios. One possibility for a future study is to have a hyperthermia treatment planning for the scenario in which the electromagnetic power deposition will be first beamformed to focus on the target and then passed on to a thermal simulator to calculate the corresponding temperature distribution. Finally, to validate the performance of the thermal imaging presented here, dielectric properties of ROI can be updated based on their temporal temperature distributions of the previous step and then fed into the reconstruction algorithm at each time step accordingly.