Abstract

Microwave thermal ablation is a cancer treatment that exploits local heating caused by a microwave electromagnetic field to induce coagulative necrosis of tumor cells. Recently, such a technique has significantly progressed in the clinical practice. However, its effectiveness would dramatically improve if paired with a noninvasive system for the real-time monitoring of the evolving dimension and shape of the thermally ablated area. In this respect, microwave imaging can be a potential candidate to monitor the overall treatment evolution in a noninvasive way, as it takes direct advantage from the dependence of the electromagnetic properties of biological tissues from temperature. This paper explores such a possibility by presenting a proof of concept validation based on accurate simulated imaging experiments, run with respect to a scenario that mimics an ex vivo experimental setup. In particular, two model-based inversion algorithms are exploited to tackle the imaging task. These methods provide independent results in real-time and their integration improves the quality of the overall tracking of the variations occurring in the target and surrounding regions.

1. Introduction

Thermal ablation is an increasingly exploited cancer treatment in which high temperatures are induced into the target zone, represented by the tumor plus a 5–10 mm safety region of surrounding normal tissue, in order to create an irreversible thermal damage. Temperatures as high as 55–60°C are needed to induce the so-called coagulative necrosis of cells [13]. The main advantages of this cancer treatment with respect to surgery interventions are related to the reduced invasiveness and costs of the treatment. Moreover, thermal ablation has minimal side effects and can be performed also in nonsurgical patients.

Depending on the heating source, it is possible to distinguish different ablation modalities. Some examples are focused ultrasound [4] and electromagnetic techniques [5]. These latter include two different modalities according to the operating frequency, which correspond to substantial different mechanisms of interaction between electromagnetic (EM) field and target tissues.

In radio frequency (RF) ablation, which is currently the most spread ablation modality, the ablation occurs in the near proximity (a few millimeters) of an electrode with an applied alternating current (in the 350–500 kHz frequency range), but a larger area is actually ablated thanks to thermal conduction, provided that involved tissues exhibit sufficient thermal conductivity.

Microwave (MW) ablation exploits as interstitial applicator an antenna typically radiating in the industrial, scientific, and medical (ISM) band (MHz or GHz). MW ablation occurs in a volume that surrounds the applicator, and not just in its proximity. Such a circumstance represents the main advantage of MW ablation with respect to RF one: the capability of treating larger volumes in a shorter time, which in turn leads to a reduced sensitivity to blood perfusion heat sink [3].

One of the main shortcomings of MW thermal ablation (MTA) is the lack of a reliable temperature monitoring system. Typically, the temperature in the region of interest is recorded with local thermometric probes. However, these have the drawback of being invasive and capable of providing only local information, so that the adoption of noninvasive monitoring techniques would be advisable. In particular, magnetic resonance [6] and ultrasound [4] imaging have been proposed in the literature (and tested in the clinics) to tackle this issue. However, besides electromagnetic compatibility issues, magnetic resonance involves high cost equipment with obvious drawbacks in terms of economic sustainability, whereas the actual effectiveness of ultrasounds when temperature increases is still an open issue, since transducers are blinded by a hyperechogenic cloud caused by water vaporization as temperature increases up to about 100°C.

The above reported limitations of existing imaging modalities motivate the interest in alternative techniques, and a possible candidate in this respect is microwave tomography (MWT), which images the electromagnetic properties of the tissue using low-cost and portable equipment. In fact, recent works have put into evidence tangible changes of both dielectric permittivity and electric conductivity of the tissue linked to the very high temperature increases (up to 120°C) achieved during the microwave thermal ablation procedures [79]. According to these measures, a decrease of both relative permittivity and electric conductivity occurs with the temperature. In particular, starting from room temperature both quantities started decreasing with an almost linear behavior. Then, in the temperature range between 60°C and 80°C a significant decrease was measured, followed by (above 90°C) a rapid drop. These changes proved to be nonreversible and can be associated at temperatures close to 60°C–80°C with proteins denaturation [10]. At temperatures close to 90°C–100°C the rapid drop of both dielectric properties can be linked to loss of water due to the vaporization phenomenon and to diffusion of the generated water vapor towards the lower-temperature surrounding tissues [11].

In principle, these changes can be imaged using MWT by probing the treated region with an array of antennas, external to the patient body. This array would record the evolution of the back-scattered field during the treatment and the proper processing of such data would allow recovering the spatial map of the complex permittivity.

In [12], a first feasibility study, concerned with an ideal configuration (i.e., a very simplified and invariant 2D scenario and ideal probes) was presented. Based on the positive outcomes of that study, this paper further explores the idea of monitoring MTA using MWT. In particular, a processing technique to form the images of the evolving scenario, based on the joint exploitation of the results achieved by two independent MWT methods, is herein developed. Both MWT methods process the same data and allow real-time results, which are necessary for the application of interest. The image fusion is based on a very simple merging of the partial results, but allows taking full advantage of the cross-validation between the two adopted imaging methods.

The overall strategy is tested against accurate, full-wave, simulated ablation experiments mimicking the typical conditions faced in ex vivo ablation studies [9, 13]. The reason for considering an ex vivo scenario is twofold. First, the success of the numerical study presented here allows us to confidently tackle a subsequent ex vivo experimental validation, which is the next foreseen step to assess the concept. Second, monitoring ex vivo experiments is itself an important issue, as this kind of studies is fundamental to properly test MTA devices and determines reliable clinical protocols. Note that some preliminary results involving one of the two methods and a similar scenario have been presented in [14].

The paper is organized as follows. Section 2 describes the adopted materials and methods, providing an overview of the basis of MWT and the details of the two adopted imaging strategies. In Section 3, the numerical test-bed adopted for the simulated assessment is described first and then the imaging results are presented. Discussion and conclusions follow.

2. Materials and Methods

2.1. Microwave Tomography: Basic Principles

MWT obtains images of an unknown target by appraising its influence on a known probing wave. In particular, the contrast in electric properties (dielectric permittivity and electric conductivity) is imaged.

The physical phenomenon on which MWT relies on is electromagnetic scattering. When an electromagnetic field interacts with an object, it induces electric currents in it, known as induced currents or contrast sources. These currents, in turn, radiate an additional electromagnetic field, known as scattered field. Notably, the scattered field depends on the morphological and electromagnetic features of the object.

Let us denote with the volume of interest and let us suppose that, at the frequency , this volume is probed by means of antennas located on a surface external to . The scattered field arising from the interaction between the probing fields and is gathered by another set of antennas located on a surface (external to ). Of course, the same antennas can be used as transmitters and receivers; in this case, coincides with and is at most equal to .

In these general conditions, the electromagnetic scattering is modeled through the following integral equation [15]:where , . In (1), denotes the scattered field, given by the difference between the total field (i.e., the solution of the Helmholtz equation in the considered scenario) and the reference/background field (i.e., the solution of the Helmholtz equation when no target/object is present in ). is dyadic Green’s function of the reference scenario (when no target is present), that is, the electric field measured at due to an elementary source located in . denotes the contrast source induced in when is probed by the source in , and it is given by the product of the total field times, the contrast function , which accounts for the electric difference existing between the object and the background scenario .

To retrieve the properties of the region encoded in the contrast function from the knowledge of the scattered field and the background field, an inverse electromagnetic scattering problem has to be solved [16]. It is worth noting that, in actual MWT experiments, the receiving antennas measure the total field, so that the scattered field used as input datum of the inverse scattering problem has to be properly extracted from such a measured quantity.

2.2. MWT for Temperature Monitoring: Concept and Formulation

The idea of adopting MWT for MTA monitoring relies on the above-mentioned experimental studies reporting the changes of dielectric properties of tissues during the ablation process [79]. In particular, these studies evidenced an irreversible temperature-dependent change in the dielectric parameters that can be characterized in terms of complex permittivity:where is the vacuum permittivity ( Fm−1), , is the relative permittivity, and (Sm−1) is the electric conductivity.

The conceptual scheme of the idea underpinning the present work is represented in Figure 1. A (limited) number of antennas are positioned above the tissue sample, while the ablation treatment is taking place. The antennas record the evolution of the scattered field, which is then processed to obtain the spatial map of the complex permittivity variation as the tissue is ablated.

To this end, since the relevant information concerns the variation with respect to the initial situation, a differential approach appears as a suitable choice. In such an approach, the electric contrast is given by the variation which occurs between the “time-zero” at which the treatment starts and the time instant of interest, that is, , with and being the complex permittivity of the scenario at time and , respectively. Accordingly, the scattered field required to image such a contrast is given by the difference between the (total) fields recorded by the antennas in the two instants, while the reference field is the (total) field in the region to be treated before the treatment starts and just after the insertion of the MTA applicator (), that is, .

According to the above, (1) can be modified to describe this differential scattering framework aswhere is the differential scattered field at , given by , and is the differential contrast source given by

The imaging problem underlying MWT monitored MTA is hence cast as retrieving, at each time of interest, the differential contrast from the measured , given the reference scenario (i.e., and ).

As can be seen, from the fact that depends on that implicitly depends on the contrast , this is a nonlinear problem. Moreover, due to the properties of , it belongs to the class of ill-posed problems.

In the next section, we discuss how to approach these difficulties and obtain reliable images of from the measured data.

2.3. Imaging Strategies

The following subsections describe the approaches adopted to cope with the imaging task. It is worth noting that the two strategies are based on completely different concepts, so that their results can be cross-validated and merged to increase the confidence in the imaging results. Moreover, it is also important to remark that they allow to achieve results in real-time, although these do not include quantitative information on the properties’ variations. While this is obviously a limitation, qualitative information can meet the first monitoring goal for ablation treatments that is to appraise the extent of the region which is being actually ablated and the boundary between treated and untreated tissue.

It is worth noting that the exploited imaging strategies are well known, but their fusion has never been exploited.

2.3.1. Truncated Singular Value Decomposition (TSVD) Inversion

The first approach relies on the assumption that the variation of the contrast function during the observed time interval is small enough to leave almost unperturbed the background field. Accordingly, the relationship between the data and the unknown can be linearized by replacing with in (4):This corresponds to the so-called distorted Born wave approximation [16]. To solve such a kind of linear and ill-posed inverse problem in a regularized fashion, simple and well-assessed tools can be exploited. In particular, the truncated singular value decomposition (TSVD) scheme [16] is adopted herein.

In the TSVD, the unknown contrast is retrieved by means of the explicit inversion formula:wherein is now the vector containing the samples of the reconstructed differential contrast function at the considered time, with being the number of the cells into which the investigation domain has been discretized, is a vector (with and denoting the number of transmitting and receiving probes, resp.) containing the values of the main component of the differential field at the considered time, and denotes the singular value decomposition (SVD) of the matrix representing discrete counterpart of the integral operator in (5), which is computed with respect to the unperturbed scenario. Finally, is the index at which the summation is truncated, which acts as a regularization parameter, being set in such a way to meet accuracy and stability of the reconstruction [16].

The main advantage of the TSVD approach is that it enables real-time imaging results. As a matter of fact, the main computational tasks are building and computing its SVD, which can be both performed offline, since, according to our formulation, is independent of the considered time instant. On the other hand, such a linear inversion scheme cannot provide quantitative information on the unknown contrast, since the underlying distorted Born approximation does not hold true in the considered scenario. As a matter of fact, as shown by the experimental results [9], the magnitude of variations of the permittivity occurring during ablation and their spatial extent violates the assumed “weak” scattering conditions. As a consequence, the image obtained using the inversion formula (6) can only provide an estimation of the location, shape, and extent of the variations in the electric properties distributions occurring in the region of interest due to the ablation process.

In the following numerical analysis, only the main component of (i.e., the one aligned along the orientation of the receiving antennas) will be considered, but the vectorial nature of the problem is preserved by considering the corresponding “row” of dyadic Green’s function and all the three components of the electric field generated in the imaging domain by the primary sources.

2.3.2. Linear Sampling Method (LSM)

The second approach considered is the linear sampling method (LSM) [17]. The LSM belongs to the class of qualitative imaging methods. The main feature of these methods is that they allow retrieving morphological information on the imaged scenario (in the case at hand: the location, shape, and extent of the variations in the electric properties due to the ablation process) through the solution of a linear inverse problem. Notably, this linear problem does not require approximations and the solution of the inverse problem can be again performed in real-time.

The LSM image is given by the plot of a suitable “indicator” function, which attains low values in points belonging to the target and large values elsewhere. In particular, the modified LSM formulation for differential imaging developed in [18] is adopted herein. In such a formulation, the value of the indicator function in a generic spatial position of the imaged region is given bywherein and , respectively, represent the th singular value and th right singular vector of the differential data matrix (where and denote the number of transmitting and receiving probes, resp.) that contains the samples of the measured component of the differential scattered field at the considered time and is a vector that contains the values of the electric field radiated by an elementary source located in and oriented in the spatial direction , as measured at the receiving probes. By exploiting reciprocity, the entries of this vector are simply given by the projection of the reference field along , that is, .

In (7), and is the regularizing coefficient that is set as [18]

It is worth noting that the orientation of the test source is a degree of freedom in LSM. Therefore, (7) represents a family of indicator functions, having as parameter. For simplicity, in this feasibility study only the orientation parallel to the main component of the measured electric field has been considered.

2.3.3. Image Fusion

As shown above, the two considered imaging methods are complementary, as they exploit the same data and quantities, but in a different fashion. Moreover, they both allow achieving real-time imaging results, which are related to the morphology of the variations occurring in the ablated tissue. This suggests that a further improvement in the quality and reliability of the imaging results can be achieved by fusing the two obtained results.

The fusion strategy herein adopted consists in summing, in a pixel-wise fashion, the TSVD, and LSM images after normalization, and clipping the obtained result to 1. This leads to images in which the reported quantity continuously varies between zero and one, with zero meaning that no variation is occurring in the dielectric properties. In particular, in the ablation treatment, larger variations occur close to the applicator, while far from the applicator the map is expected to be zero.

To facilitate the comparison between the results obtained in the different conditions, the “fusion” maps have been binarized. This also allows to directly appraise (both visually and quantitatively) the extent of the estimated ablated area. To this end, we have set a hard threshold at −3 dB. This choice is suggested by the fact the ideal result of the monitoring algorithm should be a binary image, with ones representing the pixels where a change in the dielectric properties occurred and zero representing the pixels where no changes were detected. However, MWT techniques provide a low-pass filtered image of the actual spatial profile, which, in the particular case here at hand, becomes a sinc-like decay from the area of maximum variation (close to the antenna applicator) to the area of no variation detected (outside the thermally ablated area). Accordingly, setting the threshold at −3 dB corresponds to assume that at least of the whole “energy” is confined in the selected area.

3. Simulation of an MWT Monitored Ex Vivo Ablation Experiment

In this section, the capabilities of the MWT approach described above to monitor an ablation process are assessed simulating a scenario that mimics ex vivo experiments.

3.1. Numerical Test-Bed

The considered scenario for the MWT algorithm is shown in Figure 2. The MTA applicator plays a passive role so that it was schematized as a thin metallic cylinder (radius mm). The MWT data are collected by means of dipole antennas (of mm in length, hence very small with respect to the probing wavelength, which is about cm). The antennas are evenly spaced on a -plane parallel to the major axis of the MTA antenna applicator and are assumed to be in direct contact with liver tissue. The distance between the MWT antennas array plane and the MTA antenna applicator is approximately cm.

To model the variations of electromagnetic properties of liver tissues during the treatment, results from previous experimental studies were exploited [7, 9, 13]. These results correlate the temperature changes occurring during ablation with the corresponding changes in the dielectric parameters. Since these experiments were carried out at GHz, the same frequency has been used in this work to simulate the electromagnetic scattering phenomenon. However, it is worth noting that the operating frequency adopted by the imaging system is actually a degree of freedom in the proposed technique, which can be set according to practical needs (e.g., minimizing the inference between the MWT array and MTA applicator when this latter is on, increasing wave penetration and so on). On the other hand, for the purposes of the current assessment, it is more important to model the scenario using reliable experimental values for the electric properties of tissue, so that the frequency of GHz has been considered. Moreover, the ablated region has been assumed to have a rotational symmetry around the axis of the applicator, as evidenced in [7].

On the basis of the experimental results, the 3D ablated volume of interest has been segmented in four different tissue types, in order to simulate the initial conditions,  s, when only untreated liver tissue is present, and the scenario after  s of treatment (named Ablated Livers 1, 2, and 3, resp.), when liver tissue in an increasingly large region surrounding the applicator is progressively ablated. The electromagnetic properties of the liver tissue in the 4 considered conditions are reported in Table 1, while the maps of the real part of the permittivity () of liver tissue after ablation at the three considered times are shown in Figure 3, as slices of the 3D voxel model used in the simulation procedure.

With respect to these conditions, full-wave numerical simulations using CST Microwave STUDIO® (CST Computer Simulation technology AG) have been performed to compute the quantity needed by the MWT approaches. For simplicity, in this feasibility study, the background is set as homogeneous healthy liver, neglecting the presence of the interface between tissue sample and sample holder/air. While this is obviously a simplification, it is important to note that this circumstance can be practically achieved by surrounding the sample holder with absorbing panels, which minimize the reflections from boundaries. Note that only the main component of the total measured field for the MWT (-component in Figure 2) has been considered.

3.2. Imaging Results

In this section, the results obtained by means of the two considered inversion strategies are presented. In particular, for both methods, the tomographic reconstructions obtained in a slice of the 3D volume parallel to the antenna array plane and located just above the applicator are shown.

First, the case of noiseless data has been considered to appraise the feasibility of the proposed approach. In this respect, it is worth noting that no “perfect” estimation of the ground truth is expected, even in these ideal conditions. As a matter of fact, the measurement configuration exploits only a very reduced number of antennas and is aspect limited. Moreover the electromagnetic scattering model used to cast the TSVD inversion (see (5)) assumes that the total field can be approximated with the background field (the field at ). As said, this is not the case, owing to the strong variations of the electromagnetic properties occurring among the different considered time instants. Such a model error prevents the perfect reconstruction of the unknown profile and actually avoids that the so-called inverse crime [19] is committed processing noiseless, ideal, data. To tame such an error, the inversion formula (6) does not retain all terms of the summation but only considers those corresponding to singular values not smaller than dB with respect to the leading one (i.e., ). Also for the LSM, the image is achieved using the regularized solution given in (7), so that similar considerations apply.

For the TSVD technique, the map of the normalized amplitude of the retrieved contrast at the considered time instant is reported using a threshold (see (6)), while for the LSM for each instant the map of the normalized indicator is given:The same colormap is adopted to facilitate a visual comparison of the obtained results. For each time instant, the fusion of the TSVD and LSM images is shown as well, with superimposed red contour line representing the boundary between untreated and initially ablated liver tissue (Ablated Liver 1 in Table 1).

The results are summarized in Figure 4 and show that both methods perform quite well in tracking the occurring changes. In particular, Figures 4(a), 4(d), and 4(g) provide the TSVD reconstruction, while Figures 4(b), 4(e), and 4(h) show the LSM indicator. The results of the adopted image fusion are shown in Figures 4(c), 4(f), and 4(i) and confirm that while being simple, the adopted fusion strategy provides satisfactory estimations of the extent of the ablated region. It is worth noting that a slight underestimation of the ablated area is observed at  s, which is even worse at  s. This is in agreement with the increase of the model error as the ablated area enlarges.

With the aim of considering an even simpler measurement configuration, easier to be deployed in practice, we have performed the tomographic reconstructions using a reduced number of probes. In particular, a measurement configuration made by only antennas (the antennas labeled as 1–5 in Figure 2) has been considered. Again, noiseless data are processed and in these cases the threshold of the TSVD in (6) is set to .

As can be seen, from the results shown in Figure 5, the reduction of the number of probing/sensing antennas further enhances the low-pass nature of the tomographic reconstruction process, leading to the presence of artifacts on the border of the imaging domain. Notably, these artifacts are located at a distance equal to from the main lobe, which suggests that they are possibly spurious effects (i.e., a sort of “grating” lobes) rather than actual features. In any case, also in this case the advantage of fusing the imaging results obtained from the two methods is clear.

Finally, the effect of the unavoidable presence of noise on the data has been appraised by adding Gaussian noise on the total “measured” field; in particular two signal-to-noise ratio (SNR) levels have been considered, dB and dB. It is worth noting that, in the differential framework we are considering, these SNR values do not correspond to low noise levels. As a matter of fact, the useful signal is orders of magnitude below the total (measured) field, so that it will be heavily corrupted by noise. Such an issue is well known in medical applications of MWT, wherein systems are indeed designed in order to meet similar requirements in terms of SNR and dynamic range [20, 21].

The results of such an analysis are reported in Figures 6 and 7, respectively, and confirm the importance of exploiting both methods, in order to take advantage of the image fusion. As a matter of fact, while the single approaches experience a significant degradation, the image fusion results appear to be more reliable.

To quantitatively appraise the achieved results, two error metrics have been defined. The first metric is the coverage factor (CF), defined as the percentage of pixels belonging to the actual ablated area wherein the fusion image assumes that values equal 1. The second metric is the complement to the CF (cCF), instead defined as the percentage of pixels external to the actual ablated area wherein the fusion image is . Note that the ideal figures for these metrics are for the CF and for the cCF.

The analysis of the CF and cCF for the fusion results are reported in Figure 8, where the metrics for both the considered array configurations, all time instants, and levels of noise are reported as histograms. Notably, the analysis shows the robustness against the considered level of noise, with a slight degradation when considering the reduced antennas configuration, due to the increase of cCF. It is worth nothing that both metrics have to be taken into account to appraise the quality of the image. As a matter of fact, a CF equal to , if paired to a large cCF, indicates a significant overestimation of the ablated area. This happens, for example, with the reduced measurement configuration and dB.

The quantitative metrics CF and cCF have been computed also for TSVD and LSM separately, but the achieved results were omitted for brevity. In summary they confirm that the fusion image performs better than the single methods alone.

While these are obviously preliminary results, they support the viability of the proposed concept in real conditions, wherein this number of antennas can be actually deployed and where the measurement accuracy level is expected to be achievable, considering, for instance, that typical achieved by microwave imaging devices for breast cancer [20, 21].

4. Conclusions and Discussion

In this paper, a feasibility study concerning the adoption of MWT as a real-time noninvasive tool to monitor MTA treatment was presented. To this end, experimental results on the electromagnetic properties variations obtained in an ex vivo liver sample during the ablation were exploited [79].

A differential tomographic approach has been developed to tackle the imaging problem at hand. In particular, two different (and independent) imaging methods have been exploited, both providing reliable results. This has suggested postprocessing the achieved tomographic images by combining them into a single fusion image, with the aim of further improving the quality and reliability of the final result. The presented numerical assessment confirms the possibility of tracking the complex permittivity variations of the liver tissue while ablated, by means of MWT, even using a realistic (low) number of probes.

As a preliminary study, oversimplified scenarios were considered in the performed numerical simulations. More realistic examples as well as some experimental proof of concept experiments are matter of future work. On the other hand, it is worth noting here that the presence of in-homogeneities, real antennas patterns and so on should be easily compensated in a differential scheme, like the one here considered, since such complexities will be fixed and still during the MTA treatments.

Based on these initial encouraging results, the activity is planned to progress along several directions. First, the experimental validation of the technique concerning with ex vivo scenarios will be pursued to complement and confirm the results herein achieved. In addition, further numerical studies concerned will be carried out. From one hand, scenarios closer to the clinical ones will be considered to account for the presence of other tissues surrounding the liver and assess the monitoring performances of MWT in these conditions. On the other hand, multifrequency approaches will be also assessed, with the aim of increasing the amount of available data, despite the low number of performed measurements.

Then, more complex imaging procedures, which take into account the heterogeneity of the scenario while the ablation treatment progresses, must be developed in order to further improve the accuracy of the estimation of the actual dimension of the ablated areas, even for larger ones. One possible solution is represented by a distorted approach, in which the differential data are taken with respect to a shorter time interval to reduce the magnitude of contrast variations, thus better matching the assumptions underlying the Born approximation. This case would be possible by processing the differential scattered field obtained by subtracting data collected at intermediate steps of the ablation treatment. Of course, this would require, as an additional effort, to update the background field (and the Green’s function) at each step, by solving the relevant forward scattering problem on the basis of the intermediate results.

Finally, the assessment of quantitative MWT approaches could also be considered as a way to directly have information on the temperature distribution (linked to the permittivity and conductivity map) inside the region under treatment. It is worth noting that in a quantitative analysis one cannot neglect the uncertainty between the temperature values and the corresponding value of the dielectric properties, not taken into account in the present work where only qualitative approaches were considered.

Conflicts of Interest

The authors declare that they have no conflicts of interest.

Acknowledgments

This work has been developed in the framework of COST Actions TD1301 MiMed and BM1309 EMF-MED. The authors are glad to acknowledge Professor O. M. Bucci for stimulating them to undertake this research activity.