M ULTIPHYSICS D ESIGN O PTIMIZATION OF M ICROWAVE A BLATION A NTENNAS by S HASHWAT S HARMA A thesis submitted in conformity with the requirements for the degree of Master of Applied Science Graduate Department of Electrical and Computer Engineering University of Toronto c Copyright 2016 by Shashwat Sharma ProQuest Number: 10194483 All rights reserved INFORMATION TO ALL USERS The quality of this reproduction is dependent upon the quality of the copy submitted. In the unlikely event that the author did not send a complete manuscript and there are missing pages, these will be noted. Also, if material had to be removed, a note will indicate the deletion. ProQuest 10194483 Published by ProQuest LLC (2017 ). Copyright of the Dissertation is held by the Author. All rights reserved. This work is protected against unauthorized copying under Title 17, United States Code Microform Edition © ProQuest LLC. ProQuest LLC. 789 East Eisenhower Parkway P.O. Box 1346 Ann Arbor, MI 48106 - 1346 Abstract Multiphysics Design Optimization of Microwave Ablation Antennas Shashwat Sharma Master of Applied Science Graduate Department of Electrical and Computer Engineering University of Toronto 2016 Microwave ablation is based on localized heating of biological tissues, enabled by an electromagnetic field. Antennas for ablation are commonly designed in a forward approach to generate a particular Specific Absorption Rate profile within the target. However, little attention has been dedicated to designing antennas inversely, to allow controllable synthesis of temperature profiles customized for the application. Also, most existing designs do not account for thermal inhomogeneities in tissue. An inverse, multiphysics methodology for microwave ablation antenna design is presented herein, that involves optimizing the antenna’s current distribution to synthesize a desired temperature profile, while accounting for heat diffusion. The results indicate and quantify the clear advantages of this multiphysics approach. This methodology is successfully applied towards designing an easily configurable printed dipole microwave ablation antenna that addresses several limitations of existing designs. This design provides the functionality of a phased array, and allows controlled synthesis of temperature profiles. ii To the shoulders on which I have stood. iii Acknowledgements I would like to thank several people who have played an immense role in my work and personal development over the last two years. Firstly, I want to express my gratitude to Professor Costas D. Sarris, my supervisor, for entrusting me with this project, for his patience and constant guidance and insight, and for supporting me in a way that not only culminated in this work, but has also greatly enhanced my personal growth as a researcher. Professor Sarris’ supervision was the ideal combination of providing support, yet allowing for independence, and I appreciate that very much. I also want to thank Professors Triverio and Eleftheriades for their invaluable courses, and for their willingness to be on my defense committee; as well I’d like to thank the chair of my committee, Professor Herman, for his time and interest in my research. I would also like to thank my colleagues in the Electromagnetics Group, particularly Trevor Cameron, Hans-Dieter Lang, Tony Liang, Neeraj Sood and Xingqi Zhang (in alphabetical order), for their time, advice and patience in various matters ranging from typesetting to antenna design. In particular, I’d like to acknowledge their sense of humour, their tolerance for my sense of humour, and their excellent taste in coffee. I would also like to acknowledge Alex Wong for helping me understand impedance surfaces. I’d like to thank my parents and sister for being able to radiate love, support and encouragement from a distance of 11000 km. Finally, I want to thank Adriana, for inspiring and encouraging me constantly, and for putting up with my odd work hours. iv Contents 1 Introduction 1.1 1.2 1 A Review of Literature in Microwave Ablation Technology . . . . . . . . . . . . . 2 1.1.1 MWA and Other Ablation Techniques . . . . . . . . . . . . . . . . . . . . 3 1.1.2 Monopole MWA Antennas in Literature . . . . . . . . . . . . . . . . . . . 4 1.1.3 Dipole MWA Antennas in Literature . . . . . . . . . . . . . . . . . . . . . 5 1.1.4 Other MWA Antennas in Literature . . . . . . . . . . . . . . . . . . . . . 6 1.1.5 Gaps in Literature and Motivation for our Work . . . . . . . . . . . . . . 7 Overview of the Thesis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8 2 Background and Motivation 9 2.1 The Physics of Microwave Ablation . . . . . . . . . . . . . . . . . . . . . . . . . . 11 2.2 Effects of Thermal Conductivity . . . . . . . . . . . . . . . . . . . . . . . . . . . 12 2.3 Electric Field Focusing v. Temperature Focusing . . . . . . . . . . . . . . . . . . 12 2.4 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15 3 The Multiphysics Optimization Design Methodology 3.1 3.2 3.3 16 Optimization Theory . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18 3.1.1 Standard Approach: SAR-Based Optimization . . . . . . . . . . . . . . . 18 3.1.2 Multiphysics Approach: Temperature-Based Optimization . . . . . . . . . 18 3.1.3 Optimization Method Comparison . . . . . . . . . . . . . . . . . . . . . . 19 Optimization and Computation Details . . . . . . . . . . . . . . . . . . . . . . . 20 3.2.1 Sources . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20 3.2.2 Simulation Domain . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22 3.2.3 Target Profiles . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23 Results and Comparisons . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24 3.3.1 Confinement of the Ablation Zone . . . . . . . . . . . . . . . . . . . . . . 24 3.3.2 Temperature Gradient at Target Boundary . . . . . . . . . . . . . . . . . 25 3.3.3 Optimized Currents . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25 3.4 Physical Implementation of Optimized Current Elements – A Proof of Concept . 26 3.5 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29 v 4 A Printed Antenna Concept for Microwave Ablation 30 4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30 4.2 Design Considerations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30 4.3 Design Description and Optimization . . . . . . . . . . . . . . . . . . . . . . . . . 31 4.3.1 Effects of Varying Capacitance Values . . . . . . . . . . . . . . . . . . . . 34 4.3.2 Effects of Varying the Feed Position . . . . . . . . . . . . . . . . . . . . . 34 4.3.3 Design Optimization for a Test Target . . . . . . . . . . . . . . . . . . . . 35 4.3.4 Simulation Details . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35 4.4 Results and Analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36 4.5 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39 5 Conclusions 5.1 5.2 41 Significance of this Work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41 5.1.1 Clinical Integration . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42 5.1.2 Publications . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43 Future Work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43 5.2.1 More Efficient Multiphysics Optimization . . . . . . . . . . . . . . . . . . 43 5.2.2 Uncertainty Analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44 5.2.3 Antenna Reconfigurability . . . . . . . . . . . . . . . . . . . . . . . . . . . 44 5.2.4 Further Antenna Design Optimization . . . . . . . . . . . . . . . . . . . . 44 5.2.5 Fabrication and Testing . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44 Appendices 45 A Finite Differences Approximation of the Bioheat Equation 45 A.1 Electric Field Computation in 2D . . . . . . . . . . . . . . . . . . . . . . . . . . . 45 A.1.1 Discretization of Maxwell’s Equations . . . . . . . . . . . . . . . . . . . . 46 A.2 Computation of the Bioheat Equation in 2D . . . . . . . . . . . . . . . . . . . . . 50 A.3 Validation of FDFD Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50 B Material Property Estimation by Inverse Scattering 52 B.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52 B.2 The Inverse Scattering Algorithm . . . . . . . . . . . . . . . . . . . . . . . . . . . 52 B.3 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54 B.4 Relevance to Microwave Ablation . . . . . . . . . . . . . . . . . . . . . . . . . . . 54 Bibliography 55 vi List of Tables 3.1 Material properties used in simulations. . . . . . . . . . . . . . . . . . . . . . . . 22 3.2 Relative, Euclidean and infinity norm errors for comparing optimization methods. 24 4.1 Material properties used in simulations. . . . . . . . . . . . . . . . . . . . . . . . 34 4.2 Input parameters - variable capacitance test cases. Capacitance values are in pF; reactances are in Ω. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34 4.3 Port parameters - variable capacitance test cases. . . . . . . . . . . . . . . . . . . 36 4.4 Port parameters - variable feed position test cases. . . . . . . . . . . . . . . . . . 36 4.5 Optimized port parameters for each test case. . . . . . . . . . . . . . . . . . . . . 39 vii List of Figures 1.1 Cross-section of the monopole MWA antenna concept. . . . . . . . . . . . . . . . 4 1.2 Cross-section of the triaxial MWA antenna concept. . . . . . . . . . . . . . . . . 4 1.3 Cross-section of the monopole MWA antenna concept with a choke. . . . . . . . 5 1.4 Cross-section of the dipole MWA antenna concept. . . . . . . . . . . . . . . . . . 6 1.5 Cross-section of the floating-sleeve dipole MWA antenna concept. . . . . . . . . . 7 2.1 Diros RFA probes with three tines. . . . . . . . . . . . . . . . . . . . . . . . . . . 10 2.2 Diros RFA probe temperature profiles in ◦ C, with the 50◦ C contour in black. . . 10 2.3 (a) 2D SAR profile (in W/kg) and (b) corresponding temperature profile (in ◦ C) with the 50◦ C contour. kth values in W/m/K. . . . . . . . . . . . . . . . . . . . . 12 2.4 Geometry for the field and temperature focusing study. . . . . . . . . . . . . . . 13 2.5 Temperature confinement as a function of electric field beam width and thermal conductivity. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14 2.6 (a) Electric field (at −6 dB) and (b) corresponding temperature (at 50◦ C) contours on a 2D slice through the domain, for various beam widths at a particular value of thermal conductivity. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14 3.1 Overview of the optimization-based methodology. . . . . . . . . . . . . . . . . . . 17 3.2 Temperature as a function of time at four probe points. . . . . . . . . . . . . . . 21 3.3 Source geometry for optimization with indexing. . . . . . . . . . . . . . . . . . . 21 3.4 Depiction of the simulation domain and discretization into slices for optimization and visualization. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22 3.5 Left: One xy-slice of target profiles for (a) E-based and (c) T -based optimization. Right: Points over which optimization is performed (in black) for selected xyslices for (b) E-based and (d) T -based optimization. Red curves represent target boundary. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23 3.6 Temperature profiles at selected xz-slices of the 3D domain. (a) SAR-based optimization, and (b) T -based optimization. Target region boundary shown in red; 50◦ C contour shown in black. . . . . . . . . . . . . . . . . . . . . . . . . . . 25 viii 3.7 Surface plot of T for an xz-slice at y ≈ 0 cm. (a) SAR-based optimization, and (b) T -based optimization. Target region boundary shown in red. Empty region in the centre represents position of source array. . . . . . . . . . . . . . . . . . . . 26 3.8 Surface plot of normalized |E| (in dB) for an xz-slice at y ≈ 0 cm. (a) SAR-based optimization, and (b) T -based optimization. Target region boundary shown in red. Empty region in the centre represents position of source array. . . . . . . . . 26 3.9 Top row: Optimized current source element magnitudes. (a) SAR-based optimization, and (b) T -based optimization. Bottom row: Optimized current source element phases. (d) SAR-based optimization, and (e) T -based optimization. . . . 27 3.10 Proposed proof of concept design and its surface current distribution, compared to target currents. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28 3.11 (a) Temperature profile generated by the proposed design at selected xz-slices of the 3D domain. Target region boundary shown in red; 50◦ C contour shown in black. (b) Surface plot of T generated by the proposed design for an xz-slice at y ≈ 0 cm. Target region boundary shown in red. Empty region in the centre represents position of source. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29 4.1 Proposed antenna design. (a) Perspective view with dimensions and blow-up of metallization layers. (b) Front and side views including outer teflon casing. . . . 32 4.2 Simulation geometry. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33 4.3 Temperature as a function of time at four probe points. . . . . . . . . . . . . . . 36 4.4 (a) and (c): Surface current distributions. (b) and (d): Corresponding 50 ◦ C contours at one central slice. (a) and (b): variable capacitance test cases. (c) and (d): variable feed position test cases. . . . . . . . . . . . . . . . . . . . . . . 37 4.5 Optimized current distributions on one dipole. . . . . . . . . . . . . . . . . . . . 38 4.6 Optimized temperature profiles in ◦ C. Black: 50◦ C contour. Red: target boundary. 40 A.1 2D FDTD Yee cell. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46 A.2 Equivalent 2D FDTD Yee cell. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47 A.3 (a) Ez computed via homemade FDFD. (b) Ez computed in COMSOL. (c) Ez along the vertical test line shown in (a) and (b). Blue line: computed in COMSOL. Red crosses: computed via FDFD. (d) Ez along the horizontal test line shown in (a) and (b). Blue line: computed in COMSOL. Red crosses: computed via FDFD. Ez values are in V/m. . . . . . . . . . . . . . . . . . . . . . . . . . . . 51 B.1 Geometry for 2D FDTD-based inverse scattering algorithm. Small black squares: transmitter/receiver pairs. Yellow shaded square: unknown domain being imaged. 53 B.2 (a) Actual electrical permittivity profile to be imaged. (b) Imaged permittivity profile obtained via inverse scattering algorithm. Actual profile in semitransparent red. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54 ix B.3 (a) Actual electrical conductivity profile to be imaged. (b) Imaged conductivity profile obtained via inverse scattering algorithm. Actual profile in semitransparent red. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55 x Chapter 1 Introduction There are several malignant medical conditions that are caused by, or result in, the unwanted presence of cells in localized regions of the body; in such cases, often the only promising type of treatment is to destroy these cells selectively. This form of treatment is called ablation. There are several ways in which tissue destruction can be attained; one popular method is to selectively heat unwanted tissue to deleterious temperatures. This is referred to as thermal ablation therapy. For example, cancerous tumours are mutations in healthy cells, leading to the uncontrolled growth of malignant masses of tissue. In many cases, mutated cells develop resistance to traditional treatment methods such as radiation and chemotherapy, or are not easy to surgically remove from the body. In these cases, one clinically successful alternative has been to heat the tumours to temperatures high enough to destroy the malignant cells . Another example is irregular heart rhythms (arrhythmia) caused by areas of abnormal heart tissue. In order to restore normal heart rhythms, it may be necessary to destroy the abnormal tissue via heat . Yet another example is nerve endings in areas of the body inflamed due to osteoarthritis, causing the patient chronic pain. In this case, if the osteoarthritis itself cannot be cured, the patient can be spared chronic pain by heating, or ablating, the affected areas to kill the nerve cells . In all of these cases, it is important to develop and improve upon technologies that are able to deliver heat to target regions accurately and quickly, without damaging surrounding healthy tissue and organs. Ablation therapy involves rapid delivery of heat to a targeted region in tissue, inducing temperatures in excess of 50◦ C , and is generally carried out over the span of 5−10 minutes. Heat can be delivered to the target cells in several ways, including ultrasound and electromagnetic radiation. In microwave ablation (MWA), heat delivery is enabled by electric fields radiated by interstitially-applied antennas operating at microwave frequencies. MWA has clinical applications in all of the medical conditions mentioned above, most popularly for the treatment of malignant tumours [1, 4, 5]. One major challenge with MWA is designing microwave sources to selectively apply heat to a targeted region without damaging surrounding healthy tissue. Several existing ablation antenna designs attempt to achieve this confinement by focusing an electric field within the target [4,6]. 1 Chapter 1. Introduction 2 This assumes that a focused electric field or specific absorption rate (SAR) profile will lead to a focused temperature profile. However, the relationship between an applied electric field and the consequent temperature profile is generally more complex, especially when the local electrical and thermal properties of tissue are inhomogeneous. Thus, focusing or confining the electric field does not guarantee an optimal temperature profile for ablation of a thermally inhomogeneous target region. Moreover, microwave ablation antennas are generally designed in a forward approach, and are able to generate radiation patterns, and thus temperature profiles, that are specific to a particular design [2, 7–10]. However, different applications of ablation may require different shapes and sizes of temperature profiles. One major limitation of most designs in literature is their inability to generate controllable temperature profiles customized to the shape of the region being targeted. The design of a particular ablation antenna must be informed by the application for which it is needed. Thus, a single antenna design that can be easily configured for several different applications would provide a cost-effective and versatile contribution to ablation therapy in a clinical setting. The aim of this work is to tackle the shortcomings of existing MWA source design procedures, and the limitations of existing MWA antennas, by proposing a new, more intelligent inverse methodology for designing these antennas, that accounts for thermal inhomogeneities and heat diffusion effects. Subsequently, the advantages of this methodology are demonstrated by using it to design an easily configurable printed MWA antenna that addresses several limitations of existing designs. The proposed design allows the user to synthesize a source current distribution based on a desired target to ablate, offering the functionality of a phased array implemented on the tip of an interstitial probe. In the following section, a review of existing MWA technologies is presented. Advantages of MWA compared to other ablation technologies such as High-Intensity Focused Ultrasound (HIFU) and Radio Frequency Ablation (RFA) are discussed. Various MWA designs, their results, and their limitations are discussed and used to motivate the need for a new inverse methodology and more easily configurable MWA antennas. 1.1 A Review of Literature in Microwave Ablation Technology MWA is a relatively new form of therapy. The use of an interstitial microwave-based device for liver tumour therapy was first proposed by Tubuse  in 1979. The technique was subsequently applied successfully in a clinical setting by Tubuse et al. in 1985 . Although initially this technique was mostly popular in Japan and China, as of 2003, MWA antennas have been commercially sold in the United States courtesy of Vivant Medical, Inc (Mountain View, Calif., USA). There are now several commercial outlets for MWA technology, and various advancements in MWA antenna design continue to be made. Most MWA antennas are thin, needle-like coaxial-based devices, and the most popular ones are either monopole, dipole or slot Chapter 1. Introduction 3 antennas [7–10, 13], with some helical and loop-based antennas also emerging recently . In the following subsections, we first present a comparison of MWA versus other popular ablation techniques, and subsequently summarize popular MWA antenna types and discuss their advantages and disadvantages. This analysis is then used to motivate our work. 1.1.1 MWA and Other Ablation Techniques There are several alternatives to microwave ablation in terms of destroying malignant tissue in a percutaneous or noninvasive way. One of the earliest methods developed for this purpose is High-Intensity Focused Ultrasound (HIFU), first introduced by Lynn et al.  in 1942. In this method, high intensity acoustic waves at frequencies between 0.5–10 MHz are focused using acoustic lenses into a tight target spot, resulting in the deposition of heat energy in that spot. This is completely non-invasive technique that has gained significant popularity and garnered many improvements over the years. HIFU has been successful in accurately targeting tumours in clinical trials , however its success rate is not statistically different from MWA , and modern HIFU applicators generally require the optimization of hundreds or even thousands of acoustic source elements to properly focus the energy, making it a more expensive technology. Also, successful HIFU ablation requires significantly longer treatment time than MWA . Another popular ablation technique that is closely related to MWA is Radio Frequency Ablation (RFA). RFA is similar to MWA in that thin needle-like applicators are fed with alternating currents to generate heat interstitially. However, in RFA, lower frequencies are used as compared to MWA (300–500 kHz), and the mechanism by which tissue is heated is quite different. An RFA applicator consists of a needle-like structure with a metallic section to which an alternating voltage is applied, and a ground that is either placed interstitially on the same needle, or is placed on the interface between skin and air. This forms a current path through the tissue near the applicator, which generates heat in the immediate vicinity of the needle. The heat then diffuses through the nearby target tissue to create an ablation zone . Although RFA has been successfully applied in a clinical setting and is commercially available, due to its heating mechanism it is unable to generate ablation zones as large as in MWA. It also requires a conducting path through tissue to be formed, which may create unnecessary heating in non-targeted regions, since the conducting path is difficult to control. Also, certain types of tissue are not conducive to the passage of electric current, and thus are not suitable for RFA. This limits its applications and makes resulting ablation zones unpredictable. In order to manipulate the shapes of temperature profiles, RFA applicators are often designed with sharp tips or multiple tines, which is quite invasive [3,18]. These limitations are quite easily overcome by MWA [19–22]; MWA is able to produce larger ablation zones more quickly, without relying on any conduction paths through tissue. Although each of these ablation methods has its advantages and have been in use longer than MWA, the use of microwaves offers many of the same benefits while overcoming several of the disadvantages of HIFU and RFA. 4 Chapter 1. Introduction 1.1.2 Monopole MWA Antennas in Literature These antennas are constructed from thin, semirigid coaxial cable, with the inner conductor made to extend beyond the tip of the outer conductor by stripping the outer conductor , as shown in Fig. 1.1. In , three types of monopole antennas are considered; each one is a variation in how far the dielectric extends beyond the tip, and the metallization on the tip. Simulation and experimental results were presented for each case, and in general, the temperature patterns obtained in each case were either nonuniform, or not large enough for practical application. Moreover, the temperature profiles were elongated along the tips of the antennas due to backwards current flow. Although the size of the patterns can be somewhat modified by varying the length of the extended tips or the input power, the poor return loss and nonuniform temperature profiles are major drawbacks. These designs do not allow the user control over the shape of the temperature profiles that are produced, and are not tunable to local electrical or thermal properties of tissue. Coaxial feed, outer conductor Monopole antenna Coaxial feed, inner conductor Teflon Figure 1.1 Cross-section of the monopole MWA antenna concept. More promising results were presented by Brace et al. , where a coaxial monopole antenna is placed inside a biopsy needle, to effectively yield a triaxial antenna, as shown in Fig. 1.2. This design greatly reduces the return loss and is able to generate much larger ablation zones (≈ 3 cm in diameter) at moderate input powers (≈ 50 W). However, this design also offers no control on the shapes of generated temperature profiles; the profiles are elliptical, whereas in the treatment of tumours, spherical profiles are more sought after. It cannot be configured to suit a particular application, since there are very few design parameters that can be varied, that would have an appreciable effect on the shape of the temperature profiles. Coaxial feed, outer conductor Monopole antenna Biopsy needle Coaxial feed, inner conductor Teflon Figure 1.2 Cross-section of the triaxial MWA antenna concept. A modification of the above was presented by Prakash et al.  in 2009, where rather than 5 Chapter 1. Introduction a biopsy needle, the outer conductor of a coaxial monopole antenna was fitted with a quarterwavelength metallic choke, as shown in Fig. 1.3. The quarter-wavelength nature of the choke prevents flow of current along the outer surface of the outer conductor of the coax feed, thus reducing the backward flow of current and subsequent unwanted heading along the feed line. This leads to larger and more spherical ablation zones. The authors optimize the parameters of the choke-based antenna to obtain near-spherical ablation zones. However, their simulation and experimental results both indicate somewhat elongated ablation zones, and significant unwanted heating along the feed, although less so than previous monopole designs. This design also does not offer control of the shape of the ablation zone. Coaxial feed, outer conductor Monopole antenna λ/4 Choke Coaxial feed, inner conductor Teflon Figure 1.3 Cross-section of the monopole MWA antenna concept with a choke. 1.1.3 Dipole MWA Antennas in Literature One of the earliest types of MWA antenna was a coaxial-based dipole, also built with a thin semirigid coaxial cable, where a coaxial transmission line ends into a gap, with the inner conductor connected across the gap to a metallic extension on the other side of the gap, as shown in Fig. 1.4. Jones et al.  presented this classical dipole antenna for MWA in 1989; they also investigated the use of multiple dipoles inserted to form an array. They were able to generate large (over 4 cm diameter) ablation zones using up to 8 or 9 antennas, and concluded that the complex SAR profiles generated did not seem immediately useful in the clinic, and may require further analysis with even more complex arrangements and variations in phase. These results are discouraging, since inserting more than one antenna is quite invasive, and since then, various designs have been proposed that achieve better results without having to use multiple antennas. These designs also suffer from serious drawbacks: no temperature profiles are presented and thermal properties are ignored, and the SAR or temperature profile shapes are not controllable in any way. Hürter et al.  proposed a modification of this type of antenna in 1991. The authors propose a sleeve along the coax feed line connected to the outer conductor that prevents the backward flow of current, thus ensuring that most of the power is deposited at the tip of the antenna (near the gap) rather than along the feed line. The authors compare this design to the more standard dipoles without a sleeve, and are successful in showing that the electric field is indeed more concentrated near the dipole gap, rather than along the feed line, leading to tighter 6 Chapter 1. Introduction Coaxial feed, outer conductor Coaxial feed, inner conductor Extension Dipole gap Teflon Figure 1.4 Cross-section of the dipole MWA antenna concept. ablation zones. Unlike other designs without such a sleeve, the radiation and temperature patterns of this design do not depend on the insertion depth of the antenna, which is a major advantage. However, once again all of the analysis provided by the authors is in the electrical domain; only SAR patterns (which depend on the applied electric field, electric conductivity and tissue density) are presented, with no mention of corresponding temperature profiles. Thus, thermal properties and heat diffusion effects are not accounted for. Also, as in the cases above, these designs produce a fixed SAR pattern, with little control over the shape of the generated fields or temperature profiles. More recently, a modification of Hürter’s sleeve-based dipole antenna has been presented, where the sleeve is floating rather than fixed. This was first presented by Yang et al.  in 2006. A coaxial dipole antenna, similar to the ones presented above, is encased in a teflon coating. A metallic sleeve is slid onto the teflon cylinder at a variable position, such that it is electrically isolated from the coax cable, unlike the design presented by Hürter. This is depicted in Fig. 1.5 The floating nature of the sleeve allows an added degree of freedom in terms of design parameters that can be tuned to generate a suitable temperature profile. Prakash et al.  subsequently presented an optimization procedure for the various design parameters of the floating sleeve MWA antenna to produce the largest ablation zones. Both works yielded very encouraging results, with significant reduction in the unwanted heating of the feed line, and the ablation zones were quite successfully localized near the antenna tips. However, in both studies, once again only SAR patterns are discussed. The fact that SAR patterns are well localized to antenna tips does not imply that the temperature profile will also be thus confined. Also, the various design parameters that can be adjusted in this case only modify the size and extent of the SAR patterns, not their shapes. Although these designs yield large ablation zones with reflection coefficients well below −15 dB, they do not offer control and specificity to particular applications. The generic SAR profiles may be useful in some applications, but not in others. 1.1.4 Other MWA Antennas in Literature A loop antenna was presented by Shock et al. , where a probe is inserted into the target region, and a metallic loop is deployed with the help of electrocautery (i.e. the heating up of the microwave probe) to cut through the tissue. The authors also investigated the use of multiple 7 Chapter 1. Introduction Coaxial feed, outer conductor Extension Floating sleeve Coaxial feed, inner conductor Dipole gap Teflon Figure 1.5 Cross-section of the floating-sleeve dipole MWA antenna concept. loops; experimental results indicate that these antennas are able to produce high temperatures within the ablation zones, and generate fairly spherical temperature profiles. However, this method suffers from several disadvantages; the loops are only effective if they are accurately placed such that they perfectly encircle the target region. Thus, a particular loop can only be applied to a specific target type and size. Also, deployment of the loop requires electrocautery surgery which is very invasive, and creates the risk of harming healthy, important nerve cells or surrounding tissue. Also, much like aforementioned designs in literature, there is no control over the shapes of generated temperature profiles. Luyen et al.  very recently presented a coaxial-based helical antenna, which avoids the use of a balun or choke of any kind to prevent back flow of current. However, the results indicate elongated, uncontrollable temperature profiles with significant tail-end heating of the feed line. The advantages of this design include a return loss of −30 dB and the fact that the lack of a balun allows the antenna to have a lower diameter than other designs; however, the disadvantages in ablation zone shapes do not sufficiently justify the absence of a balun. 1.1.5 Gaps in Literature and Motivation for our Work The various designs that have been proposed in literature have shown great improvement over the years, and have collectively proved that MWA is a promising and viable clinical technique for various applications. These designs provide several advantages over other ablation methods such as RF and HIFU; however, several issues exist among these designs, and there is significant room of improvement in the field. As we claimed earlier, one common lacking feature among these designs is that they produce generic temperature profiles whose shapes cannot be controlled or modified; this is a fundamental limitation of the fact that these antennas are designed in a forward approach with no specific target profile in mind. Another common flaw is that many of these designs focus on the electric field or SAR profiles generated, but do not consider heat diffusion around the antennas, or the local thermal properties of tissue. These are the two primary areas of improvement that are targeted in this thesis. Chapter 1. Introduction 1.2 8 Overview of the Thesis The thesis is organized as follows: first, we provide background regarding our initial work in electromagnetic ablation. The physics of MWA is discussed, and the relations between the thermal and electrical domains are studied through simulations to motivate our investigation into multiphysics methods and the need for an inverse antenna design approach. A methodology for the optimization of MWA sources is then presented, and the standard electric field-based approach is compared with the proposed multiphysics approach. The results obtained via each approach are then analyzed and discussed. A proof-of-concept MWA source design is discussed. The multiphysics methodology is then put to practice in order to design a printed microwave ablation antenna, that can be configured to synthesize a range of temperature profile shapes and sizes. The antenna provides the functionality of a phased array in allowing control over its current distribution, which can be optimized through simple physical and electrical parameters to generate a desired temperature profile. The performance of the antenna is analyzed to assess its feasibility in a clinical setting. The work is concluded with a discussion of the significance and relevance of the proposed ideas, as well as future directions that may be taken to bring these ideas into a clinical setting. Chapter 2 Background and Motivation Our exploration into electromagnetic ablation methods began with performing multiphysics simulations of RFA probes designed by Diros Technology Inc., Markham, ON, Canada. During the collaboration we simulated the electric field and temperature profiles of new RF probes that were experimentally validated at Diros, and are now commercially available . The simulations were set up in the commercial finite element method solver, COMSOL Multiphysics (COMSOL Inc., Burlington, MA). The relation between the electric and thermal domains is modeled by the Pennes Bioheat Equation  (see below for details), and the probes were modeled in the time domain to evaluate thermal profiles after 10 minutes of application. Fig. 2.1 shows one of the new RFA probes proposed by Diros. The reasoning behind using three sharpedged tines is based on the conventional wisdom that sharp corners lead to higher local electric field intensities, leading to higher temperatures around the probes. The simulated temperature profiles, along with the 50◦ C contours, are shown in Fig. 2.2. The temperature profile is extended near the tips due to the three tines. The validity of the simulations was confirmed by experimental results provided by Diros. The probes were inserted into bovine liver, which was sliced post-therapy to study the ablation zones. The fundamental design ideology behind these, and several other RFA probes is that the physical shapes of the probe tips are manipulated in order to control the shape and extent of the ablation zone. However, the sharp-edged nature of these designs makes them quite invasive. Additionally, if the ablation zone shapes and sizes were to be manipulated via the electrical properties of the probe rather than its physical characteristics, it could allow for more versatile, simpler, smaller and less invasive designs. This was one of our initial motivators for looking into a more intelligent, inverse way to design ablation devices. Upon our investigation into electromagnetic ablation methods, the advantages of microwavebased techniques over RF probes (outlined in the previous chapter) led us to apply the idea of an inverse design procedure to MWA antennas, rather than RFA probes. In order to develop a more intelligent methodology for the design of MWA antennas, it is necessary to study the physics of MWA, and the factors that relate the applied electric field to the consequent temperature profile. In the following sections, through simulations, we study the relations between the 9 10 Chapter 2. Background and Motivation (a) (b) (c) Figure 2.1 Diros RFA probes with three tines. (a) (b) Figure 2.2 Diros RFA probe temperature profiles in ◦ C, with the 50◦ C contour in black. Chapter 2. Background and Motivation 11 electrical and thermal domains, and understand the influence of thermal properties and heat diffusion. 2.1 The Physics of Microwave Ablation In MWA, the ablation zone is generated via a combination of heat energy deposited by the applied electric field, and the diffusion of that deposited energy through tissue . The electric field, oscillating at microwave frequencies, causes the water molecules to constantly rotate to align their polarity with the field. This increases the local kinetic energy, thus generating heat and resulting in ablation if the temperature is high enough [1,19]. Since the relative permeability of biological tissue is generally 1, the complex dielectric permittivity determines the transmission of microwave energy into tissue. The complex permittivity is treated as a relative permittivity that represents the real part, and an electric conductivity that represents the imaginary part. In microwave ablation, electrical energy is primarily transmitted into tissue via displacement currents generated by the applied field. This is different from RF ablation, where the energy is primarily deposited in an ohmic fashion via a conductive path formed through tissue between the source and ground. Thus, in RF ablation, the applicator is essentially a probe, while in microwave ablation, the applicator is an antenna . The mathematical relation between an applied electric field E and resulting temperature profile T in biological tissue is well approximated by the Pennes Bioheat Equation . ρCp ∂T 1 = σ|E|2 + ∇ · (k∇T ) + Wb Cb (Tb − T ) ∂t 2 (2.1) where ρ, Cp , k and σ are the density, heat capacity, thermal conductivity and electrical conductivity of tissue, and Wb and Cb are the perfusion rate and heat capacity of blood. The first term on the right represents electrical energy deposition in tissue; the second term on the right represents heat diffusion through tissue; and the third term represents the cooling effect of blood perfusion through tissue. Another quantity relevant to MWA is the specific absorption rate (SAR), which is defined as follows: SAR = σ |E|2 2ρ (2.2) The SAR is an indication of energy deposition in tissue due to an applied electric field, that takes into account tissue density. It is clear that the temperature profile depends on both, the applied electric field, and the diffusion of heat through the tissue, which in turn depends on the local thermal conductivity. Higher values of k imply increased diffusion of heat around the region of deposition. Thus, it is clear that the shape and size of the ablation zone depend on two factors: heat deposition (via the strength and pattern of the applied electric field); and heat diffusion (via the thermal properties and inhomogeneities of the target region). 12 Chapter 2. Background and Motivation 0.01 kth = 0.9 60 kth = 0.9 kth = 0.1 1500 0 y (m) 2000 0.005 y (m) 0.01 2500 kth = 0.1 0.005 55 0 50 1000 45 -0.005 -0.005 th th -0.01 -0.01 -0.005 500 k = 0.7 k = 0.3 0 0.005 0.01 k = 0.7 k = 0.3 th th -0.01 -0.01 -0.005 x (m) (a) 0 0.005 40 0.01 x (m) (b) Figure 2.3 (a) 2D SAR profile (in W/kg) and (b) corresponding temperature profile (in ◦ C) with the 50◦ C contour. kth values in W/m/K. 2.2 Effects of Thermal Conductivity The influence of thermal conductivity on a temperature profile resulting out of an applied electric field can be demonstrated by a simple simulation of (2.1) in a domain with varying thermal conductivities. A 2D square domain with four quadrants of exactly the same electrical and thermal properties, except thermal conductivities that range from 0.1 to 0.9 W/m/K, is simulated. The source is assumed to be a circular surface current in the centre of the domain, flowing in the direction perpendicular to the 2D plane (ẑ). The resulting SAR and temperature profiles are shown in Fig. 2.3. As expected, the SAR profile in Fig. 2.3a is symmetric about the source current, since the electrical properties and densities of each of the four materials is the same. However, it is clear from Fig. 2.3b that the corresponding temperature profile is significantly distorted due to varying thermal conductivities. This occurs, as expected, because the rates of thermal conduction through each of the four materials is significantly different. This clearly indicates the need to design and study MWA antennas in the context of the temperature profiles and not just the SAR profiles they generate. This is important to allow accurate ablation of the targets, and to minimize damage to healthy surrounding tissue. 2.3 Electric Field Focusing v. Temperature Focusing In order to develop a more intelligent design methodology for MWA antennas, it is necessary not only to study and quantify the roles of heat deposition and heat diffusion, but also to study the fundamental limits on the ability to control and shape temperature profiles. To this end, and to motivate our investigation into multiphysics methods, a study was conducted to analyze 13 Chapter 2. Background and Motivation the extent to which focusing of an electric field allows focusing of the temperature profile, and the limits beyond which heat diffusion effects take over. The purpose of this initial study was to, through simulations, understand and quantify the extent to which each of these factors (heat deposition and heat diffusion) contributes towards the temperature distribution. In this study, a dipole antenna array with 5 × 5 elements placed near human liver tissue was simulated to focus electric fields within the tissue to varying degrees, and the resulting temperature profiles were computed using (2.1). The computations were carried out using COMSOL Multiphysics. The geometry is depicted in Fig. 2.4. The excitation of each element was chosen via the Dolph-Tschebyscheff method  to control and vary the maximum beam width of the generated electric field profile. The beamwidth is controlled by specifying the desired side lobe level, RdB , in dB. The excitations an are chosen using the following equations : an = M +1 X (−1)M −q+1 (z0 )2(q−1) q=n 1 z0 = 2 " (q + M − 2)!(2M ) εn (q − n)!(q + n − 2)!(M − q + 1)! # q 1/P q 2 R0 + R0 − 1 + R0 − (2.3a) 1/P R02 − 1 R0 = 10RdB /20 (2.3b) (2.3c) If N is the number of array elements along one direction, then M = (N − 1) / (2), n = 1, 2, . . . M + 1, ε equals 2 if n = 1 and 1 otherwise, and P = N − 1. In this way we can find the excitations required along each dimension. The excitation vectors can then be matrix-multiplied to obtain the excitation of each element in the 2D array. Liver Air Dipole array Figure 2.4 Geometry for the field and temperature focusing study. 14 Chapter 2. Background and Motivation T profile (kth = 0.10) Diameter of 500C isotherm (cm) 3.6 T profile (kth = 0.30) T profile (kth = 0.57) 3.55 T profile (kth = 0.75) T profile (kth = 0.95) 3.5 3.45 3.4 3.35 3.3 3.25 2.95 3 3.05 3.1 3.15 3.2 3.25 Diameter of -6 dB electric field contour (cm) 2 2 1.5 1.5 z (cm) z (cm) Figure 2.5 Temperature confinement as a function of electric field beam width and thermal conductivity. 1 0.5 1 0.5 0 0 -2 -1 0 1 2 -2 -1 0 x (cm) x (cm) Sources Sources (a) 1 2 (b) Figure 2.6 (a) Electric field (at −6 dB) and (b) corresponding temperature (at 50◦ C) contours on a 2D slice through the domain, for various beam widths at a particular value of thermal conductivity. The maximum diameter of the 50◦ C isotherm at steady state is reported in Fig. 2.5 as a function of the maximum electric field beam width at −6 dB, for different thermal conductivities of tissue, in order to understand the contributions of deposited electric energy and heat diffusion Chapter 2. Background and Motivation 15 to the obtained temperature profiles. Fig. 2.6a and Fig. 2.6b are examples of a slice of the electric field and corresponding temperature contours (at −6 dB and 50◦ C) obtained for a particular value of thermal conductivity, for different beam widths. Two main conclusions can be drawn from Fig. 2.5 and Fig. 2.6: • There is a fundamental limit to which the temperature profile can be focused or confined, regardless of how focused the electric field is, for a given amount of energy deposited. This limit arises due to heat diffusion around the area where the electric field is confined; the energy deposited via an electric field will diffuse around the region even if the electric field itself is confined to a small region. This must be taken into account when attempting to ablate a given target region. • The 50◦ C isotherm diameter is inversely related to the thermal conductivity. This can be understood as follows: although a lower thermal conductivity implies less heat diffusion, this reduced diffusion causes an accumulation of heat near the source, thus raising the overall temperature around the source. Given that biological tissue may be thermally inhomogeneous, and that tumours can exhibit thermal conductivities ranging from 0.3 to 0.7 Wm−1 K−1 , , the significant relation between temperature and thermal conductivity must be taken into account in designing antennas for ablation. The observations made above clearly demonstrate the complex relation between applied fields and consequent temperature profiles, and indicate the necessity for taking these complexities into account when designing MWA antennas. 2.4 Summary The purpose of this chapter was to establish the need for an inverse multiphysics approach to designing MWA antennas. Through simulations, the relation between the electric field and temperature profiles were studied, as well as the effect of thermal conductivity. The results of these studied showed that the thermal profile depends significantly on both heat deposition in the tissue, and heat diffusion through tissue. It was also established that simply focusing or shaping the electric field or SAR profile is not sufficient to guarantee optimal generation of a desired ablation zone. In the next chapter, an inverse multiphysics design approach is described and analyzed. Chapter 3 The Multiphysics Optimization Design Methodology With the results of the aforementioned studies and the understanding of heat deposition and heat diffusion, we are now in a position to propose an inverse, optimization-based methodology for MWA antenna design. One heat-based therapeutic technique that has been used successfully in clinical applications is hyperthermia, where an antenna array placed outside the body is used to focus the electric field at a target within the body . This is achieved by appropriately picking the excitation phase and magnitude of each element of the array. Thus, it is an inverse technique, where the source is optimized in order to confine the electric field or specific absorption rate (SAR) at the target. In , an antenna array was proposed that allows sub-wavelength focusing and shaping of the magnetic field in the near-field of the array via source optimization. This technique can also be applied to shaping the electric field in the near-field. However, little attention has been given to applying this inverse approach of source optimization in designing interstitial antennas for microwave ablation. The purpose of this part of the work is to apply the concept of source optimization to the design of interstitial MWA antennas, and subsequently to develop a methodology for the design of these antennas to synthesize a desired temperature profile based on some target. Since the generated temperature profile depends on the applied electric field, it also indirectly depends on the current distribution on the antenna that produces the electric field. Thus, the approach taken in this work is to optimize the current distribution on the antenna based on a desired electric field or temperature profile. Fig. 3.1 outlines the methodology of this approach: given a target to ablate, a corresponding target temperature profile, Ttarget , and target electric field profile, Etarget , are defined. The ablation antenna is modeled as an array of surface current elements placed inside the target region, such that the current distribution on the elements generates an electric field E and corresponding temperature profile T . The current element magnitudes and phases are then 16 Chapter 3. The Multiphysics Optimization Design Methodology 17 optimized such that the generated field or temperature profile best approximates the target profile. Rather than manipulating the physical shape of the antenna, we use the properties of antenna arrays to, through their near fields, control the generated temperature profiles. This allows for the possibility of designing antennas that can generate customized temperature profiles for different applications. As mentioned previously, the standard approach in designing MWA antennas or optimizing arrays for hyperthermia involves focusing or confining the electric field or SAR within the target, which does not account for the dynamics of heat diffusion and thermal inhomogeneities in the target region. Thus, to optimize the current elements, a multiphysics approach is employed that accounts for thermal inhomogeneities in tissue. The results are compared to the standard approach of simply generating strong and confined electric fields or SAR profiles, to confirm the claim that a purely field-based optimization procedure does not guarantee an optimal temperature distribution for a given target in biological tissue. Although it is intuitively apparent that a temperature-based multiphysics optimization approach would work better than a purely field- or SAR-based one, it is useful to quantify the differences between the two approaches and to ascertain whether the multiphysics approach is worth the added computational costs. In many cases where ablation is necessary, particularly for the treatment of tumours, it is important for MWA antennas to generate symmetric, spherical ablation zones . Thus, for the purpose of this work, we focus on designs that are symmetric about the axis of the antenna. The current distribution on the ablation antenna can be thought of as a linear array of discrete current elements with magnitudes I (i) and phases β (i) , where i is an index to each element. Our aim is to optimize I (i) and β (i) such that they produce an electric field profile E, and consequently a temperature profile T , that best approximates a given target temperature profile, Ttarget . In optimizing the current elements, we compare two optimization approaches: the standard electric field-based approach where the elements are optimized such that the SAR is confined to within the target (SAR-based optimization), and a multiphysics approach based on synthesizing Ttarget directly (T -based optimization). Study the relation between E and T in tissue Identify Ttarget , and corresponding Etarget Model microwave ablation antenna as discrete surface current elements Study and compare E-based and T -based optimization methods I Optimize surface current elements to generate Etarget and Ttarget Pick an optimization method based on results Realize optimized current distribution by designing an appropriate antenna Figure 3.1 Overview of the optimization-based methodology. Chapter 3. The Multiphysics Optimization Design Methodology 3.1 3.1.1 18 Optimization Theory Standard Approach: SAR-Based Optimization In this approach, a locally high electric field and thus SAR is generated within the target, such that energy deposition is minimized outside the target. The SAR is defined as follows: SAR = σ |E|2 2ρ (3.1) Since we assume no prior information about the actual field intensities required to induce ablative temperatures, the optimization is performed in two steps. • First, we consider normalized SAR only, and attempt to shape the pattern such that the −6 dB SAR isoline lies within the target region, and decays below −10 dB near the target boundary. Based on these constraints, we define a normalized target SAR distribution, SARtarget . I (i) and β (i) are then optimized by minimizing the following objective function: fE = m,n SARm,n norm − SARtarget X 2 (3.2) m,n∈K where K is the set of indices over which the optimization is performed and SARm,n norm is the normalized SAR distribution generated by the source that consists of I (i) and β (i) . • At this point, an appropriate scaling of the optimized current magnitudes is required to induce ablative temperatures within the target. Thus, we introduce a scaling factor b that is multiplied with the optimized current magnitudes. The scaling factor itself is optimized such that the resulting temperature at the boundary of the target approaches 50◦ C. This is achieved by minimizing the following cost function: fb = X (T m,n − 50)2 (3.3) m,n∈B where B is the boundary of the target region. The updated optimized source current magnitudes are bI (i) , and the optimized phases are still β (i) . 3.1.2 Multiphysics Approach: Temperature-Based Optimization In this approach, a target temperature profile, Ttarget , is defined such that the temperature is above 50◦ C within the target region, close to 50◦ C at the target boundary, and decays to normal body temperature (37◦ C) elsewhere. The goal is to optimize the source currents to generate a temperature profile that approaches Ttarget . To achieve this, I (i) and β (i) are optimized by Chapter 3. The Multiphysics Optimization Design Methodology 19 minimizing the following objective function: fT = X m,n T m,n − Ttarget 2 (3.4) m,n∈K where K is the set of indices over which the optimization is performed. No assumptions or conditions are placed on the corresponding electric field required to generate Ttarget . This method involves computing first the field profile, and then the consequent temperature profile via (2.1) at each iteration of the optimization. 3.1.3 Optimization Method Comparison The results obtained via the two optimization methods described above are to be compared on the basis of two sets of metrics: Confinement of the Ablation Zone To quantify how closely the 50◦ C contour matches the target boundary, we define three comparison metrics: the L1 norm f1 , the Euclidean norm f2 and the infinity norm f∞ with respect to the achieved and target temperatures at the target boundary: f1 = X T m,n − T m,n target (3.5a) m,n∈B f2 = X m,n T m,n − Ttarget 2 (3.5b) m,n∈B m,n f∞ = max T m,n − Ttarget m,n∈B (3.5c) where B refers to the target boundary. Lower values of f1 , f2 and f∞ indicate better optimization results. Since the errors are computed only on the target boundaries, the number of points over which the errors are calculated is the same for each optimization method, thus making for a meaningful comparison. Temperature Gradient at Target Boundary An important consideration in ablation is to confine the heat damage to the target region, and minimize the damage to surrounding tissue. Thus, it is important to quantify the temperature gradient at the target boundary and the rate at which it decays to normal body temperature. To quantify the rate at which T drops below ablative levels outside the target, we define gB as the x-component of the gradient of the optimized temperature profile, evaluated at a single point on the target boundary. We also compare gmax , the maximum of the x-component of the 20 Chapter 3. The Multiphysics Optimization Design Methodology gradient of the optimized temperature profile on a line along x̂ in the center of the target. gB = ∇T · x̂|(xb ,zb ) (3.6a) gmax = max (∇T · x̂) (3.6b) line where (xb , yb ) ∈ B and B is the set of points corresponding to the target boundary. 3.2 Optimization and Computation Details COMSOL Multiphysics was used to solve the forward problem of computing electric fields and temperature profiles. This was interfaced with Matlab (The MathWorks Inc., Natick, MA) for the optimization procedure and for visualization of results. In both optimization approaches, a sequential quadratic programming (SQP) algorithm was employed. The function space is assumed to be convex, to allow the use of a constrained, convex optimizer. However, since the function space has 2i dimensions, convexity is not guaranteed, and it may be better to use a global optimization method. However, since the purpose here is to compare two different objective functions for optimization, a convex optimizer is used to conserve computational resources. A frequency of 2.45 GHz was used in all simulations; this is one of the clinical standards for microwave ablation . Electric fields are computed in the frequency domain. In order to reduce computational cost, it is beneficial to solve (2.1) at steady state using a stationary solver. To establish the validity of this approach, the sources are simulated to compute temperatures in the time domain. The temperatures were probed at four points at varying distances (r) from the antenna, and plotted as a function of time (t) to study the time domain characteristics. The results are shown in Fig. 3.2. It is clear that steady state in the thermal domain is achieved at t ≈ 5 minutes, and most ablation treatments last at least that long. Thus, computing temperatures at steady state is justified for the purposes of source optimization. Although tissue properties exhibit nonlinear changes at high temperatures, especially above of 100◦ C, these effects are not considered in this work since regions with such high temperatures are relatively small and confined only to the core of the target. Thus the accuracy of the results obtained are sufficient for the purposes of demonstrating the inverse design procedure and for comparing optimization methods. 3.2.1 Sources The ablation antenna is divided into six independent surface current source elements. Each surface current element is a cylindrical surface of radius a = 1 mm and height 2 mm, and all six elements are arranged in a linear array as shown in Fig. 3.3. The current elements are not included as part of the domain in which fields and temperatures are computed. The surface Chapter 3. The Multiphysics Optimization Design Methodology 21 160 r = 0.4 cm r = 0.7 cm r = 1.1 cm r = 2.1 cm 140 T ( °C) 120 100 80 60 40 2 4 6 8 10 t (min) Figure 3.2 Temperature as a function of time at four probe points. (i) current density of the ith element, Js , is defined as Js(i) = I (i) jβ (i) e ẑ 2πa (3.7) where I (i) and β (i) are the magnitudes and phases of each current element as defined earlier. The current elements are encapsulated by an outer layer of teflon (PTFE), a biocompatible material commonly used as the outer sheath in microwave ablation antennas , . i=6 i=5 i=4 2.0 mm i=3 2.6 mm i=2 3 mm i=1 Figure 3.3 Source geometry for optimization with indexing. 22 Chapter 3. The Multiphysics Optimization Design Methodology 3.2.2 Simulation Domain A cylindrical three-dimensional inhomogeneous domain was used for all simulations, as shown in Fig. 3.4a. The cylinder has a radius of 3.7 cm, and a total height of 5.6 cm. The target region to be ablated is a tumour in human liver, which is approximately a distorted sphere of radius 1 cm. The electrical and thermal material parameters used for all simulations are given in Table 3.1 [3, 38]. The domain is terminated on all sides with a perfectly matched layer (PML) . Table 3.1 Material properties used in simulations. Liver Tumour Teflon Blood Relative Relative Electrical Thermal permittivity permeability conductivity conductivity εr µr σ (S/m) kth (W/m/K) Heat capacity C (J/kg/K) Density ρ (kg/m3 ) 43.0 54.9 2.1 3540 3540 1143 3617 1079 1079 1806 1050 1 1 1 1.69 1.99 1 × 10−25 0.52 0.54 0.27 Since a 3D domain is used, the visualization of data is done on a slice-by-slice basis. The electric fields and temperature profiles are interpolated into a 3D grid of 150 × 150 × 150 cells. In other words, we consider 150 slices of 150 × 150 cells in the xz-plane. The size of each grid Liver Target 4 cm .. . .. . Sources (a) (b) Figure 3.4 Depiction of the simulation domain and discretization into slices for optimization and visualization. 23 Chapter 3. The Multiphysics Optimization Design Methodology 2 y (cm) normalized |E| 1 0.5 1 0 -1 -2 0 2 1 0 -1 -2 5 4 3 2 2 x (cm) z (cm) x (cm) 0 (a) 5 4 3 5 4 3 2 1 z (cm) (b) 2 y (cm) 55 50 T ( °C) -2 45 1 0 -1 40 -2 2 1 0 -1 x (cm) (c) -2 5 4 3 2 z (cm) 2 0 x (cm) -2 2 1 z (cm) (d) Figure 3.5 Left: One xy-slice of target profiles for (a) E-based and (c) T -based optimization. Right: Points over which optimization is performed (in black) for selected xy-slices for (b) E-based and (d) T -based optimization. Red curves represent target boundary. cell is approximately 0.35 × 0.35 × 0.37 mm. This is depicted in Fig. 3.4b. 3.2.3 Target Profiles The target field profile for E-based optimization, Etarget , and the target temperature profile for T -based optimization, Ttarget , are defined over the entire domain. However, optimization is only performed over a set of coordinates in the vicinity of the target boundary, as described below. E-based Optimization We define a normalized Etarget such that it has a value of 1 within the target, and 0 outside. Fig. 3.5a shows Etarget and Fig. 3.5b shows the points over which optimization is performed. These points are picked concentric to the target, both inside and outside of it. Chapter 3. The Multiphysics Optimization Design Methodology 24 T -based Optimization In the ideal case, the optimized temperature profile would be at or above 50◦ C within the target, and would decay to 37◦ C just outside. The rate at which temperature decays outside the target is a function of diffusion as well as energy deposition via the electric field. In the ideal case, there would be no electromagnetic energy deposition outside the target, and so the temperature profile at and outside the target boundary would be dictated purely by diffusion. Thus, Ttarget is defined such that it has a value of 50◦ C at the target boundary, and decays to 37◦ C outside due to diffusion only, in accordance with (2.1) with |E| set to 0. Fig. 3.5c shows Ttarget and Fig. 3.5d shows the points over which optimization is performed. These points are picked concentric to the target, at the boundary and outside of it. 3.3 Results and Comparisons The results are presented and discussed here from three points of view: the confinement of the 50◦ C contour to the target boundary, the T gradients achieved at the target boundary, and the optimized current distributions obtained via each optimization method. 3.3.1 Confinement of the Ablation Zone Fig. 3.6 is a 3D depiction of the optimized temperature profiles for selected slices on the xzplane, for both optimization methods. In each case, the 50◦ C contour is well confined to the boundary of the target at all slices. However, in the T -based case the 50◦ C contour is closer to the desired profile. Table 3.2 lists the comparison metrics f1 , f2 and f∞ for each method. As expected, the T -based multiphysics method performs better than the purely SAR-based method. In terms of the relative error, Euclidean norm and infinity norm, the T -based method provides a 39.3%, 70.1% and 67.9% improvement over the SAR-based method. Table 3.2 Relative, Euclidean and infinity norm errors for comparing optimization methods. f1 f2 f∞ gB gmax SAR-based optimization T -based optimization 4.96 × 104 2.30 × 105 18.2 30.8 ◦ C/cm 8.06 × 103 ◦ C/cm 3.01 × 104 6.87 × 104 5.85 44.8 ◦ C/cm 2.22 × 104 ◦ C/cm Chapter 3. The Multiphysics Optimization Design Methodology (a) 25 (b) Figure 3.6 Temperature profiles at selected xz-slices of the 3D domain. (a) SAR-based optimization, and (b) T -based optimization. Target region boundary shown in red; 50◦ C contour shown in black. 3.3.2 Temperature Gradient at Target Boundary Fig. 3.7 shows a surface plot of the optimized T profiles for each method for one xz slice at y ≈ 0 cm. It is clear that the temperature falls off with a much sharper gradient in the T -based case than in the SAR-based case, and thus is less likely to cause damage to surrounding healthy tissue. Sharper T gradients allow for more precise targeting. Table 3.2 lists the comparison metrics gB and gmax to quantify these gradients. At the boundary of the target, along the x̂ direction, T -based optimization provides a 45.5% sharper gradient. The maximum slope along the the x̂ direction is 2.75 times steeper for the T -based case. Fig. 3.8 shows the normalized SAR profile in dB for the same slice as above. It is interesting to note that the electric field is more strongly confined within the target in the T -based case than the SAR-based case. This contributes to the fact that sharper T gradients are achieved in the former, since the T profile outside the target is influenced more by diffusion of heat than deposition of energy via the electric field. This is a direct result of the target profiles that were chosen in each optimization method. In fact, the importance of heat diffusion is part of the reason that thermal inhomogeneities in tissue play an important role in the temperature distribution, and must be taken into account in addition to the electric field distribution. 3.3.3 Optimized Currents Fig. 3.9 shows the optimized current source magnitudes (top row) and corresponding phases (bottom row) for each element, for each optimization method. The optimized phases are fairly constant along the antenna, which indicates that controlling the current distribution magnitude alone provides sufficient degrees of freedom in synthesizing the target temperature profiles. The current magnitudes obtained are realistic, and follow a distribution that can be implemented 26 Chapter 3. The Multiphysics Optimization Design Methodology 150 T ( °C) T ( °C) 150 90 70 100 130 110 90 100 70 50 50 50 2 1 0 -1 -2 5 3 4 50 2 2 1 z (cm) x (cm) 0 -1 -2 3 4 5 z (cm) x (cm) (a) 2 (b) Figure 3.7 Surface plot of T for an xz-slice at y ≈ 0 cm. (a) SAR-based optimization, and (b) T -based optimization. Target region boundary shown in red. Empty region in the centre represents position of source array. -3 -6 -10 -10 -20 -30 2 1 0 -1 -2 x (cm) (a) -3 -6 -10 0 |E| (dB) |E| (dB) 0 5 4 3 -10 -20 2 -30 2 z (cm) 1 0 -1 -2 x (cm) 5 4 3 2 z (cm) (b) Figure 3.8 Surface plot of normalized |E| (in dB) for an xz-slice at y ≈ 0 cm. (a) SAR-based optimization, and (b) T -based optimization. Target region boundary shown in red. Empty region in the centre represents position of source array. in practice, as shown in the following section, where their realization is discussed. In terms of computational cost, the SAR-based method took approximately six hours to converge, while the T -based approach took approximately eight hours. The added cost is justified by the significant advantages provided by the T -based method. 3.4 Physical Implementation of Optimized Current Elements – A Proof of Concept In this section, a possible implementation of the optimized current distributions is presented. This design is a proof of concept demonstrating that ablation antennas can be designed in an inverse way by synthesizing desired current distributions. A more comprehensive antenna design is presented in Chapter 3; the purpose of this section in simply to establish that the 27 Chapter 3. The Multiphysics Optimization Design Methodology 2 1 Surface current magnitude (A) Surface current magnitude (A) 1.2 0.8 0.6 0.4 0.2 0 1.5 1 0.5 0 1 2 3 4 5 6 1 2 Source current element (a) 4 5 6 (b) : Surface current phase (rad) : Surface current phase (rad) 3 Source current element 3:/4 :/2 :/4 0 3:/4 :/2 :/4 0 1 2 3 4 5 6 Source element (c) 1 2 3 4 5 6 Source element (d) Figure 3.9 Top row: Optimized current source element magnitudes. (a) SAR-based optimization, and (b) T -based optimization. Bottom row: Optimized current source element phases. (d) SAR-based optimization, and (e) T -based optimization. optimized current distributions are realistic. Since the optimized surface current magnitudes have a central peak and a smooth decay on each side, we can implement the distributions via a single cylindrical dipole antenna with a varying radius along its axis. Since the optimized phases are nearly constant, we need only deal with the current magnitudes. The cylinder can be divided into six sections, each section corresponding to one of the sources in the optimization procedure described above. Based on the continuity of the normal component of the current density at an interface, a simple method to vary the current magnitude is to modulate the cross-sectional area of the antenna. Hence, the radius of each section is chosen to be proportional to the surface current magnitude required in that section. The antenna was simulated inside the same tissue structure that was used for the optimization procedure, to enable direct comparison of the resulting temperature profiles. The antenna cylinder radii are chosen such that the resulting current in each section corresponds to the T -based optimized currents obtained earlier. The length of each section is equal to the Chapter 3. The Multiphysics Optimization Design Methodology 2.50 mm 1.25 mm 0.41 mm Feed point 1.26 mm 0.44 mm0.09 mm 2 Surface current (A) 28 Physical implementation Target 1.5 1 0.5 0 0 0.2 0.4 0.6 0.8 1 1.2 Distance along probe (cm) Figure 3.10 Proposed proof of concept design and its surface current distribution, compared to target currents. length of each source that was used for optimization (2 mm). Copper was used as the antenna material. The feed point was placed in the middle of the largest section, with a 250 V applied voltage, which provides the correct scaling of currents to match the target currents. The surface currents are sampled at the midpoint of each cylindrical section for comparison with the target currents that we are trying to synthesize. Fig. 3.10 shows the structure that was simulated and the resulting surface current distribution along its length. The distribution is a fairly accurate replication of the optimized currents desired in each section. Fig. 3.11a and Fig. 3.11b show the temperature profile generated by this design. It is clear that the temperature profile is in good agreement with the profile generated via T -based optimization of six sources as in Fig. 3.6b and Fig. 3.7b. Although it is important to design a suitable feeding network for the proposed dipole antenna, this proof of concept demonstrates that current distributions obtained via optimization can be implemented physically in principle, to generate controllable temperature profiles. The presented design is fairly straightforward, reconfigurable (cylinders of different radii can be assembled in any configuration), and resembles commercial antennas while achieving targeted ablation in a more intelligent and controllable way. This design is presented as a straightfor- 29 Chapter 3. The Multiphysics Optimization Design Methodology T ( °C) 150 130 110 90 100 70 50 50 2 1 0 -1 x (cm) (a) -2 5 4 3 2 z (cm) (b) Figure 3.11 (a) Temperature profile generated by the proposed design at selected xz-slices of the 3D domain. Target region boundary shown in red; 50◦ C contour shown in black. (b) Surface plot of T generated by the proposed design for an xz-slice at y ≈ 0 cm. Target region boundary shown in red. Empty region in the centre represents position of source. ward way of achieving desired current distributions on the surface of an antenna. However, there several ways of manipulating source currents, and a more detailed and complete antenna design is presented in the next chapter. 3.5 Summary A multiphysics inverse design methodology for the design of MWA sources was presented in this chapter. The multiphysics approach was compared with the standard approach of considering only the electrical domain without accounting for thermal properties of tissue. The results clearly indicated the advantages of the multiphysics approach in terms of generating more accurate ablation zones based on a given target. The multiphysics approach provides steeper temperature gradients at the target boundary, thus reducing damage to surrounding healthy tissue. A proof-of-concept design consisting of a dipole with varying radii along its axis was presented as a physical implementation of the optimized current distributions yielded by the optimization methodology. A more detailed, printed antenna design is presented in the next chapter. Chapter 4 A Printed Antenna Concept for Microwave Ablation 4.1 Introduction The purpose of this chapter is to design a realistic MWA antenna concept that can be optimized using the multiphysics procedure described above. The antenna is designed such that its current distribution can be controlled and optimized based on a desired target temperature profile. This is achieved by specifying simple design parameters to manipulate the current distribution on the antenna. The antenna is easily configurable to allow ablating targets of different shapes and sizes, and allows the possibility of reconfigurability. This antenna concept has not previously been applied to microwave ablation, to the best of our knowledge. 4.2 Design Considerations There are several constraints within which such an antenna must operate. The relative electrical permittivity in biological tissue commonly lies in the range 40 − 65 . At a frequency of 2.45 GHz this implies that the wavelength in tissue is approximately in the range 1.5 – 2.0 cm. In general, from a clinical perspective, the diameter of an interstitial device such as an ablation antenna or probe of any sort must not exceed 3 mm, for safe percutaneous insertion. Also, since MWA is generally applied to ailments lying closer the skin, and for target regions no smaller than ≈ 1 cm, the total length of the antenna should typically not exceed ≈ 2 cm. Since the antenna must be placed within a biocompatible casing, the entire antenna must fit within a cylinder of diameter 2.5 mm and height 2 cm for it to comply with the above constraints and be comparable to existing clinically acceptable designs. In electrical lengths, this corresponds to a maximum diameter on the order of ≈ λ/7, and maximum height on the order of ≈ 2λ. In addition, another design consideration was to attempt a printed antenna design, as it would be easier to manufacture, potentially cheaper and has not been attempted before in the context of 30 Chapter 4. A Printed Antenna Concept for Microwave Ablation 31 MWA, to the best of our knowledge. Although the optimization of a set of current element magnitudes and phases brings to mind the properties of phased arrays, the space and manufacturing constraints make it difficult to implement phased arrays here. Even if four to six antenna elements were able to fit, it would be very challenging to design a feed network that satisfies these space constraints while allowing sufficient control over individual element excitations without excessive coupling between elements and the feed lines. With the above considerations in mind, it is more reasonable to attempt to manipulate the current distribution on a single antenna with a simple feed, rather than attempt to fit a phased array with a complex feed network. The long but narrow nature of a typical dipole antenna lends itself quite well to the space constraints described above. Also, simple manipulations to the arms and feed of a dipole antenna can allow the current distribution to be manipulated to some extent. Two straightforward ways of achieving this are loading the arms of the dipole reactively, and varying the position of the feed point along the dipole. Thus, in this case, the reactive loading on the arms, and the feed position along the dipole, are the optimization parameters, rather than current magnitudes and phases. Given that the reactive loading of the dipole arms and position of the feed are to be picked based on an optimization procedure for a particular target temperature profile, there is no fixed input impedance to which a feed can be matched. During the optimization procedure, the entire antenna structure must be simulated at each iteration. Thus the widths of the lines that feed the dipole cannot be picked beforehand to have a particular impedance; instead, it is more meaningful to study the input impedance at the port that connects to the entire structure, rather than at the dipole gap. The widths of the feed line are thus chosen based on space constraints, and not their impedance. These geometric characteristics can also be included as optimization variables (in addition to the reactive loading and feed position), however this would greatly increase the computational cost. For the same reasons cited earlier, we focus on designing the antenna such that it is capable of generating symmetric, spherical ablation zones, while allowing the possibility of generating non-symmetric profiles as well. 4.3 Design Description and Optimization The geometry of the proposed antenna is shown in Fig. 4.1. It consists of two back-to-back dipoles printed on either side of a PCB strip. This configuration is necessary such that the structure is symmetric about the plane that cuts through the PCB strip, to enable the generation of symmetric ablation zones. The dipoles are loaded along the arms with lumped capacitors, which can be surface-mounts or printed. The printed arms of the dipole have gaps to allow for this loading. There are two capacitors per arm; this number can be increased or decreased as per specific requirements. Including more capacitors would allow finer control of the current 32 Chapter 4. A Printed Antenna Concept for Microwave Ablation Port 1 Port 2 λ = 18.67 mm C1 C2 C3 1.270 mm 1.500 mm C4 C1 C2 C3 0.381 mm C4 0.127 mm λ = 18.67 mm 0.254 mm 0.127 mm 0.381 mm Ground (a) Copper Shield Teflon Water FR4 (b) Figure 4.1 Proposed antenna design. (a) Perspective view with dimensions and blow-up of metallization layers. (b) Front and side views including outer teflon casing. distribution on the arms, but would increase the number of variables to optimize for, and thus increase computational complexity. The two dipoles are each one wavelength in length, where the wavelength is calculated with respect to the permittivity of the biological tissue being targeted. Although the input impedance of a single-wavelength dipole in infinite, since the arms are to be reactively loaded Chapter 4. A Printed Antenna Concept for Microwave Ablation 33 5.0 cm Liver Ablation Antenna Target 4.5 cm Figure 4.2 Simulation geometry. and the feed position is not fixed at the centre, this limitation does not apply. The PCB is encased in teflon (PTFE) in accordance with existing ablation antenna designs [9, 10], and the space in between the PCB strip and the teflon casing is assumed to be filled with flowing water to maintain the inner temperature at 20◦ C . Most surface mount capacitors available on the market are rated for temperatures up to ≈ 150◦ C, and FR4, which is used as the PCB material, is rated for temperatures up to ≈ 140◦ C; thus, the high temperatures reached during ablation (≈ 100◦ C) are not expected to be of concern. Additionally, since cool flowing water is pumped through the device, the temperatures near the PCB will be significantly lower than 100◦ C. The dipoles are fed by a dual stripline, where the signal lines and ground planes are embedded within the multilayer board, and connected to the dipole arms through vias. The ground planes are tapered gradually (over 3λ/4) to provide a balanced signal to the dipoles. Each dipole and its feed are treated as a separate (but coupled) single-port network. Each port is treated as a lumped port with a prescribed input voltage. It is assumed that the ports are fed by a 50 Ω transmission line. The region to be ablated is assumed to be in human liver, and the simulation domain is a cylinder as shown in Fig. 4.2. The antenna is designed for an operating frequency of 2.45 GHz. The electrical and thermal material properties used for simulating the antenna are given in Table 4.1 [3, 38]. The resulting physical dimensions of the antenna are shown in Fig. 4.1. The width of each arm is chosen to be 1 mm, the feed line width is 0.2 mm, and the total PCB length is 3.25 cm. To control the shape of the generated temperature profile, there are several sets of parameters that can be varied. The values of the capacitors, the positions of the capacitors along the arms, and the feed point along the length of the dipole can all be picked to modulate the current distribution on the arms. These parameters can also be optimized to obtain a desired temperature profile. We first study the effects of varying the capacitance values and the feed positions independently, to establish that the current distribution on the dipole can be controlled to an 34 Chapter 4. A Printed Antenna Concept for Microwave Ablation extent. We then present three optimization case studies to demonstrate the performance of the design. 4.3.1 Effects of Varying Capacitance Values The current through each capacitor should be inversely related to the reactance of that element. Thus, by varying the values of the capacitors, and reactance and thus the current flow can be controlled. Five test cases were simulated. The capacitance and reactance values, as well as the input voltage Vin applied to each port, are outlined in Table 4.2. The two dipoles are identical in each case. For these tests, the capacitors are placed at 1/6, 2/6, 4/6 and 5/6 of the way along the length of the dipoles. In each case, the feed is at the center of each dipole. 4.3.2 Effects of Varying the Feed Position In this case, the lumped capacitor positions are fixed at 0.15d, 0.25d, 0.75d and 0.85d along the dipoles, where d is the total dipole length measured from the end of the dipole closer to the feed. In this case, d = λ. The value of each capacitor is fixed at 65 pF. Three feed positions (as measured from the end of the dipole closer to the feed) were studied: in case 1, the feed is at 0.5d; in case 2, it is at 0.33d; and in case 3, it is at 0.67d. In cases 1 and 2, an input voltage of 10 V was applied to each port. In case 3, an input voltage of 14 V was used. In all cases, the Table 4.1 Material properties used in simulations. Liver Water Teflon FR4 Relative Relative Electrical Thermal permittivity permeability conductivity conductivity εr µr σ (S/m) kth (W/m/K) Heat capacity C (J/kg/K) Density ρ (kg/m3 ) 43.0 80.0 2.1 4.5 3540 4187 1143 1369 1079 1000 1806 1900 1 1 1 1 1.69 5.5 × 10−6 1.0 × 10−25 0.004 0.52 0.59 0.27 0.3 Table 4.2 Input parameters - variable capacitance test cases. Capacitance values are in pF; reactances are in Ω. Case Case Case Case Case 1 2 3 4 5 Vin (V) C1 (X1 ) C2 (X2 ) C3 (X3 ) C4 (X4 ) 14 5 12 4 15 65.0 0.11 65.0 65.0 65.0 1.30 0.13 65.0 65.0 65.0 1.30 65.0 65.0 65.0 0.13 65.0 65.0 65.0 0.11 0.11 (−1) (−600) (−1) (−1) (−1) (−50) (−500) (−1) (−1) (−1) (−50) (−1) (−1) (−1) (−500) (−1) (−1) (−1) (−600) (−600) Chapter 4. A Printed Antenna Concept for Microwave Ablation 35 two dipoles are identical. 4.3.3 Design Optimization for a Test Target In order to demonstrate the optimization of the design, we take three test cases, each consisting of a target with a different shape. The multiphysics temperature-based optimization procedure described earlier is employed here. Each target is a closed surface resembling an ellipsoid in human liver tissue, as depicted in Fig. 4.2. A corresponding target temperature profile, Ttarget , is defined in each case such that the temperature is above 50◦ C within the target region, at 50◦ C at the target boundary, and decays to normal body temperature (37◦ C) elsewhere. The goal is to optimize the capacitor values, input voltage at each port, and the feed position, to generate a temperature profile T that approaches Ttarget . To achieve this, the parameters stated above are optimized to minimize the following objective function: fT = X m,n T m,n − Ttarget 2 (4.1) m,n∈B where B is the set of coordinates corresponding to the target boundary. The capacitor positions are fixed at 0.15d, 0.25d, 0.75d and 0.85d along the dipoles. These positions can also be set as optimization variables, which may lead to more accurate results but a greater computational cost. The capacitance values are constrained between 0.1 pF and 6.5 nF; the input voltage is constrained within 0.1 V and 1000 V; and the feed position is constrained within 0.3d and 0.7d, measured from the end of the dipole closer to the feed. 4.3.4 Simulation Details The forward problems of computing electric fields and temperature profiles were solved using COMSOL Multiphysics. This was interfaced with Matlab for the optimization procedure and for visualization of results. As described previously, in order to reduce computational cost, it is beneficial to solve (2.1) at steady state using a stationary solver. To establish the validity of this approach, the antenna was simulated to compute temperatures in the time domain. The temperatures were probed at four points at varying distances (r) from the antenna, and plotted as a function of time (t) to study the time domain characteristics. The results are shown in Fig. 4.3. It is clear that steady state temperatures were achieved after approximately 7 minutes; since microwave ablation is generally applied for 5 − 10 minutes, it is confirmed that results can be computed at steady state to reduce computational cost. For optimization, the sequential quadratic programming (SQP) algorithm was used for objective minimization. The simulation domain is terminated on all sides with a scattering boundary condition: the COMSOL equivalent of a radiation boundary condition. For each test case, including optimization, the port impedance, input power, delivered power and reflection coefficient are reported. Chapter 4. A Printed Antenna Concept for Microwave Ablation 36 120 r = 0.4 cm r = 0.7 cm r = 1.1 cm r = 2.1 cm T ( °C) 100 80 60 40 2 4 6 8 10 t (min) Figure 4.3 Temperature as a function of time at four probe points. 4.4 Results and Analysis The dipole current distributions, temperature profiles and lumped port parameters for each of the test cases are discussed in this section. The following port parameters are reported: input power, Pin , total power delivered, Pdel , input impedance, Zin and reflection coefficient S11 . Fig. 4.4a and Fig. 4.4b show the current distribution and temperature profiles for each Table 4.3 Port parameters - variable capacitance test cases. Case Case Case Case Case 1 2 3 4 5 Pin (W) Pdel (W) Zin (Ω) S11 (dB) 13.9 10.0 29.7 17.0 17.6 12.0 9.70 28.3 16.8 15.4 42.3 − 36.4j 54.2 − 16.3j 51.9 + 22.9j 48.0 − 9.35j 60.7 − 40.8j −8.51 −16.0 −13.2 −20.3 −8.94 Table 4.4 Port parameters - variable feed position test cases. Case 1 Case 2 Case 3 Pin (W) Pdel (W) Zin (Ω) S11 (dB) 25.6 16.8 15.6 24.6 15.8 13.7 52.8 + 20.6j 70.2 + 21.9j 84.4 + 35.2j −14.1 −12.3 −9.02 37 Chapter 4. A Printed Antenna Concept for Microwave Ablation test case with a varying capacitance. Results are only reported for one of the two ports, as the two dipoles are identical in this case. The correspondence between reactance values, the current peaks and temperature profiles can clearly be seen. For example, it is clear that a higher reactance on the lower arm of each dipole causes the current peak to be pushed towards the upper arm, and causes the temperature to flatten towards the tip of the antenna. Thus, manipulating the arm reactance at various points allows controlling the current and temperature distributions to an extent. Table 4.3 shows the resulting port parameters. The range of Zin values are fairly easy to match, and the input powers required are realistic compared to existing ablation antennas. Despite being unmatched at the ports, the S11 in each case is fairly low. Case 1 Case 2 Case 3 Case 4 Case 5 30 1 Case 1 Case 2 Case 3 Case 4 Case 5 I (A) 0.6 20 z (mm) 0.8 10 0.4 0 0.2 -10 0 -5 -20 0 5 10 15 -20 20 -10 0 10 20 x (mm) z (mm) (a) (b) 30 1 Case 1 Case 2 Case 3 20 z (mm) 0.8 I (A) 0.6 10 0.4 0 0.2 -10 0 -5 Case 1 Case 2 Case 3 -20 0 5 10 z (mm) (c) 15 20 -20 -10 0 10 20 x (mm) (d) Figure 4.4 (a) and (c): Surface current distributions. (b) and (d): Corresponding 50 ◦ C contours at one central slice. (a) and (b): variable capacitance test cases. (c) and (d): variable feed position test cases. Chapter 4. A Printed Antenna Concept for Microwave Ablation 38 Fig. 4.4c and Fig. 4.4d show results for varying feed positions. Results are only reported for one of the two ports, as the two dipoles are identical in this case. It is clear that the current distribution and shape of the temperature profile depend on the feed position. Pushing the feed higher or lower on the dipole allows one to create different distortions in the ellipsoidal profiles. The port parameters are once again realistic in terms of input powers required, and the ports are fairly straightforward to match to an input line. Fig. 4.6 shows the resulting temperature profiles in the three optimization test cases. The plots on the left show selected slices through the 3D domain, while the plots on the right show a corresponding slice in the centre of the domain. The target boundary is represented by red dots; the 50◦ C contour is a black line. It is clear that in each case, the optimized temperature profile very closely approximates the desired temperature profile (i.e. the 50◦ C contour closely follows the target boundary). The three target regions differ in their shapes, and the antenna is able to accurately ablate each type of target. Case 1 Case 2 Case 3 0.8 I (A) 0.6 0.4 0.2 0 0 5 10 15 z (mm) Figure 4.5 Optimized current distributions on one dipole. Fig. 4.5 shows the current distribution on one of the dipoles for each optimization test case. It is clear that the current distribution is successfully manipulated over the three cases to synthesize the required target temperatures, in accordance with the initial goal of this work. Table 4.5 lists the post-optimization port parameters for each port, for each optimization test case: required input power (Pin ), the input impedance (Zin ), reflection coefficient (S11 ), and feed point. The required input powers are realistic from a laboratory or clinical point of view. Despite the fact that the dipoles are not matched to the feed, and the physical dimensions of the design have not been optimized, the S11 is at most −13.3 dB, which indicates a fairly efficient and realistic design. The range of input impedances can quite easily be matched by a matching 39 Chapter 4. A Printed Antenna Concept for Microwave Ablation Table 4.5 Optimized port parameters for each test case. Case I Case II Case III Pin (W) Port 1 Zin (Ω) S11 (dB) 7.33 16.8 14.5 43.3 − 4.33j 59.5 + 16.1j 58.5 + 22.3j −21.4 −15.5 −13.3 Pin (W) Port 2 Zin (Ω) S11 (dB) Feed point 19.8 17.2 15.7 45.2 − 0.83j 60.0 + 15.5j 59.5 + 21.0j −25.7 −15.6 −13.7 0.39d 0.40d 0.40d network placed before the ports. The optimized capacitances lie in the range [1.07, 272] pF, which are physically realizable within the space constraints. 4.5 Summary The concept of multiphysics optimization was applied to the design of a printed dipole antenna for MWA. By appropriately loading the dipole arms with capacitors, and varying the feed position, the current distribution on the antenna can be manipulated and optimized to generate a desired target temperature profile. Simulations of the antenna demonstrate that it is able to accurately ablate targets with realistic power requirements and low return losses. Thus, the antenna is a low-profile and simple design that can be incorporated into a surgical teflon needle, and provides control over its current distribution without the need for a phased array. 40 Chapter 4. A Printed Antenna Concept for Microwave Ablation 90 2 z (cm) 80 70 1 60 0 50 40 -1 30 -2 -1 0 1 2 x (cm) (a) (b) 100 z (cm) 2 80 1 0 60 -1 40 -2 -1 0 1 2 x (cm) (c) (d) 100 2 90 z (cm) 80 1 70 60 0 50 40 -1 30 -2 -1 0 1 2 x (cm) (e) (f ) Figure 4.6 Optimized temperature profiles in ◦ C. Black: 50◦ C contour. Red: target boundary. Chapter 5 Conclusions It has been shown qualitatively and quantitatively that in order to achieve greater accuracy in the treatment of malignant tissue via MWA, it is important to account for the thermal properties and inhomogeneities of the target tissue, and thus the diffusion of deposited heat energy around the region of deposition. The need for a multiphysics approach was established through simulations, and a novel multiphysics optimization technique for the design of MWA probes was presented, and compared with the standard electric field-based approach. The multiphysics methodology was applied to design a novel printed MWA antenna concept. Simulation results of the antenna indicate that it is feasible and advantageous to design ablation antennas in an inverse way to generate controllable, customizable temperature profiles by manipulating the current distribution on the antenna. 5.1 Significance of this Work To the best of our knowledge, the principle of source optimization to generate desired temperature profiles has not previously been applied to MWA. The shapes, sizes, and electrical and thermal properties of regions that are to be ablated can vary significantly from patient to patient; the inverse multiphysics design methodology can be used to design sources that can target ablation zones more accurately while causing less damage to healthy tissue than existing SAR-based methods. The printed dipole MWA antenna presented is a manifestation of the multiphysics methodology. One of the most significant advantages of this design is that it achieves, through just two dipoles, the temperature profile patterning that one would normally obtain using a phased array. The antenna is a low-profile and simple design that can be incorporated into a surgical teflon needle, and provides control over its current distribution without the need for a phased array. Additionally, the optimization of the antenna parameters presented here is significantly simpler than the optimization of tens, hundreds or thousands of elements that is required for HIFU or phased arrays that are used for hyperthermia treatment . Yet this design achieves more accurate, customizable and easily configurable ablation zones than existing MWA designs. 41 Chapter 5. Conclusions 42 In addition, another advantage of this design is that several existing MWA antennas operate at powers as high as 120 W [9, 24]. However, the proposed design can generate comparably large ablation zones at a total power in the range approximately between 30–60 W, which allows for significant power savings. Since this design is printed and configured quite easily, one can envision an “ablation toolkit” consisting of a series of printed MWA antennas, each one configured differently to generate a different canonical temperature profile based on the typical shapes of ablation zones generally required for a certain application. For example, the toolkit could consist of an antenna that generates spherical profiles for ablating tumours, and another antenna that generates elliptical profiles for arrhythmia. The relative ease of manufacturing several antennas on a single PCB, which can then be cut into strips, would allow these toolkits to be more cost-effective than manufacturing an equal number of separate non-printed antennas that currently exist in literature. Additionally, this design methodology can quite easily be integrated into a clinical setting, as described below. 5.1.1 Clinical Integration From a clinical perspective, an ablation antenna with easily configurable current element magnitudes and phases can be incorporated into treatment planning as follows: • A stack of images of the target region can be obtained via MRI, PET, or other imaging techniques, and the different tissues and structures are identified by contouring techniques. • The electrical and thermal properties of the tissues and structures are obtained as documented in literature. • The image stack is taken as the geometry to be used in the optimization procedure described above, with the MWA antenna positioned appropriately inside the target region. • A target temperature profile is created with respect to the image stack, such that ablative temperatures are achieved within the target. • The parameters of the antenna are optimized using the multiphysics techniques described above. If an ablation toolkit consisting of several such printed antennas is available, the antenna configuration that most closely matches the required configuration can be picked to perform the therapy. Ablation antennas designed in this way would allow for more predictable, controllable and customizable ablation of targets while reducing damage to surrounding healthy tissue. Chapter 5. Conclusions 5.1.2 43 Publications Parts of this work have been presented at one conference and been submitted to two journals, currently being written and under review. • S. Sharma and C. D. Sarris, “Novel Antenna Design Concept for Temperature Profile Synthesis in Microwave Ablation,” Physical Review Letters, 2016 (in development). • S. Sharma and C. D. Sarris, “A Novel Multiphysics Optimization-Driven Methodology for the Design of Microwave Ablation Antennas,” IEEE Journal on Multiscale and Multiphysics Computational Techniques, 2016 (under review). • S. Sharma, A. Ludwig and C. D. Sarris, “Design and Optimization of Microwave Ablation Probes for Tumour and Pain Therapy,” 2015 IEEE International Symposium on Antennas and Propagation and North American Radio Science Meeting, July 2015. 5.2 Future Work The work presented above should be seen only as the foundations of an exciting new way to think about the design of MWA antennas. There is much room for improvement in the presented methodology and design, and we believe that there is immense scope for the adaptation of this work in a clinical context. In order to achieve this, there is a significant amount of work yet to be done. This is discussed in the following subsections. 5.2.1 More Efficient Multiphysics Optimization The focus of the multiphysics optimization methodology presented here was not the optimizer itself, but rather the idea of synthesizing a source current distribution based on a target temperature profile. However, the optimization procedure can be improved significantly by making use of impulse response functions of each source element to allow faster field computations. This would require homemade field computation code, as well as a method to couple the field computation with either a commercial solver of the Bioheat equation, or a homemade computation of the same. A finite differences in frequency domain (FDFD) solver was developed for 2D electric field impulse response computation and was validated with comparisons against COMSOL; as well, a 2D Bioheat equation solver was developed to be coupled with the FDFD solver (see Appendix A for details). However, sufficient computational resources were not available to allow solving for the electric field and temperature profiles in 3D in a similar non-iterative way; instead, a commercial solver had to be used, which was unable to store and reuse impulse response functions to speed up computation. Chapter 5. Conclusions 5.2.2 44 Uncertainty Analysis The local tissue geometry realistically consists of uncertainties in its dimensions (given that the geometry is like obtained via an imaging method); as well, the material properties associated with various parts of the tissue consist of uncertainties, as tissue properties are only known to a limited extent. The multiphysics optimization methodology can be enhanced by accounting for these uncertainties when optimizing for the current distributions. 5.2.3 Antenna Reconfigurability The antenna design presented above was proposed with the possibility of reconfigurability in mind. Since the lumped reactive elements on the dipole arms consist of capacitors, it should be possible to instead use variable capacitors (varactors) and bias them appropriately to control their capacitance. The biasing network would consist of DC signals, and thus would not cause interference with the microwave frequency currents through the feed lines. However, the biasing will have to be designed carefully such that it does not distort the local electric fields significantly. However, minor distortions can simply be accounted for in the optimization procedure. This would likely require an additional layer of metallization in the PCB strip. Additionally, reconfigurability may also be possible in the feed postion by designing a switching network. The feed can be connected via switches to several discrete points along the dipole arms, and the actual via that feeds the dipole can be controlled by way of biasing the switches. Although this would be a challenging component to incorporate on an already crowded board, it is an avenue worth exploring for the significant amount of flexibility it would provide. 5.2.4 Further Antenna Design Optimization Several of the design parameters, such as the line widths, the dipole widths, the number of capacitors used, and the positions of the capacitors, were chosen empirically. However, the design can be further optimized by studying the effects of varying these parameters and optimizing them to minimize the return loss. 5.2.5 Fabrication and Testing The permormance of the proposed antenna must be tested experimentally by fabricating it and applying it to a test mass of tissue. Several fabrication challenges need to be overcome to do this: a water pump must be attached to the teflon casing as well as a water reservoir, to allow water to flow in and out of the casing during operation. The PCB must be attached to the teflon casing to keep it in place; this can likely be achieved with the use of epoxy. A matching network must be designed to allow feeding the antenna ports with a standard 50 Ω transmission line. The ablation zones generated by the antenna must be studied to validate the simulation results. Appendix A Finite Differences Approximation of the Bioheat Equation A.1 Electric Field Computation in 2D The electric field can be computed via finite differences in the frequency domain . The advantage of a homemade electric field solver in the context of the optimization methodology presented above, is that one can compute and store the electric field impulse response (i.e. the electric field due to a source current of magnitude 1 A and phase of 0 rad) for each source element. These impulse responses can then be reused algebraically via the principle of superposition, during each iteration of the optimization procedure, rather than having a commercial solver evaluate the fields during each iteration. This approach is described below for the case of a 2D domain. A similar approach can be followed for 3D, given sufficient computational resources. To evaluate the electric field due to a current source in 2D, the equations to be discretized, with current sources, are Maxwell’s equations in the frequency domain, in an anisotropic medium (to account for the use of a perfectly matched layer (PML, ): ∇ × E = −jωµs̄¯H + Jm (A.1a) ∇ × H = jωεs̄¯E + Je (A.1b) where E is the electric field intensity, and H is the magnetic field intensity. Je and Jm are electric and magnetic current sources respectively. Also, sy sz sx s̄¯ = 0 0 0 sx sz sy 0 45 0 0 sx sy sz (A.2) Appendix A. Finite Differences Approximation of the Bioheat Equation 46 where σx jωε σy sy = 1 + jωε σz sz = 1 + jωε sx = 1 + (A.3a) (A.3b) (A.3c) σx , σy and σz are the conductivities in the x̂, ŷ and ẑ directions respectively and ε is the permittivity of the medium. If we now split (A.1a) and (A.1b) up into components, we obtain: ∂Ey ∂Ez − = −jωµsxx Hx + Jm,x ∂y ∂z ∂Ex ∂Ez − = −jωµsyy Hy + Jm,y ∂z ∂x ∂Ey ∂Ex − = −jωµszz Hz + Jm,z ∂x ∂y ∂Hy ∂Hz − = jωεsxx Ex + Je,x ∂y ∂z ∂Hx ∂Hz − = jωεsyy Ey + Je,y ∂z ∂x ∂Hy ∂Hx − = jωεszz Ez + Je,z ∂x ∂y A.1.1 (A.4a) (A.4b) (A.4c) (A.5a) (A.5b) (A.5c) Discretization of Maxwell’s Equations (A.4a) - (A.5c) can be discretized in space using the 2D Yee cell  shown in Fig. A.1. Since the Ez (i, j) Ex , Hy (i, j + Ey , Hx Hz 1 2 , j) 1 2, j (i + (i + Ez 1 2) (i, j + 1) Ey , Hx + 1 2) (i + 12 , j + 1) ŷ Ez (i + 1, j) x̂ Ex , Hy (i + 1, j + 12 ) Ez (i + 1, j + 1) Figure A.1 2D FDTD Yee cell. Appendix A. Finite Differences Approximation of the Bioheat Equation 47 computation will take place in matrix form in Matlab, we can make the following transformation to get integer indices to reference each field point. 1 i→i+ k 2 1 j→j+ k 2 (A.6a) (A.6b) where k = [0, 1, 2 . . .] and the length of k is equal to the number of nodes. This would result in the equivalent Yee cell shown in Fig. 2. Now we can discretize (A.4a) - (A.5c), keeping in Ez (i, j) Ey , Hx (i + 1, j) Ex , Hy (i, j + 1) Hz (i + 1, j + 1) Ez (i, j + 2) Ey , Hx (i + 1, j + 2) ŷ Ez (i + 2, j) x̂ Ex , Hy (i + 2, j + 1) Ez (i + 2, j + 2) Figure A.2 Equivalent 2D FDTD Yee cell. mind that derivatives in z vanish, since we’re working in 2D. Ezi+2,j − Ezi,j i+1,j i+1,j = −jωµi+1,j si+1,j + Jm,x xx Hx ∆y Ezi,j+2 − Ezi,j i,j+1 i,j+1 = −jωµi,j+1 si,j+1 + Jm,y yy Hy ∆x Exi+2,j+1 − Exi,j+1 i+1,j+1 − = −jωµi+1,j+1 si+1,j+1 Hzi+1,j+1 + Jm,z zz ∆y − Eyi+1,j+2 − Eyi+1,j ∆x (A.7a) (A.7b) (A.7c) Hzi+3,j+1 − Hzi+1,j+1 i+2,j+1 = jωεi+2,j+1 si+2,j+1 Exi+2,j+1 + Je,x (A.8a) xx ∆y Hzi+1,j+3 − Hzi+1,j+1 i+1,j+2 = jωεi+1,j+2 si+1,j+2 Eyi+1,j+2 + Je,y (A.8b) yy ∆x Hxi+3,j+2 − Hxi+1,j+2 i+2,j+2 − = jωεi+2,j+2 si+2,j+2 Ezi+2,j+2 + Je,z (A.8c) zz ∆y − Hyi+2,j+3 − Hyi+2,j+1 ∆x The fields above have been defined at grid points shown in Fig. A.2. However, to reduce memory usage and redundancy, we can define separate indices for each of the field components. If the Appendix A. Finite Differences Approximation of the Bioheat Equation 48 total number of Yee cells is (Nx × Ny ), then the field components have the following indices: Ex : i = [1, Nx + 1] j = [1, Ny ] (A.9a) Ey : i = [1, Nx ] j = [1, Ny + 1] (A.9b) Ez : i = [1, Nx + 1] j = [1, Ny + 1] (A.9c) Hx : i = [1, Nx ] j = [1, Ny + 1] (A.9d) Hy : i = [1, Nx + 1] j = [1, Ny ] (A.9e) Hz : i = [1, Nx ] j = [1, Ny ] (A.9f) The grid points for material properties such as permittivities and permeabilities should still be defined with reference to Fig. A.2 and are hereafter referred to with the subscript ‘m’. (A.7a) - (A.8c) now become: Ezi+1,j − Ezi,j im +1,jm m +1,jm Hxi,j + Jm,x = −jωµim +1,jm sixx ∆y Ezi,j+1 − Ezi,j im ,jm +1 m ,jm +1 = −jωµim ,jm +1 siyy Hyi,j + Jm,y ∆x Exi+1,j − Exi,j im +1,jm +1 − = −jωµim +1,jm +1 sizzm +1,jm +1 Hzi,j + Jm,z ∆y − Eyi,j+1 − Eyi,j ∆x (A.10a) (A.10b) (A.10c) Hzi+1,j − Hzi,j im +2,jm +1 m +2,jm +1 = jωεim +2,jm +1 sixx Exi+1,j + Je,x ∆y (A.11a) Hzi,j+1 − Hzi,j im +1,jm +2 m +1,jm +2 = jωεim +1,jm +2 siyy Eyi,j+1 + Je,y − ∆x (A.11b) Hyi+1,j+1 − Hyi+1,j Hxi+1,j+1 − Hxi,j+1 im +2,jm +2 − = jωεim +2,jm +2 sizzm +2,jm +2 Ezi+1,j+1 + Je,z ∆x ∆y (A.11c) ∂Je,x ∂Je,y 1 1 2 2 Lx,2D + 2 Ly,2D + γ I Hz = ∂y − ∂x (∆x) (∆y) where 2 2 γ = ω εµ Lx,2D sx sy sz (A.12) 2 L 0 0 ... 0 0 0 L 0 . . . 0 0 = . .. .. .. . . 0 0 0 ... 0 L (A.13) (A.14) Appendix A. Finite Differences Approximation of the Bioheat Equation 49 where 0 is a nx × nx matrix of zeros, and the size of Lx,2D is nx ny × nx ny . Also, L is a matrix of size nx × nx defined as: −2 1 0 0 ... 0 . . . 1 −2 1 1 −2 1 . . . 0 L= . .. . . . 0 1 −2 1 ... 0 0 1 −2 (A.15) Also, Ly,2D −2I I 0 ... −2I I I .. .. = . . 0 0 0 0 0 0 ... 0 0 0 I −2I I 0 I −2I 0 0 0 .. . (A.16) where I is an nx × nx identity matrix and there are ny matrices along the diagonal. Now, if we have a y-directed source Jy defined at a point (i, j), we can treat it as two sources for Hz as in Inan, Marshall (2011) page 350 , and write the equation for Hz as follows: 1 1 2 Lx,2D + Ly,2D + γ I Hz = − (∇ × J)z (∆x)2 (∆y)2 (A.17) with (∇ × J)z |i,j+1 = − Jy |i,j 2∆x (∇ × J)z |i,j = 0 (∇ × J)z |i,j−1 = Jy |i,j 2∆x (A.18a) (A.18b) (A.18c) (A.18d) Now we can calculate Ex and Ey using Maxwell’s equations: Ex = 1 ∂Hz jωεS̄¯(1, 1) ∂y Ey = − 1 ∂Hz ¯ jωεS̄ (2, 2) ∂x (A.19a) (A.19b) 50 Appendix A. Finite Differences Approximation of the Bioheat Equation A.2 Computation of the Bioheat Equation in 2D Given the electric fields computed above, the temperature profile is computed by discretizing the Bioheat equation : ρCp ∂T σ = |E|2 + ∇ · (k∇T ) + Wb Cb (Tb − T ) ∂t 2 (A.20) where Cp is the heat capacity per unit mass, ρ is the mass density, k is the thermal conductivity, Wb is the mass flow density for blood, Cb is the specific heat for blood, and Tb is the normal blood temperature, 37◦ C. In the frequency domain, equation 1 becomes: jωρCp T (ω) = σ |E(ω)|2 + ∇ · (k∇T (ω)) + Wb Cb (Tb − T (ω)) 2 (A.21) This can be discretized in 2D space as done before. First we define the following variables: Ti,j+1 − 2Ti,j + Ti,j−1 ∆x2 Ti+1,j − 2Ti,j + Ti−1,j Ỹ = ∆y 2 X̃ = Z̃ = 0 (A.22a) (A.22b) (A.22c) Then, −k X̃ + Ỹ + (jωρCp + Wb Cb ) Ti,j = 0.5σ (Ei,j )2 + Wb Cb Tb (A.23) Now, following the same discretization approach as in Inan, Marshall (2011), page 349 , we can write this as: −k −k 2 2 Lx,2D + 2 Ly,2D + γI T = 0.5σ (E) + Wb Cb Tb (∆x) (∆y) (A.24) where T and E are vectorized, and γ = jωρCp + Wb Cb (A.25) Now T can be calculated via matrix inversion. A.3 Validation of FDFD Results The electric field computation via FDFD was compared with COMSOL to confirm the correct functioning of the PML. The z component of the electric field, Ez , computed by the homemade FDFD code and by COMSOL are shown respectively in Fig. A.3a and Fig. A.3b. For a better comparison, Fig. A.3c and Fig. A.3d show Ez along the lines drawn in Fig. A.3a and Fig. A.3b. Fig. A.3c shows Ez along the vertical lines. Fig. A.3d shows Ez along the horizontal lines. It is clear that the finite differences solver computes the electric fields accurately, and if implemented Appendix A. Finite Differences Approximation of the Bioheat Equation 51 in 3D, it is a suitable method for calculating the electric field impulse responses for each source element, and can be stored and used in the optimization procedure to speed up the process. (a) (b) (c) (d) Figure A.3 (a) Ez computed via homemade FDFD. (b) Ez computed in COMSOL. (c) Ez along the vertical test line shown in (a) and (b). Blue line: computed in COMSOL. Red crosses: computed via FDFD. (d) Ez along the horizontal test line shown in (a) and (b). Blue line: computed in COMSOL. Red crosses: computed via FDFD. Ez values are in V/m. Appendix B Material Property Estimation by Inverse Scattering B.1 Introduction In this section, a separate but related project is discussed. It is often necessary to determine the electrical, thermal or other material properties of an unknown object or domain. This is particularly necessary for the optimization methodology discussed in this work; it is important to estimate as accurately as possible the electrical and thermal properties of the target region in tissue in order to optimize the MWA antenna. This is particularly difficult if the electrical and thermal discontinuities are to be accounted for, which is important for accurate targeting, as shown in Chapter 2. A method to estimate and reconstruct the electrical properties of an unknown 2D domain was presented by Rekanos in 2002  (subsequently published in further detail in 2012 ), and is presented and discussed here in the context of MWA. If a set of electromagnetic transmitters and receivers is placed outside the region of the interest to reconstruct, the scattered fields can be analyzed to reconstruct the electrical properties of the target region. Rekanos presents a finite-differences in time-domain algorithm to perform this reconstruction, and the algorithm, referred to as an inverse scattering algorithm, is tested via simulations. We implement the algorithm in Matlab. B.2 The Inverse Scattering Algorithm We assume that the unknown domain of interest is in 2D, and is surrounded by eight transmitter/receiver pairs evenly spaced around the target domain in a circle. This setup is illustrated in Fig. B.1. The entire domain is divided into 70 × 70 Yee cells for the FDTD procedure , and electromagnetic fields are evaluated in the entire domain. This includes an outer padding of 6 PML cells  along the boundary. The target domain, in the centre of the whole domain, 52 53 Appendix B. Material Property Estimation by Inverse Scattering consists of 26 × 26 cells. The size of each cell is 6.25 mm. Each transmitter/receiver is a point current source in the z direction, which is normal to the plane of the page, facing outwards. 70 60 50 40 30 20 10 0 0 20 40 60 80 Figure B.1 Geometry for 2D FDTD-based inverse scattering algorithm. Small black squares: transmitter/receiver pairs. Yellow shaded square: unknown domain being imaged. The scattered fields are sampled at the receiver points. The transmitters are iteratively excited one-by-one, and the scattered fields are sampled at each of the other seven receiver points, other than the one that corresponds to the transmitter for a particular iteration. An objective functional is defined based on the differences between the measured field and the calculated field based on the current estimate of the properties of the unknown domain. This functional is minimized using a conjugate gradients algorithm. The algorithm is presented in detail in  and ; the purpose of this section is to outline the computational differences in our implementation. In our implementation, each transmitter provides a modulated Gaussian excitation as follows: Jz = exp − (t − t0 )2 t2w ! sin (ωt) (B.1) where Jz is the source current density, t0 is the time point corresponding to the peak of the modulated Gaussian, tw is related to the temporal standard deviation of the pulse, and ω is the angular frequency of the modulating sinusoid. In this case, the frequency is chosen to be 2.31 GHz. 301 time steps are used, ranging from t = 0 to t = 2.17 ns, with a time step of Appendix B. Material Property Estimation by Inverse Scattering 54 7.22 ps. This excitation is chosen such that it is temporally and spatially smooth enough to prevent numerical errors in the FDTD implementation. The electric fields generated by each source are calculated via homemade FDTD code based on , accounting for the PML based on . For the line-search minimization step in the optimization algorithm described in  and , we used Matlab’s patternsearch function to evaluate the step length that minimizes the described functional for each iteration. The optimization procedure attempts to image the permittivity, permeability and conductivity of the target region. For applications in MWA where targeted tissues are generally non-magnetic, it is sufficient to image only for the permittivity and conductivity. B.3 Results To test the algorithm, we assume that the target region contains two objects, with a relative permittivity profile as shown in Fig. B.2a, and an electrical conductivity profile as shown in Fig. B.3a. The optimization procedure was run for 40 iterations. The resulting imaged relative permittivity profile is shown in Fig. B.2b; the imaged conductivity profile is shown in Fig. B.3b. It is clear that the imaged profiles are fairly good approximations of the targets. The infinity norm (maximum deviation) in relative permittivity is 0.59; the infinity norm in electrical conductivity is 0.12 S/m. (a) (b) Figure B.2 (a) Actual electrical permittivity profile to be imaged. (b) Imaged permittivity profile obtained via inverse scattering algorithm. Actual profile in semi-transparent red. B.4 Relevance to Microwave Ablation Given an algorithm that is able to reconstruct the electrical properties of a target domain, a set of transmitter/receiver pairs can be used to image these properties in the target region to be ablated. Knowledge of these properties is necessary in order to implement the MWA Appendix B. Material Property Estimation by Inverse Scattering (a) 55 (b) Figure B.3 (a) Actual electrical conductivity profile to be imaged. (b) Imaged conductivity profile obtained via inverse scattering algorithm. Actual profile in semi-transparent red. antenna optimization methodology described in this thesis. In addition, there are further future directions that must be taken in order to make full use of the inverse scattering algorithm: • The algorithm must be extended to 3D in order to allow full compatibility with 3D ablation targets. • An equivalent algorithm must be developed to reconstruct the thermal properties of the target region to enable the use of the multiphysics optimization methodology. • Uncertainty quantification should be incorporated into the inverse scattering procedure to account for the finite resolution of medical imaging methods, and allow for more accurate targeting in designing MWA antennas. Bibliography  C. L. Brace, “Radiofrequency and microwave ablation of the liver, lung, kidney, and bone: What are the differences?” Current Problems in Diagnostic Radiology, vol. 38, no. 3, pp. 135–143, May 2009.  H.-M. Chiu, A. S. Mohan, A. R. Weily, D. J. R. Guy, and D. L. Ross, “Analysis of a novel expanded tip wire (etw) antenna for microwave ablation of cardiac arrhythmias,” IEEE Transactions on Biomedical Engineering, vol. 50, no. 7, pp. 890–899, July 2003.  E. R. C. Jr. and E. R. C. Sr., “Electric and thermal field effects in tissue around radiofrequency electrodes,” Pain Medicine, vol. 6, no. 6, pp. 405–424, December 2005.  A. Karampatzakis, S. Khn, G. Tsanidis, E. Neufeld, T. Samaras, and N. Kuster, “Antenna design and tissue parameters considerations for an improved modelling of microwave ablation in the liver,” Physics in Medicine and Biology, vol. 58, no. 10, pp. 3191–3206, May 2013.  N. C. Yu, D. S. K. Lu, S. S. Raman, D. E. Dupuy, C. J. Simon, C. Lassman, B. I. Aswad, D. Ianniti, , and R. W. Busuttil, “Hepatocellular carcinoma: Microwave ablation with multiple straight and loop antenna clusters – pilot comparison with pathologic findings,” Radiology, vol. 239, no. 1, pp. 269–275, April 2006.  G. Schaller, J. Erb, and R. Engelbrecht, “Field simulation of dipole antennas for interstitial microwave hyperthermia,” IEEE Transactions on Microwave Theory and Techniques, vol. 44, no. 6, pp. 887–895, June 1996.  C. L. Brace, P. F. Laeseke, D. W. van der Weide, and F. T. L. Jr., “Microwave ablation with a triaxial antenna: Results in ex vivo bovine liver,” IEEE Transactions on Microwave Theory and Techniques, vol. 53, no. 1, pp. 215–220, January 2005.  Z. Ibitoyea, E. Nwoyeb, M. Awedaa, A. Oremosuc, C. Annunobid, and O. Akanmu, “Optimization of dual slot antenna using floating metallic sleeve for microwave ablation,” Medical Engineering and Physics, vol. 37, no. 4, pp. 384–391, April 2015.  D. Yang, J. M. Bertram, M. C. Converse, A. P. ORourke, J. G. Webster, S. C. Hagness, J. A. Will, and D. M. Mahvi, “A floating sleeve antenna yields localized hepatic microwave 56 Bibliography 57 ablation,” IEEE Transactions on Biomedical Engineering, vol. 53, no. 3, pp. 533–537, March 2006.  J. M. Bertram, D. Yang, M. C. Converse, J. G. Webster, and D. M. Mahvi, “A review of coaxial-based interstitial antennas for hepatic microwave ablation,” Critical Reviews in Biomedical Engineering, vol. 34, no. 3, pp. 187–213, 2006.  K. Tabuse, “A new operative procedure of hepatic surgery using a microwave tissue coagulator,” Nihon Geka Hokan, vol. 48, no. 2, pp. 160–172, March 1979.  K. Tabuse, M. Katsumi, Y. Kobayashi, H. Noguchi, H. Egawa, O. Aoyama, H. Kim, Y. Nagai, H. Yamaue, K. Mori, Y. Azuma, and T. Tsuji, “Microwave surgery: Hepatectomy using a microwave tissue coagulator,” World Journal of Surgery, vol. 9, no. 1, pp. 136–143, February 1985.  P. Liang and Y. Wang, “Microwave ablation of hepatocellular carcinoma,” Oncology, vol. 71, no. 1, pp. 124–131, December 2007.  S. A. Shock, K. Meredith, T. F. Warner, L. A. Sampson, A. S. Wright, T. C. W. III, D. M. Mahvi, J. P. Fine, and F. T. Lee, “Microwave ablation with loop antenna: In vivo porcine liver model,” Radiology, vol. 231, no. 1, pp. 143–149, April 2004.  J. G. Lynn, R. L. Zwemer, A. J. Chick, and A. E. Miller, “A new method for the generation and use of focused ultrasound in experimental biology,” Journal of General Physiology, vol. 26, no. 2, pp. 179–193, November 1942.  H. Yu and C. T. Burke, “Comparison of percutaneous ablation technologies in the treatment of malignant liver tumors,” Seminars in Interventional Radiology, vol. 31, no. 2, pp. 129–137, June 2014.  W.-P. Zhao, Z.-Y. Han, J. Zhang, and P. Liang, “A retrospective comparison of microwave ablation and high intensity focused ultrasound for treating symptomatic uterine fibroids,” Radiology, vol. 84, no. 3, pp. 413–417, March 2015.  D. Haemmerich, S. Tungjitkusolmun, S. T. Staelin, J. Fred T. Lee, D. M. Mahvi, and J. G. Webster, “Finite-element analysis of hepatic multiple probe radio-frequency ablation,” IEEE Transactions on Biomedical Engineering, vol. 49, no. 7, pp. 836–842, July 2002.  C. L. Brace, “Microwave tissue ablation: Biophysics, technology, and applications,” Critical Reviews in Biomedical Engineering, vol. 38, no. 1, pp. 65–78, 2010.  M. G. Lubner, C. L. Brace, J. L. Hinshaw, and F. T. L. Jr., “Microwave tumor ablation: Mechanism of action, clinical results and devices,” Journal of Vascular and Interventional Radiology, vol. 21, no. 8 Suppl, pp. S192–203, August 2010. Bibliography 58  A. S. Wright, L. A. Sampson, T. F. Warner, D. M. Mahvi, and F. T. Lee, “Radiofrequency versus microwave ablation in a hepatic porcine models,” Radiology, vol. 236, no. 1, pp. 132–139, July 2005.  C. L. Brace, “Microwave ablation technology: What every user should know,” Current Problems in Diagnostic Radiology, vol. 38, no. 2, pp. 61–67, April 2009.  S. Labonte, A. Blais, S. R. Legault, H. O. Ali, and L. Roy, “Monopole antennas catheter ablation,” IEEE Transactions on Microwave Theory and Techniques, vol. 44, no. 10, pp. 1832–1840, October 1996.  P. Prakash, M. C. Converse, J. G. Webster, and D. M. Mahvi, “An optimal sliding choke antenna for hepatic microwave ablation,” IEEE Transactions on Biomedical Engineering, vol. 56, no. 10, pp. 2470–2476, October 2009.  K. M. Jones, J. A. Mechling, J. W. Strohbehn, and B. S. Trembly, “Theoretical and experimental sar distributions for interstitial dipole antenna arrays used in hyperthermia,” IEEE Transactions on Microwave Theory and Techniques, vol. 37, no. 8, pp. 1200–1209, August 1989.  W. Hurter, F. Reinbold, and W. J. Lorenz, “A dipole antenna for interstitial microwave hyperthermia,” IEEE Transactions on Microwave Theory and Techniques, vol. 39, no. 6, pp. 1048–1054, June 1991.  P. Prakash, G. Deng, M. C. Converse, J. G. Webster, D. M. Mahvi, and M. C. Ferris, “Design optimization of a robust sleeve antenna for hepatic microwave ablation,” Physics in Medicine and Biology, vol. 53, pp. 1057–1069, January 2008.  H. Luyen, S. C. Hagness, and N. Behdad, “A balun-free helical antenna for minimally invasive microwave ablation,” IEEE Transactions on Antennas and Propagation, vol. 63, no. 3, pp. 959–965, March 2015.  D. T. Inc., “Diros OWL Trident Cannulae,” http://dirostech.com/new-innovative/#!, [Online; accessed 10-August-2016].  H. H. Pennes, “Analysis of tissue and arterial blood temperatures in the resting human forearm,” Journal of Applied Physiology, vol. 1, no. 2, pp. 93–122, August 1948.  M. Ahmed, Z. Liu, S. Humphries, and S. N. Goldberg, “Computer modeling of the combined effects of perfusion, electrical conductivity, and thermal conductivity on tissue heating patterns in radiofrequency tumor ablation,” International Journal of Hyperthermia, vol. 24, no. 7, pp. 577–588, November 2008.  C. L. Dolph, “A current distribution for broadside arrays which optimizes the relationship between beamwidth and side-lobe level,” Proceedings of the IRE, vol. 34, no. 6, pp. 335–348, June 1946. Bibliography 59  D. Barbiere, “A method for calculating the current distribution of tschebyscheff arrays,” Proceedings of the IRE, vol. 40, no. 1, pp. 78–82, January 1952.  V. Ekstrand, H. Wiksell, I. Schultz, B. Sandstedt, S. Rotstein, and A. Eriksson, “Influence of electrical and thermal properties on rf ablation of breast cancer: is the tumour preferentially heated?” BioMedical Engineering OnLine, July 2005.  M. Ahmeda, Z. Liua, S. Humphriesb, and S. N. Goldberga, “Computer modeling of the combined effects of perfusion, electrical conductivity, and thermal conductivity on tissue heating patterns in radiofrequency tumor ablation,” International Journal of Hyperthermia, vol. 24, no. 7, pp. 577–588, November 2008.  S. K. Das, S. T. Clegg, and T. V. Samulski, “Computational techniques for fast hyperthermia temperature optimization,” Medical Physics, vol. 26, no. 2, pp. 319–328, February 1999.  A. Ludwig, C. D. Sarris, and G. V. Eleftheriades, “Near-field antenna arrays for steerable sub-wavelength magnetic-field beams,” IEEE Transactions on Antennas and Propagation, vol. 62, no. 7, pp. 3543–3556, July 2014.  P. Hasgall, E. Neufeld, M.-C. Gosselin, A. Klingenbck, and N. Kuster, “It’is database for thermal and electromagnetic parameters of biological tissues,” September 2015.  J.-P. Berenger, “A perfectly matched layer for the absorption of electromagnetic waves,” Journal of Computational Physics, vol. 114, no. 2, pp. 185–200, October 1994.  A. P. O’Rourke, M. Lazebnik, J. M. Bertram, M. C. Converse, S. C. Hagness, J. G. Webster, and D. M. Mahvi, “Dielectric properties of human normal, malignant and cirrhotic liver tissue: in vivo and ex vivo measurements from 0.5 to 20 ghz using a precision open-ended coaxial probe,” Physics in Medicine and Biology, vol. 52, no. 15, pp. 4707–4719, August 2007.  B. T. McWilliams, E. E. Schnell, S. Curto, T. M. Fahrbach, and P. Prakash, “A directional interstitial antenna for microwave tissue ablation: Theoretical and experimental investigation,” IEEE Transactions on Biomedical Engineering, vol. 62, no. 9, pp. 2144–2150, September 2015.  U. S. Inan and R. A. Marshall, Numerical Electromagnetics: the FDTD Method, 1st ed. Cambridge University Press, 2014.  I. T. Rekanos, “Inverse scattering in the time domain: An iterative method using an FDTD sensitivity analysis scheme,” IEEE Transactions on Magnetics, vol. 38, no. 2, pp. 1117–1120, March 2002. 60 Bibliography  ——, “Time-domain inverse scattering using lagrange multipliers: An iterative FDTDbased optimization technique,” Journal of Electromagnetic Waves and Applications, vol. 17, no. 2, pp. 271–289, April 2012.  A. Taflove and S. C. Hagness, Computational Electrodynamics: The Finite-Difference Time-Domain Method, 3rd ed. Artech House Publishers, 2005.  S. D. Gedney, “An anisotropic perfectly matched layer-absorbing medium for the truncation of FDTD lattices,” IEEE Transactions on Antennas and Propagation, vol. 44, no. 12, pp. 1630–1639, December 1996.