1246

IEEE TRANSACTIONS ON BIOMEDICAL ENGINEERING. VOL. 38, NO. 12. DECEMBER 1991

Numerical Optimization of 3-D SAR Distributions in Cylindrical Models for Electromagnetic Hyperthermia Dipakbin Q. Chowdhury, Student Member, IEEE, and Steven C . Hill, Member, IEEE

Abstract-Numerically optimized SAR (specific absorption rate) distributions in a source free 3-D multilayered concentric cylindrical (MCC) model are presented. The fields were expanded in the modes of the MCC. Cost functions which specify mathematically the relative weight assigned to differences between an SAR distribution and a desired SAR distribution were

close a particular SAR distribution is to the desired distribution (an example of a simple cost function is the velOf the 'quared difference between the SAR and the desired SAR). Once such an ideal SAR distribution is computed it can be compared with the optimal SAR

SAR is largest. A five-layered model, including the outer water layer for cooling and improved matching with the source, was used. The frequency was 70 MHz. The current and charge distributions computed on a perfectly conducting cylindrical surface just outside the model are also shown. The surface current and charge distributions depend strongly on the relative importance of the cost for acute heat and systemic heat. A technique is developed for generating a new set of basis functions for reducing the number of unknowns to be optimized. We suggest that the approach shown could be useful in designing hyperthermia applicators.

tors. There are a number of previous reports of optimization of numerically calculated SAR or temperature distributions for EM hyperthermia. For example, Archangeli et al. [4] optimized the SAR inside a 2-D inhomogeneous model of the human chest. Strohbehn et al. [5] optimized the SAR distribution generated by an annular phased array in a 2-D discrete inhomogeneous numerical model. Rappaport and Morgenthaler optimized power at the center of a lossy homogeneous sphere 161, and optimized the phase of axially polarized sources in order to focus power at the center of a homogeneous cylinder [7]. Morita et al. [8] optimized the SAR generated by a set of line sources in a 2-D inhomogeneous elliptical model. De Wagter [9] optimized the transient temperature distribution generated by multiple EM applicators in a numerically discretized 2-D inhomogeneous model. Knudsen and Hartmann [ 101 optimized the temperature distribution in a homogeneous 2-D cylindrical model. In this paper our objective is to determine the best achievable, nonaxisymmetric 3-D SAR distributions for one type of inhomogeneous body model, a multilayered concentric cylindrical model (MCC). Two types of basis functions are employed: 1) the TE and TM eigenmodes of the MCC system, and 2) a set of basis functions which are combinations of the TE and TM modes. Expansion of the fields in the lower order eigenmodes of the MCC provides the maximum variation in the field. The coefficients of the fields are determined using gradient search optimization which minimizes the differences between the desired and achievable SAR. The cost function used can emphasize either the systemic or the acute heating. The desired SAR distribution is one which we assume could have been determined, using bioheat models, to be best for heating and would not necessarily overlap the region to

I. INTRODUCTION

T

HERE are fundamental limits on the extent to which electromagnetic energy can be localized in the noninvasive generation of hyperthermia for cancer therapy [1]-[3]. Energy deposited into a part of the body to be heated is always accompanied by energy deposition into other regions, and this problem is particularly apparent with noninvasive hyperthermia. A potentially useful concept in the evaluation of hyperthermia applicators is that of an ideal SAR distribution. We define an ideal specific absorption rate (SAR) distribution as the distribution, consistent with Maxwell's equations, that is closest to a particular desired SAR distribution. The definition has meaning only if the following are specified: a body or a model of the body, a desired SAR distribution, and a cost function for assessing how Manuscript received April 1, 1991; revised August 30, 1991. This work was supported by the National Science Foundation Grant EET-8708687. D. Q. Chowdhury is with the Department of Applied Physics and Center for Laser Diagnostics, Yale University, New Haven, CT 06520. S . C. Hill is with the Department of Electrical and Computer Engineering, New Mexico State University, Las Cruces, NM 88003, and Atmospheric Sciences Laboratory, White Sands Missile Range, NM 88002. IEEE Log Number 9103957.

0018-9294/91$01.00

1

0 1991 IEEE

CHOWDHURY AND HILL: NUMERICAL OPTIMIZATION OF 3-D SAR DISTRIBUTIONS

be heated. In a patient the desired SAR distribution would need to be determined using thermal models. In this paper we model only the EM problem. Finite element [ l l ] or finite difference time domain [ 121 methods would be more appropriate for modeling more realistic body models. This paper first presents the method of computation of the fields in MCC models, and then presents the cost function and basis functions used in the optimization. Then the results of optimizing the SAR distributions and the current and charge distribution associated with the optimal SAR distributions are presented and discussed. Then limitations of the approach and work presented are discussed. OF FIELDSI N MULTILAYERED 11. COMPUTATION CONCENTRIC CYLINDERS The fields in a multilayered cylindrical model have been described previously [ 131, [ 141. In each layer the TM and the TE potentials are expressed as

1247

0.005%. When k, is large the fields near the surface are much larger than they are further into the model. Consequently, the contribution of high-k, fields to a desired SAR distribution is small if not negligible unless the region where power is to be deposited is near the surface. As n increases the acceptable kz decreases, and vice versa. Our objective is to optimize the fields and SAR inside the cyclinder in a region of length L. The external fields are ignored during the optimization. Since we try to minimize the SAR near the ends of the cylinder, the currents leaving the ends of the cylinder are kept relatively small. The SAR is defined as: (1 /2) ( o l d ) IEI2 where U is the conductivity and d is the density of the medium. Since the SAR is proportional to EE", a combination of orders n is required to obtain a &dependent SAR. Similarly, multiple kz's are required for a z-dependent SAR. The TE and TM fields are coupled for the distributions considered here since both n and k, must take nonzero values. The total potential is ti=m

$?&I

=

nm

[A?,(~,)J,,(~,,P)+ ~ ? , ( k , ~) , , ( k , / ~ ) ~ e ~ ~ ' e ~ ' - ~ (1)

$:,(k,) = [Af,,,(k)J,(k,/P) + B?.,(k,) yn(k,/p)leJn"OJk-:

(2) where J, (k,,p) and Y,, ( k , / p ) are cylindrical Bessel and Neumann functions of order n ; A;", B;", A;, and B ; are expansion coefficients for the TE and the TM fields in the Ith layer; k, is the wave propagation constant in the z direction, and k,/ is the radial propagation constant of the Ith layer. These propagation constants are related to the medium propagation constants by k:

= w2p1eI =

ki

+ k:.

(3)

where p, is the permeability and e, is the complex permittivity in layer 1. An eJw' time variation is assumed. In the central region (I = 1 ) the coefficients of the Y , , ( k , p ) are zero. The electric fields are obtained from the potentials as [I51

v

x (i$;(k,)),

(4) where i is the unit vector in the z direction. The relations between the A, and B, for different layers are obtained by matching the tangential components of the electric and magnetic fields ( E z , H,, E,, H , ) at the boundaries between layers. The Bessel and Neumann functions were computed using the method and computer code described by Mason [ 161. To avoid the effects of round-off errors that occurred in computing the matrices and fields for large kZ and n , the boundary conditions for the tangential E and normal displacement vector were checked at each layer for each pair of k, and n. The fields were not included whenever the boundary conditions were not matched to within -

B

where $ and $ are defined in (1) and (2), a, ( k , ) are complex coefficients for the modes that are TM in the core, and b,(k,) are complex coefficients for the modes that are TE in the core. The n's take on only integer values. We assume the k, are real. Hence the fields are valid for an infinite cylinder or for a finite length cylinder with no sources at the ends. The E fields obtained [15] from $( p , z , 4 ) of ( 5 ) are

where ETMis the electric field of a mode that is TM in the core and ETEis the electric field of a mode that is TE in the core, and the integral in ( 5 ) has been approximated as a discrete sum.

111. OPTIMIZATION

A. Basis Functions Two types of basis functions are employed: first, the eigenmodes described by (1) and (2), and second, linear combinations of these modes which confine the fields axially, and usually azimuthally and radially. The eigenmodes would form a complete set of basis functions if all were used. However, only a few n and k, are used for the optimizations employing eigenmodes. Modes having high k, and/or n deposit most of their power close to the outer surface [13] and would contribute little to a best SAR distribution for the models considered here. Computations (not shown) indicate that a maximum k, of 15 cm-I, a maximum n of 4, and a Ak, of 5 cm-l (a total

1248

IEEE TRANSACTIONS ON BIOMEDICAL ENGINEERING. VOL. 38. NO. 12. DECEMBER 1991

of 20 TE and 20 TM modes) are sufficient to express the best SAR distributions considered here to within a few percent. If we were to try to deposit the maximum power in a small region nearer to the surface then larger k; and n would be required. The generation and use of the second set of basis functions is motivated by a desire to reduce the number of unknowns required for optimization. These combined basis functions are obtained by adding the eigenmode basis functions in such a way that the total E field is localized, usually in all three directions. Axial localization is achieved by combining modes having different k,’s. Similarly, azimuthal localization is achieved by combining modes with different I Z ’ S . The E-field amplitude distributions in the axial and the azimuthal directions are chosen to be Gaussian on a cylindrical surface outside the model, i.e., the amplitudes are assumed to vary as e-(z’2‘)2 e where 2, and +tz are constants which determine the widths of the distributions. Other distributions of amplitude would require broader range of 11’s and kZ’s [13]. Radial confinement is achieved primarily by choosing the phases of the modes so that the fields add constructively at the region of interest. The Gaussian E-field amplitude distributions are chosen to help confine the fields in the axial and azimuthal directions. The phases of the fields on the surface are estimated (assuming rectilinear propagation except at the interfaces) so that the associated fields add in phase at the desired location inside the cylinder. To estimate the phase information at the surface outside the model we initiate rays from the center of the region to be heated and let them propagate out to the source surface. Then the total phase shift that each ray undergoes on its way to the surface is computed. The rays obeys Snell’s law of refraction at each interface. Diffraction of the rays is neglected. The ray picture of EM waves is applicable only when the dimensions involved are much larger than the wavelength and so it is not as accurate as desired in this application where the objects are smaller than the wavelength. Here the ray method is used only to estimate the phase. To estimate the phase of the E field over the whole surface, the rays are initiated in enough different directions in z and to cover the required region on the surface. Rays are totally reflected for some incident angles at an interface (those exceeding the critical angle) where the inner layer has higher permittivity than the outer layer. As B result it is not possible to trace the rays from the center of the region to be heated to all points on the outer surface. Phases are assigned to such points with an aim to avoid any abrupt change in phase. That is, a smoothing algorithm is used to keep the Fourier spectrum relatively narrow. A similar approach for determining the phase of axially polarized annular sources was used [7] to focus power on the axis at the center of a homogeneous cylinder. Only two of the three field components at the surface can be set to a desired distribution of phase and ampli-

tude. In the work presented here the and Ez1,,,,are specified as functions of and z at a certain radial distance ( p o ) . Then Fourier transformation is used to determine the an,k,and btl,k;.The transformation is written:

+

E ~( po, ~ Z , ~4) ,= ~C ~ Ci,kreJ(nd ~ +‘A n.k:

(7) where C i , k and Ct,L.are the Fourier coefficients for transand E,lo,, respectively. From (6) and (7) forms of Ez,c,,a, the relation between Fourier coefficients and an.k.,bn,k is

ci,L.= a n , k

(PO)

+ bn,k:E:,E~

( PO)

So, once E; and E, are specified at a radial distance p O , the C&, and C!,k:are obtained from (7), then the an,k, and b,,+ are found using (8) and finally E,, the E field of the ith basis function, can be computed with (6). The computed SAR is determined from the total E field which is the weighted sum of all the basis functions, i.e., R p , 2,

4)

=

ccImP,

Z,

4)

(9)

where ci are the coefficients that are determined by minimizing the cost and E,( p , z , +) is the ithbasis function.

B. Cost Function and Method of Optimization The cost function used here is a weighted sum of the following: 1) the systemic cost, S , which is the average SAR outside the region to be heated divided by the average SAR inside the region, and 2) the acute cost, A , which increases rapidly when the SAR outside the region to be heated increases above the average SAR in the region to be heated. That is, the cost function is cost = h,yS

+ h,A

(10)

where h,Tand h, are adjustable parameters. The systemic cost is

+

n

where the subscripts o and i on both the SAR and the integrals indicate outside or inside the region to be heated, respectively, and the bar above the SAR indicates volume average. The acute cost is A =

where

P ( P ,z ,

4) dz!

(12)

1249

CHOWDHURY AND HILL: NUMERICAL OPTIMIZATION OF 3-D SAR DISTRIBUTIONS

-

whenever SAR ( p , z , 4) > SARi and is zero otherwise. By varying h, and ha the relative cost of systemic and acute heating can be varied. Problems due to systemic heating [ 1714191 and problems due to acute heating [20], [21] are often treatment limiting. The cost function is minimized by using a quasi-Newton algorithm [22] which computed the optimized coefficients for the basis functions. The weights of the basis function, ci's (9), are complex multipliers of the individual basis functions. By keeping one of the weights real, (2k - 1) variables are optimized where k is the number of basis functions. IV. RESULTSAND DISCUSSION A , Optimization with Eigenmode Basis Functions Optimized SAR distributions in a MCC model are shown in Figs. 1-5 for a frequency of 70 MHz. The model geometry and constitution is shown in Table I. The permittivity, density, and specific heat data are obtained from [23] and [24]. The water layer in the model simulates the bolus often used for cooling and for achieving better coupling with the source. In Figs. 1-5 the fields are expanded in 40 eigenmode basis functions. Results for three different cost functions are computed. The kz's in the eigenmode basis functions range from 0 to 15 m-' with an interval of 5 m-'. For each k,, the azimuthal index n runs from 0 to 4 giving a total of 20 basis functions for each type of mode (TE or TM at the core of the cylinder). Computed distributions are shown for two different regions of desired maximum SAR. The first region is centered at p = 5 cm, 4 = 0" and z = 0 cm. The region is enclosed by the surfaces p = 2.5 cm and p = 7.5 cm in the radial direction, z = -2.5 cm and z = 2.5 cm in the axial direction and q5 = - 18" and q5 = 18" in the azimuthal direction. The second region is centered at 2.5 cm and is enclosed by the surfaces p = 0 cm and p = 5 cm in the radial direction, z = -2.5 cm and z = 2.5 cm in the axial direction and q5 = - 18" and q5 = 18" in the azimuthal direction. In each of the surface plots of Figs. 1-5, (a) shows the optimized SAR distribution on one of the two symmetric halves of the z = 0 plane and (b) shows the distribution on one of the two symmetric halves of the y = 0 plane. The discontinuities in Figs. 1-5 are a result of two factors: 1) the discontinuity of the radial field at the interface between materials of two different permittivities [l], [13], and 2) the discontinuity of the conductivity. Fig. l(a) and (b) shows the optimized SAR distribution for the first case where the center of the region to be heated is 5 cm from the center. No cost is assigned to localized heating (h, = 0). The optimized SAR distribution is localized in all three directions. At the outer surface of the cylinder the variation of the SAR appears to be approximately Gaussian in the axial and the azimuthal directions. Distributions that appear Gaussian are characteristic of all the optimized results shown here.

r--

L

@ ... ..

_--

...... . . ... . .:

, .

(b) Fig. 1. Optimized SAR distribution with 40 basis functions. The desired region of large SAR was centered at x = 5 cm, y = 0 cm, and z = 0 cm. (a) Shows one symmetric half of the z = 0 plane and (b) shows one symmetric half of the y = 0 plane of the model. The cost is calculated with h, = 0 and h , = 1. The range of kLis from 0 to 15 m-'with a step of 5 m - ' . For each k:, the n goes from 0 to 4 with unit step for both TE and TM modes in the core. The left insert shows the cross section of the model that applies to all the results. The arrow shows the point of intended focus. As in all the figures the frequency is 70 MHz.

The optimized SAR distribution for the second case, i.e., the region centered at 2.5 cm off axis is shown in Fig. 2. The cost function and basis functions were the same as those used for Fig. 1. Figs. 1 and 2 indicate that when no cost is assigned to acute heating, then a source which has a significant amplitude around the body is not ideal, even when the center of the region to be heated is only 2.5 cm from the axis. However, Figs. 1 and 2 also show that where the location of the region to be heated is closer to the axis, the optimized source distribution on a surface outside the cylinder spreads in both the axial and azimuthal directions. This gradual spread of the surface fields as the center of the region to be heated moves from the axis is more pronounced in Figs. 4 and 5 where the cost for the regional acute heat is much larger. When the region is on the axis, the optimal source distribution is uniform in q5 and is more spread axially. As the cost for acute heat is raised, by increasing h,,

1

1250

IEEE TRANSACTIONS ON BIOMEDICAL ENGINEERING, VOL. 38, NO. 12, DECEMBER 1991

z

t

(a) z

t

(b) Fig. 2. Optimized SAR distribution with basis functions and cost function as described in Fig. 1, but with the center of the desired region of maximum SAR at x = 2.5 cm, y = 0 cm, and z = 0 cm.

the optimized SAR distribution spreads axially as shown in Figs. 3 and 4. In these figures the region of heating is centered at 5 cm from the axis as was the case in Fig. 1. The weighting term for the acute heat, h,, is 0 in Fig. 1, 15 in Fig. 3, and 30 in Fig. 4. The parameters used in optimizing the SAR shown in Fig. 5 are identical to those used in Fig. 4, except that the region where power is to be deposited is centered 2.5 cm from the cylinder axis. As the center of the region shifts from 5 cm (Fig. 4) to 2.5 cm (Fig. 5) the fields at the surface spread, primarily in the azimuthal direction. The current J, on a perfectly conducting cylindrical surface 12.5 cm from the cylinder axis was obtained from the tangential components of the magnetic fields, Hz and H,, at the surface by applying the boundary condition. The charge distribution on the surface was computed using the boundary condition for the normal component of E at a perfect conductor. The current and charge distributions associated with the SAR distributions shown in Figs. 4 and 5 are shown in Figs. 6 and 7, respectively. The distributions are shown for z = 0 to I = 30 cm and for 4 = 0 to 4 = 180". In each of these figures, (a) shows the current distribution and (b) shows the charge distribution at time t = 0. Both

n

(b) Fig. 3. Optimized SAR distribution with basis functions and the desired region of maximum SAR described as in Fig. 1, but the cost is calculated with ha = 15 and h, = 1.

the current and charge distributions are symmetric about the plane 4 = 0 to 4 = 180". The current distribution is symmetric above and below the z = 0 plane. The phases for the charge distribution are almost uniform for z > 0 and z < 0 but there is a phase difference of T between these two regions. In Figs. 6(a) and 7(a) the directions and magnitudes of the currents at ut = 0 are shown by the lines with an arrow head where w is the angular frequency. The lines with a circle at the end show the directions and magnitudes of the currents at ut = 7r/2. The locus traced each period by the tip of the current vector is also shown in the plot. For example, in Fig. 6(a) all the current vectors of z = 0.0 cm or 4 = 0" are linearly polarized, i.e., they oscillate in one direction only. However, the current at most other points is elliptically polarized. In cases shown in Figs. 6 and 7, where the maximal desired SAR is in the 4 = 0" plane, the maximum current and charge do not occur at 4 = 0", but occur closer to 4 = 60" - 70" when the center of the region where power is to be deposited is 5 cm from the axis, and occur in the range 4 = 60" to 180" when the center of the region is 2.5 cm from the axis. The electric fields (not shown) at the p = 12.5

1251

CHOWDHURY A N D HILL: NUMERICAL OPTIMIZATION OF 3-D SAR DISTRIBUTIONS

(b) Fig. 3. Optimized SAR distribution with basis functions and the desired region of maximum SAR as described in Fig. 1, but the cost is calculated with h, = 30 and h, = 1.

cm surface are however largest at q5 = 0" when the region is centered at p = 5 cm and are relatively uniform from 4 = 0" to 90" when the region is at p = 2 . 5 cm. The total E field is a sum of the fields generated by the charges and currents and depends on the phase relations between the sources. The results, shown in Figs. 6 and 7 , that the current and charge are smaller closest to the region of desired power deposition, may appear counter-intuitive. The results illustrate that optimization of 3-D SAR distributions in layered structures is not simple. Even in a homogeneous problem a particular source may contribute more to the generation of fields in the region of desired power deposition or may contribute more to canceling unwanted fields in a region where no power deposition is desired. The layered problem is more complex in that power may be reflected or guided by the interfaces, and often will not take the shortest distance between two points. So it is even more difficult to obtain a simple intuitive understanding of the relation between the optimized SAR and the currents and charges at a surface outside the body. The results of Figs. 4-7 indicate that if the SAR is to be optimized, both the phases and the amplitudes of the

1

(b) Fig. 5. Optimized SAR distribution with basis functions used in Fig. 1, but the desired region of maximum SAR is centered at x = 2.5 cm, y = 0 cm, and z = 0 cm and the cost is calculated with h , = 30 and h,7 = 1.

sources need to be controlled. Similar results have been observed in 2-D numerical studies [ 5 ] , [8], [4], [25] and in experimental measurements [26], [27].

B. Optimization with Combined Basis Functions The optimizations described in Section IV-A are computationally intensive because many eigenfunctions were used. In order to reduce the number of variables, we generated combined basis functions in which the fields are confined axially and usually azimuthally and/or radially (see Section 111-A). Expressions for the amplitude distributions at p = 12.5 cm are given in Table I1 for the combined basis functions used. The phase of the function at the p = 12.5 cm surface is chosen so that the fields add in phase at a point 5 cm from the axis. The optimal SAR obtained with ten combined basis functions is shown in Fig. 8 in which the cost function is the same as that of Fig. 4, i.e., h, = 30 and h, = 1.0. The SAR distributions of Figs. 4 and 8 appear similar in most respects, although the SAR in the outer fat layer is much larger in Fig. 8 than it is in Fig. 4 . To further compare the optimized distributions obtained with the two different types of basis functions, contour plots of the distributions of Figs. 4 and 8 are shown in Fig. 9. The dotted lines show the optimized SAR dis-

I

IEEE TRANSACTIONS ON BIOMEDICAL ENGINEERING. VOL. 38. NO. 12. DECEMBER 1991

1252

TABLE 1 A N D CONSTANTS USEDFOR THE MCC MODELAT 70 MHz (M = MUSCLE, F DESCRIPTION 1 /2M + 1 / 2 F 3 VISCERA A N D 1 /3M + 2 / 3 F = SKIN) Radial Width (cm)

Constitution

1/2M

+ 1/2F

0-2.5 2.5-10 10-1 1 11-12 12-12.5

M F

1/3M + 2 / 3 F Water

Relative Permittivity

t

47.3 83.0 11.5 35.3 78.0

FAT,

Conductivity U S . m-'

Density g . cm '

Spec. heat C,,caI . g . " C - '

0.425 0.800 0.050 0.300 0.000

0.96 1.02 0.90 0.94 1 .oo

0.55 0.86 0.24 0.45

'

I .oo

. 4

i

4 4

4 I

I

0

90

I 180

q5 in degrees (a)

(b)

(b)

Fig. 6. Distribution of (a) current and (b) charge on one quarter of a perfectly conducting cylindrical surface at p = 12.5 cm associated with the SAR distribution shown in Fig. 4. In (a) the lines with arrowheads show the direction and magnitude of current at t = 0 and the lines with circles at the end show the direction and magnitude of current at w t = ?r/2 where w is the angular frequency. In (b) the charge distribution is shown at f = 0.

Fig. 7. Distribution of (a) current and (b) charge on the sample plane as in Fig. 6, but are the distributions associated with the SAR distribution shown in Fig. 5. The lines with arrowhead and circles have the same meaning as in Fig. 7. In (b) the charge distribution is shown at t = 0.

tribution with ten basis functions, corresponding to Fig. 8, and the solid lines show the optimized SAR distribution with 40 eigenmodes, corresponding to Fig. 4. The SAR for each case is normalized to unity at p = 5 cm, z = 0

1

cm and $I = 0 " . In Fig. 9(a), which shows the SAR in the z = 0 plane, the contours are similar, especially for the larger values of the SAR. Fig. 9(b) shows the contours for the SAR in the $I = 0" plane which are also very similar, except for the outer layer of fat. The relatively close agreement between the two results suggests the possibility

1253

CHOWDHURY AND HILL: NUMERICAL OPTIMIZATION OF 3-D SAR DISTRIBUTIONS

TABLE I1 MATHEMATICAL EXPRESSIONS FOR THE AMPLITUDE DISTRIBUTION AT po = 12.5 cm FOR THE COMBINED WHERE BASISFUNCTIONS USEDI N THE OPTIMIZATION z IS IN cm A N D 6 I S I N DEGREES. THEPHASEFOR ALL THE BASISFUNCTIONS AT po = 12.5 cm PLANF IS CHOSEN so AS TO MAKETHE FIELDS ADDI N PHASEAT x = 5 cm, z = 0 cm, A N D y = cm Basis

12

c 6 x

.d

Amplitude Distribution in z and 4

1

n

0

-6

v-12

2 3 4 5 6 7 8

6

12

x in cm (a)

9 10

24 I '

I !

8

18 -

.*c N

12 -

II 4 II

2

x in cm (b) Fig. 9. Contour plots for the SAR distributions shown in Figs. 4 and 8 are superimposed on top of each other. The dotted lines show the contour plot for the SAR distribution shown in Fig. 8 and the solid lines show the contour plot corresponding to SAR distribution shown in Fig. 4. The contours are shown (a) in one symmetric half of the z = 0 plane and (b) in one symmetric half of the y = 0 plane.

v.

(b) Fig. 8. Optimized SAR distribution with ten basis functions as defined in Table 11. The desired region of maximum SAR is centered at x = 5 cm, y = 0 cm, and z = 0 cm. SAR distribution is shown (a) in one symmetric half of the z = 0 plane, and (b) in one symmetric half of the y = 0 plane of the cylinder. The cost is calculated with h,, = 30 and h, = 1. The arrow shows the point of intended focus.

of using the combined basis functions for preliminary investigations. However, the significant difference between SAR's in the outer fat layer indicates that the ten combined basis functions used here cannot expand the fields as well as the 40 eigenmodes.

1

DISCUSSION OF THE APPROACH AND LIMITATIONS The objective of the work was to determine the best SAR distributions in one type of model, the MCC model. Some limitations of the work presented are as follows. 1) MCC models are not particularly good approximations to the human body or even the torso. They have no inhomogeneities in the axial direction. Finite element [ 113 or finite difference time domain [ 121 methods would be more appropriate for modeling more realistic body models. 2) The cost function and desired SAR distribution were only estimated. There are two problems here. First, we only model the SAR, not the temperature. But the cost function should be defined in terms of the temperature distributions and possibly other factors. As stated in the introduction we assume that the desired SAR distributions are ones that could have been determined using bioheat models, and we assume that the cost functions could have been approximated from cost functions based on temperature. But even if the cost function and desired SAR distribution were very consistent with those obtained using thermal models, the

1254

IEEE TRANSACTIONS ON BIOMEDICAL ENGINEERING, VOL. 38. NO. 12. DECEMBER 1991

results could not be as good as if the temperature distribution was optimized; optimization of subsystems is not as good as optimization of the whole system. Second, the data required to determine an accurate cost function, based on temperature and possibly other parameters, have not yet been obtained and may not be obtained for years. Nevertheless, the approach presented and illustrated here could be useful for evaluating and designing hyperthermia applicators. The desired SAR distributions and cost functions are not perfect, but they are consistent with what is attempted in hyperthermia treatments, i.e., in order to generate a temperature distribution, energy is deposited in some regions and not others. The models are less than ideal, but they are inhomogeneous, they could be readily fabricated for testing, and in these models it is possible to compute 3-D nonaxisymmetric SAR distributions with a high degree of accuracy. MCC models should also be useful for verifying the accuracy of computations made using finite methods. Two questions are raised by the results shown here. First, could existing hyperthermia applicators generate the optimized SAR distributions shown? This question could be addresssed either with mathematical or physical models. We suspect that even a 16-element dipole array would not have the flexibility to generate the optimized SAR distributions shown here, but without further modeling we cannot say if the differences are large enough to be important. Further work is required. Second, is it possible, and if so how difficult is it, to generate the current and charge distributions associated with these ideal SAR distributions? Our speculation is that the current and charge distributions could be generated with multielement arrays. We see nothing in the field, current, or charge distributions that make us think otherwise. Elements of the arrays may include pairs of dipoles oriented perpendicular to each other so that the polarization as well as the amplitude and phase can be controlled. Elements whose main purpose is to hold charge, such as capacitor plate applicators, may also be required. VI. SUMMARY This paper addresses the following question: given a MCC model of the body, and given a nonaxisymmetric desired SAR distribution and a cost function for assessing how close a particular SAR distribution is to the desired distribution, what is the best possible 3-D SAR distribution consistent with Maxwell’s equations? Gradient search optimization was used to find the coefficients of the fields for the optimized SAR. The computed charge and current distributions on a conducting surface outside the model, i.e., on a surface that might contain the sources of the optimized fields, are dependent upon the radial position of the region to be heated and upon the cost function used. As the region to be heated moves from the center, the variations in phase and amplitude distributions in the azimuthal direction are more pronounced. Some limitations of the results pre-

sented here are: 1) the models of the body are less realistic than those that could be modeled using finite element or finite difference time domain methods, 2) no thermal modeling was done, i.e., only the SAR was computed, and 3) the cost functions were only estimated. Further work to see how close existing applicators could come to generating the SAR distributions presented here is warranted.

ACKNOWLEDGMENT We thank P. Barber and P. Fessenden for helpful discussions. A major part of this research was conducted using the Northern Parallel Architectures Center at Syracuse. Computations were also performed at the Cornel1 National Supercomputer Facility, a resource of the Corne11 Theory Center, which is funded in part by the National Science Foundation, New York State, the IBM corporation and members of the Center’s Corporate Research Institute.

REFERENCES [ l ] D. A. Christensen and C. H. Durney, “Hyperthermia production for cancer therapy: A review of fundamentals and methods,” J . Microwavel‘ower, vol. 16, no. 2, pp. 89-105, 1981. [2] D. Q . Chowdhury and S . C. Hill, “Limits on focused power deposition for electromagnetic hyperthermia,” Int. J . Hyperthermia, vol. 7, pp. 185-196, 1991. [3] J . B. Anderson, “Theoretical limitations on radiation into muscle tissue,” Int. J . Hyperthermia, vol. I , pp. 45-55, 1985. [4] G . Arcangeli, P. P. Lombardini et a l . , “Focusing of 915 MHz electomagnetic power on deep human tissue: A mathematical model study,” IEEE Trans. Biomed. Eng., vol. BME-31, pp. 47-52, 1984. [ 5 ] J. W. Strohbehn, E. H. Curtis, K. D. Paulsen, and D. R. Lynch, “Optimization of the absorbed power distribution for an annular phased array hyperthermia system,” Int. J . Radiation Oncol. Biol. Phys., vol. 16, pp. 589-599, 1989. [6] C. M. Rappaport and F. R. Morganthaler, “Optimal source distribution for hyperthermia at the center of a sphere of muscle tissue,” IEEE Trans. Microwave Theory Tech., vol. MTT-35, pp. 1322-1327, 1987. “Localized hyperthermia with electromagnetic arrays and leaky[7] -, wave troughguide aoolicator.” IEEE Trans. Microwave Theory Tech., vol. MTT-34, pp. 636-642, 1986. N. Morita, T . Hamasaki, and N. Kumagi, “An optimal method in multi-applicator systems for forming a hot zone inside the human body,” IEEE Trans. Microwave Theory Tech., vol. MTT-34, pp. 532-538, 1986. C . DeWagter, “Optimization of simulated two-dimensional temperature distributions induced by multiple electromagnetic applicators,” IEEE Trans. Microwave Theory Tech., vol. MTT-34, pp. 589-596, 1986. M. Knudsen and U. Hartmann, “Optimal temperature control with phased array hyperthermia system,” IEEE Trans. Microwave Theory Tech., vol. MTT-34, pp. 597-603, 1986. D. R. Lynch and K. D. Paulsen, “Time-domain integration of the Maxwell equations on finite elements, ” IEEE Trans. Antennas Propagar., vol. 38, pp. 1933-1942, 1990. D. Sullivan, “Three-dimensional computer simulation in deep regional hyperthermia using the finite-difference time-domain method,” IEEE Trans. Microwave Theory Tech., vol. MTT-38, pp. 204-21 1 , 1990. S . C. Hill and D. Q . Chowdhury, “Three-dimensional SAR distribution computed in a multilayered cylindrical model for electromagnetic hyperthermia,” IEEE Trans. Microwave Theory Tech., vol. MTT-37, pp. 1192-1199, 1989. P. G . Cottis, G. E. Chatzarakis, and N. K. Uzunoglu, “Electromagnetic energy deposition inside a three-layer cylindrical human body model caused by near-field radiators,” IEEE Trans. Microwave Theory Tech., vol. MTT-38, pp. 990-999, 1990.

CHOWDHURY AND HILL: NUMERICAL OPTIMIZATION OF 3-D SAR DISTRlBUTlONS

[IS] R. Harrington, Time-Harmonic Electromagnetic Fields. New York: McGraw-Hill, 1961, ch. 5 . [I61 J . P. Mason, “Cylindrical Bessel functions for a large range of complex arguments,” Computer Phys. Commun., vol. 30, pp. 1-1 1, 1983. [17] D. S . Shimm, T. C. Cetas, J. R. Oleson, J . R . Cassady, and D. A. Sim, “Clinical evaluation of hyperthermia equipment: The University of Arizona institutional report for the NCI hyperthermia equipment evaluation contract,” Int. J . Hyperthermia, vol. 4, no. 1 , pp, 39-51, 1988. 1181 D. S . Kapp, P. Fessenden, T . V . Samulski, M . A. Bagshaw, R. S. Cox, E. R. Lee, A. W . Lohrbach, J. L. Meyer, and S.D. Prionas, “Stanford University Institutional report. Phase 1 evaluation of equipment for hyperthermia treatment of cancer,” Int. J . Hjperrhermia, vol. 4 , no. 1, pp. 75-115, 1988. [I91 M . D. Sapozink, F. A. Gibbs Jr., P. Gibbs, and J. R. Stewart, “Phase I evaluation of hyperthermia equipment- University of Utah iristitutional report,” Int. J . Hyperfhermia, vol. 4 , no. 1, pp. 117-132, 1988. (201 M. J . Hagmann and R. L. Levin, “Aberrant heating: A problem in regional hyperthermia,” IEEE Trans. Biomed. Eng., vol. BME-33, pp. 405-41 1 , 1986. [21] P. F. Turner, “Hyperthermia and inhomogeneous tissue effects using an annular phased array,’’ IEEE Trans. Microwave Theory Tech., vol. MTT-32, pp. 874-882, 1984. [22] R. B. Schnabel, J . E. Koontz, and B . E. Weis, “A modular system of algorithms for unconstrained minimization,” ACM Trans. Mathemat. Software, vol. 11, pp. 419-440, 1985. [23] C.H. Durney et al., Radiofrequency Radiation Dosimetry Handbook, 2nd ed. Salt Lake City, UT: Univ. Utah, 1978. [24] D. R . Lynch, K. D. Paulsen, and J . H. Strohben, “Finite element solution of Maxwell’s equations for hyperthermia treatment planning,” J . Comput. Phys.. vol. 58, pp. 246-269, 1985. [25] V. Sathiaseelan, M. F. Iskander, G. C. W. Howard, and N. M. Bleehen, “Theoretical analysis and clinical demonstration of the effect of power pattem control using the annular phased-array hyperthermia system,” IEEE Trans. Microwave Theory Tech., vol. MTT-34, pp. 514-519, 1986. [26] T. V. Samulski, D.S . Kapp, P. Fessenden, and A. Lohrbach, “Heating deep seated eccentrically located tumors with an annular phased array system: A comparative clinical study using two annular array

R

1255

operating configurations,” Int. J . Radiation Oncol. Biol. Phys., vol. 13, pp. 83-94, 1987. [27] P. F. Turner, “Regional hyperthermia with an annular phased array,” IEEE Trans. Biomed. Eng., vol. BME-31, pp. 106-114, 1984.

Dipakbin Q. Chowdhury (S’89) was born in Bangladesh on January I , 1964. He received the B.Sc. degree in electrical and electronic engineering from Bangladesh University of Engineering and Technology, Dhaka, Bangladesh, in 1986, and the M.S. and Ph.D. degrees in electrical and computer engineering from Clarkson University, Potsdam, NY, in 1989 and 1991, respectively. He is now a Research Associate with the Department of Applied Physics and Center for Laser Diagnostics, Yale University, New Haven, CT. His interest is in computational electromagnetics and applications to biomedical engineering, light scattering by small particles, and nonlinear optics. His present research area is nonlinear optics in small particles.

Steven C. Hill (S’82-M’85) received the B.S. degree in chemistry and the M.S. and Ph.D. degrees in electrical engineering from the University of Utah, Salt Lake City, in 1973, 1982, and 1985, respectively. He also received the M.S. degree in pharmacology from Yale University, New Haven, in 1977. He was an assistant professor at Clarkson University, Potsdam, NY from 1985 to 1991. He is now an Adjunct Associate Professor of Electrical and Computer Engineering at New Mexico State University, Las Cruces, NM, and a National Research Council Associate with the Army Atmospheric Sciences Laboratory at White Sands, N M . His interests are in numerical electromagnetics and optics, nonlinear optics, and in their application to biomedical problems.