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

: Integration of an adaptive finite element method (AFEM) with a conventional least squares method has been presented. As a 3D full-wave forward solver, CST Microwave Studio has been used to model and extract both electric field distribution in the region of interest (ROI) and S-parameters of a circular array consisting of 16 monopole antennas. The data has then been fed into a differential inversion scheme to get a qualitative indicator of how the temperature distribution evolves over a course of the cooling process of a heated object. Different regularization techniques within the Tikhonov framework are also discussed, and a balancing principle for optimal choice of the regularization parameter was used to improve the image reconstruction quality of every 2D slice of the final image. Targets are successfully imaged via proposed numerical methods.


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 [10,14,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 [8,12,16,25]. 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][29][30]. 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.

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.  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.

Differential image reconstruction
Throughout this section, e −iωt is assumed to be the time convention and suppressed. Let us suppose we have a bi-static pair (i, j) of antennas located on the scan line Γ, i.e. r i , r j ∈ Γ, and for the sake of simplicity, we use the notation E s (r j , r i ) = E s ji . Then, using the Lorentz reciprocity theorem and under Born approximation, the scattered electric field between the pair at angular frequency of ω can be written as where Ω is the imaging domain, I(ω) is the excitation current of the transmitter, k b is the lossless background wavenumber, and G and O = ∆ϵ r + i ωϵ b ∆σ are the dyadic Green function and the object function of the scenario, respectively. Equation (2.1) is the standard Fredholm integral equation of the first kind in which the dyadic Green function relates the fields in the object domain (vector) to the scattered fields (vector). However, in heterogeneous complex scenarios, there is no closed-form analytic expression for G. Moreover, routine measurements performed with vector network analyzers (VNAs) gather S-parameters of the antennas, which are scalar instead of the vector quantities of the scattered field at the location of the probes. Suggested by [15], an efficient way to address both issues is to integrate a full-wave commercial software, like CST in this study, into the abovementioned formula which obviates the need for a closed-form Green function. Finally, modifications, brought about by considering the scattered fields replacement with their corresponding S-parameters as well as the input power and characteristic impedance of the ports, change equation (2.1) into the following: where C = −k 2 b /(4iωμ) and E CST inc,i is the exported E-field from CST under irradiation of the i-th antenna. Using equation (2.2) and considering the object function O b as the baseline, we can then proceed to relate according to the following formula [16]: After discretizing ROI and considering A as the linear forward operator, equation (2.3) can be written as Am = d with m = ∆O and d = ∆S. As discussed earlier, the problem is ill-posed and hence requires a regularization scheme. To this end, the Tikhonov regularization functional is defined as follows [28,29]: where λ is the regularization parameter. The optimal value of m is found via the solution of the following optimization problem which is basically the Fréchet derivative of functional (2.4) [4]: Looking closely at equation (2.5), one can recognize that system (2.5) is the regularized system of normal equations with the matrix A corresponding to the operator A, which has the following solution if normal method is to be used: Although the method of normal equations is the fastest method for solving the linear least-squares problems (LLSP), it is not as accurate as QR or SVD decompositions. For instance, the method of normal equations is applied to the solution of LLSP when the condition number of the matrix A is not large [5]. Alternatively, using SVD of the matrix A = UΣV T and substituting it into (2.6), one can obtain a much more robust solution for m as follows: Finally, it is not always trivial to determine the optimal value of λ, and it differs vastly from scenario to scenario. This has been a matter of debate in the relevant literature and has been widely discussed elsewhere [1,2,18,20,26]. Here, we just try to briefly touch upon the issue. The goal of regularization is to construct and m * is the "ideal" exact solution which cannot be found in real experiments but we assume that it exists; see details in [1,2,18,20]. In [1,5], it is shown that condition (2.8) is satisfied when the following iterative update of the regularization parameters λ k is considered: where k is the iteration number and λ 0 is the initial guess. So, starting from λ 0 = 1 and using (2.9), optimum λ is found to be 0.9 for this scenario.

Adaptive finite element method with balancing principle
For improvement of the solution m 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], m is the solution obtained via (2.7) and ‖ ⋅ ‖ 2 is the L 2 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 u which minimizes functional (3.1), we look for a stationary point of this functional with respect to u such that J λ (u)(b) = 0, where J λ (u) is the Fréchet derivative of functional (3.1) acting on b; 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 z. Then every slice is discretized by the mesh K h which consists of non-overlapping triangles K satisfying the mesh regularity assumption [19] with the mesh function h which is defined as a piecewise-constant function such that h = h K , where h K is the diameter of K which we define as a longest side of K.
Further, we introduce the finite element space V h consisting of piecewise linear functions in space. The finite element method for (3.1) reads: find u h ∈ V h such that, for all v ∈ V h , where J λ (u h ) 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 [18,28], for the Tikhonov functional 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 C 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].
(1) Start with the initial approximations λ 0 , and compute the sequence of λ k in the following steps.
(3) Update the regularization parameter λ := λ k+1 as such that the following a posteriori error estimate for the regularized solution u λ and for all u h ∈ V h holds:

Here, C I is the interpolation constant and R(u h ) is the residual R(u h
Proof. Let u h be the minimizer of Tikhonov functional (3.1). Then functional (3.1) is strongly convex on the space H 1 with the strong convexity constant λ, and thus where J λ (u h ), J λ (u λ ) are the Fréchet derivatives of functional (3.1). Since u λ is the minimizer, where P h is the orthogonal projection of u λ , together with the Galerkin orthogonality principle The right-hand side of (3.8) is estimated via (3.5) as

Using this equation in (3.8), we obtain ‖u
. By the property of orthogonal projection, ‖P h u λ − u λ ‖ H 1 (Ω) ≤ C I ‖hu λ ‖ H 2 (Ω) , with the interpolation constant C I , we get the a posteriori error estimate The term ‖hu λ ‖ H 2 (Ω) can be estimated for all u h ∈ V h similarly to [23, Theorem 5.1]: Using estimate (3.10) in the right-hand side of (3.9), we get a posteriori estimate (3.6).
According to the above theorem, the finite element mesh K l := (K h ) 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(u h ) := |hu h | attains its maximal values. To get the function m l := (m h ) l on the new refined mesh K l , l = 1, . . . , M, we use interpolation of m l−1 from the mesh K l−1 to the new refined mesh K l using the procedure described in [3]. For ease of access and more clarity, the whole process of AFEM is tabulated in Algorithm 2.
(1) Generate the initial mesh K 0 in Ω z . Obtain m 0 via (2.7). Compute the sequence of reconstructions u l := (u h ) l , l ≥ 0, on the refined meshes, as follows. (2) Compute the finite element solution u l on K l using the finite element method (3.2).

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.  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  Here, l is the number of mesh refinement, nno is the number of nodes and nel is the number of elements in the mesh K l . Right figures: convergence of AFEM on adaptive locally refined meshes. All results are presented for reconstructions shown in Figure 4.
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.

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.