The role of organic carbon, iron, and aluminium oxyhydroxides as trace metal carriers: Comparison between the Trinity River and the Trinity River Estuary (Galveston Bay, Texas)

Phys. Med. Biol. 48 (2003) 1173–1190 PII: S0031-9155(03)56193-0

Treatment time reduction for large thermal lesions by

using a multiple 1D ultrasound phased array system

Hao-Li Liu1_{, Yung-Yaw Chen}1_{, Jia-Yush Yen}2_{and Win-Li Lin}3 1_{Department of Electrical Engineering, National Taiwan University, Taipei, Taiwan,}

Republic of China

2_{Department of Mechanical Engineering, National Taiwan University, Taipei, Taiwan,}

Republic of China

3_{Institute of Biomedical Engineering, National Taiwan University, Taipei, Taiwan,}

Republic of China E-mail: yychen@cc.ee.ntu.edu.tw Received 11 November 2002 Published 22 April 2003 Online atstacks.iop.org/PMB/48/1173 Abstract

To generate large thermal lesions in ultrasound thermal therapy, cooling intermissions are usually introduced during the treatment to prevent near-field heating, which leads to a long treatment time. A possible strategy to shorten the total treatment time is to eliminate the cooling intermissions. In this study, the two methods, power optimization and acoustic window enlargement, for reducing power accumulation in the near field are combined to investigate the feasibility of continuously heating a large target region (maximally 3.2× 3.2 × 3.2 cm3_{). A multiple 1D ultrasound phased array system generates} the foci to scan the target region. Simulations show that the target region can be successfully heated without cooling and no near-field heating occurs. Moreover, due to the fact that there is no cooling time during the heating sessions, the total treatment time is significantly reduced to only several minutes, compared to the existing several hours.

1. Introduction

Over 60% of cancer patients do not discover breast or liver tumours until they have reached a size of 2–3 cm in diameter, or more (Kahn et al 1983,Leibel and Phillips 1998). From a treatment standpoint, it is desirable to treat the entire tumour in a single treatment session. One promising treatment method is ultrasound thermal therapy, whose focused ultrasound applicators generate a single lesion that is usually very small (approximately 1–2 mm in diameter). Hence, to generate large thermal lesions, the focus is scanned, mechanically or electrically, over the tumour region (Malcolm and ter Haar 1996, Wan et al 1996). The alternation of focus areas usually requires cooling intermissions to avoid excessive power

accumulation in front of the target region, or what is commonly referred to as a ‘near-field heating problem’ (Damianou and Hynynen 1993). Due mainly to such cooling intervals, the time required to heat a 3 × 3 × 3 cm3 volume is several hours (Fan and Hynynen 1996). Similar cases of equivalent treatment time have also been reported in the literatures (Daum and Hynynen 1999,Wu and Sherar 2002).

Obviously, one way to reduce the time required to heat the tumour in the previous method would be to lessen, or even eliminate, the cooling intervals. Indeed, two methods have been proposed. The first method minimizes the ultrasound energy focused on the target region to reduce the energy accumulation in the near-field (Wan et al 1996, Daum and Hynynen 1999). A second method enlarges the geometrical gain of the ultrasound beam (defined as the ratio between the acoustic window size and the target region size) (Lin et al 1999, Hutchinson and Hynynen 1996). Although these two methods are both effective for treatment-time reduction, neither approach can be used to continuously heat a lesion with a volume larger than 5 cm3_{. However, we believe a combination of these two strategies would allow treatment} of a large target region, without cooling, while retaining the advantage of a reduced treatment time.

In this study, we use a multiple 1D ultrasound phased array system to investigate the feasibility of combining the two methods described above to continuously heat and generate a large thermal lesion (maximally 3.2× 3.2 × 3.2 cm3). The transducer system is designed to change the geometrical gain of the ultrasound beam easily by altering the array number. Moreover, compared to spherical transducers, the focus generated by a single 1D phased array is relatively large (Hutchinson and Hynynen 1996), which in turn reduces the total number of scanned foci. We also developed an algorithm to determine the power weightings of scanned foci as well as the total scanning time, so as to optimize thermal dose conformability. Results show that we can generate a thermal lesion reaching 3.2× 3.2 × 3.2 cm3_{, and at the same} time significantly reduce the total treatment time to only several minutes. Parametric studies were also performed, including treatment effects for different target volumes, influence of the blood perfusion rate and the switching speed between foci.

2. Methods

2.1. Pressure field calculations

To calculate the ultrasonic pressure field, the transducer is modelled as a grid of point sources. Then, utilize the Rayleigh–Sommerfeld integral to sum up the contribution of each point source,r
, to the point of the field, at r (O’Nell 1949):

p(x, y, z) =iρck 2π
S u e−(α+ik)(r−r
₎ r − r
dS (1)

whereρ is the tissue density (1050 kg m−3), c is the speed of sound (1500 m s−1), k is the wave number (2π/λ, where λ is the wavelength), u is the complex surface velocity of the source andα is the attenuation (4.1 Np m−1MHz−1). In an attenuated medium, the power deposition,

q, for the desired volume is given as (Nyborg 1981)

q = 2αI I = p2

2ρc (2)

where I is the power intensity of the field. The computations of the acoustic pressure and the resulting power deposition were performed with the source grid size ofλ/6. The nonlinear, refraction and scattering effects on the tissue/water interface were not included due to the similar acoustic properties of water and tissue (Lu et al 1996).

Target volume 16 cm 11 cm Target volume 1.5 cm 1 5 8 9 12 4 45˚ R_x= 15 cm R_y= 20 cm 13 cm 11 cm

Y

Z

o

X

Z

o

Figure 1. Geometry of the multiple 1D phased array system in (a) side view and (b) top view.

2.2. Temperature and thermal dose calculations

The tissue temperature response, T, can be simulated using the well-known bio-heat transfer equation (Pennes 1948)

ρct∂T

∂t = k∇2T − wbcb(T − Tar) + q (3)

where ctand cbare the specific heats of tissue and blood (both set to be 3770 J kg−1 ◦C−1), k is the thermal conductivity (0.56 W m−1 ◦C−1), wbis the blood perfusion rate (1 kg m−3s−1 in our simulation unless otherwise mentioned) and Tar is the arterial blood temperature (37◦C). This equation is solved using a numerical finite difference method with all boundary and initial conditions set to 37 ◦C. The time step and the grid spacing in the x-, y- and

z-directions are 50 ms, 0.5 mm, 1 mm and 1 mm, respectively.

The thermal dose (TD), in terms of equivalent minutes at 43◦C, is used to estimate the necrosed tissue volume, and is calculated numerically using the following equation (Dewey 1994,Pearce and Thomsen 1995,Sapareto and Dewey 1984):

TD=
_tf t0 R(T −43)_{dt ≈}tf t0R (T −43)_{t (4)}

where R= 2 for T
43◦C, and R= 4 for 37◦C< T < 43◦C,t is the time step, and t0and

tfrepresent the initial and final times, respectively. The TD value required for total necrosis

ranges from 25 to 240 min for brain to muscle tissues (Damianou et al 1995,Dewey 1994). In this study, a value of TD= 240 min is chosen as the threshold for complete heating.

2.3. Multiple 1D ultrasound phased array system and the heating strategy

A system of multiple 1D ultrasound phased arrays, arranged as shown in figure1, is used to examine the feasibility of continuously heating the target regions. This arrangement is advantageous in its ease of increasing the number of phased arrays to enlarge the acoustic window of the system. The test system consists of 12 identical 80-element, 1D phased arrays, with curvature radius set to 15 cm. The driving frequency is assumed to be 800 kHz, which is used to provide high penetration into tissue suggested by Hutchinson and Hynynen

Depth direction (cm) 8 9 10 11 12 13 Transverse direction (cm) -3 -2 -1 0 1 2 3 (b) Depth direction (cm) (a) 8 9 10 11 12 13 -3 -2 -1 0 1 2 3 Longitudinal direction (cm) start

switch to next array

Figure 2. (a) Side view and (b) top view of focal positions of a single array from group 2 (denoted as dots) and corresponding sample relative ultrasound intensity fields (50% and 70% contours). The focal positions are set with a 0.2 cm spacing in the longitudinal direction and a 0.5 cm spacing in the depth direction. Arrows indicate the scanning direction of this 1D phased array.

(1996). Another advantage is that, compared to higher driving frequency, the generated focal dimension is relatively larger, and this is advantageous to reduce the total foci number to scan a large PTV. Because the 1D phased arrays lack focusing ability in the y-direction, they are arranged to project overlapping power at the position of the target region, as shown in figure1(b).

During heating, the 1D phased arrays are alternately enabled in a circular pattern, and the phase signals of the 80 elements in one array are activated to generate the focus to scan the target region. The phase differences of these 80 driving signals are determined using the ‘in-phase’ method proposed in a previous study (Hutchinson and Hynynen 1996). A typical focal position and scanning sequence for a single 1D phased array is shown in figure2. Power depositions of scanned foci are also shown in figure2(for clarity of the figures, only 45 out of 135 foci are shown in figure2(a), and 5 out of 9 foci are shown in figure2(b)). The –3 dB dimension of focus along the longitudinal direction and depth is about 0.2 cm and 0.5 cm on average. Therefore, the focal spacing is set to 0.2 cm (17 points total) along the longitudinal direction, and 0.5 cm (9 points total) along the depth, to provide enough lesion overlap between foci (Damianou et al 1993).

In order to shorten the overall scanning period, the nth and the (13 – n)th arrays are triggered with the same driving signals simultaneously, taking advantage of x- and y-axis symmetry in the system. The dwelling time for each focal pair is set to 50 ms (i.e. 20 Hz) unless otherwise mentioned. To simplify the calculation, delay and power loss caused by switching of the 1D phased arrays are not considered. The overall heating process is performed continuously, without introducing any cooling periods.

2.4. Optimization of thermal dose conformability

To optimize ultrasound power, metrics for evaluating the treatment results are defined, and an optimization algorithm is developed. Conformability between the thermal lesion and the target region is the optimization goal in this study. Here, the conformal index used in 3D conformal radiation therapy, defined as the difference between a coverage index (CI) and an

external volume index (EI), is utilized (Saw and Suntharalingam 1990): Conformal index= CI − EI CI= PTVTD
240 PTV × 100% EI= OTV_TD
240 PTV × 100% (5) where PTV_TD
240is the volume inside the planning target region (PTV) with TD
240 min, and OTVTD
₂₄₀is the volume outside the PTV with TD
240 min. A high conformal index value implies more heating, and thus more thermal damage, is incurred within the target region than outside the target region.

To maximize the conformal index value, we must determine the relative power weightings of scanned foci at different positions u1, u2, . . . , uN, and the total heating time, ton. However,

calculation of the thermal dose is quite time-consuming, so it is nearly impossible to directly alter u1, u2,. . . , uNto evaluate the resulting performance. Instead, the heating effect of an

averaged power deposition, Q, will be close to the effect of sequentially scanning when the switching frequency is high, which can then be used to simplify this optimization problem:

Q ∼= 1

N(u1h1+u2h2+· · · + uNhN) (6)

where h1, h2,. . . , hNare the power depositions of different foci.

Thus, the optimization problem becomes simply determining Q and ton to obtain a maximal conformal index value. The number of variables is then greatly reduced. To solve this simplified problem, a simple direct search algorithm is used (Relkaitis et al 1983). In each iteration, when Q and tonare adopted, the corresponding u1, u2,. . . , uNfor the given

Q are solved to calculate the resulting conformal index value. To solve this, a constrained

least-squares optimization method (Golub and Van Loan 1996) is applied. Details of this adaptation method are described in the appendix. The overall optimization flow is shown in figure3.

To narrow down the search space in this optimization process, Q is uniformly distributed inside the target volume and is relatively higher on the boundary (Lin et al 2001, Ocheltree and Frizzell 1988,Wan et al 1999). The power depositions at the boundary and inside the target region are both tunable in iterations. Moreover, after finding u1, u2,. . . , uN,

they are scaled to set the maximum instantaneous power intensity at 500 W cm−2. This increases the allowable maximum absorbed power density under the cavitation threshold (Hynynen 1991).

3. Simulation results

3.1. Effect of adjusting foci weightings

To examine the effect of optimizing the ultrasound energy, we compare the power depositions with and without adjusting the power weightings. For the case of adjusted power weightings, the weightings of the foci are determined according to the optimization flow in figure3and described above. For the non-optimized case, the power weightings are all identical. In both cases, 12 1D phased arrays are used, and a volume of 3.2× 3.2 × 3.2 cm3_{, centred at}

z= 11 cm, is taken as the target region. Figures4(a)–(d) show the resulting averaged power depositions, and figure5shows the number of foci used for each 1D phased array. Results show that adjusting the foci weightings effectively reduces the averaged power accumulation in front of the target region. The power distribution can also be regulated to approach the target volume, and the power near the boundary can be effectively raised in the x-direction, as shown in figure4(d). The raising effect of power deposition in the y-direction is not apparent

Select target volume size

•Give initial power deposition and

sinitial scanning time ton

•Calculate initial u1, u2, ...., uN from susing constrained least squares method •Calculate initial Conformal index

•Change boundary/centre power ratio of by step size •Change scanning time ton by t

•Calculate new u1, u2, ...., uN and omit foci with negative sor zero weightings

•Calculate (Conformal index)new

• t = t/2

•Divide power ratio step by 2

Yes Yes

No

No (Conformal index)new

> (Conformal index)?

t and power ratio step size less than exit criteria?

Start

Finish Qˆ

Qˆ

Qˆ

Replace original Conformal index, , and tQˆ onwith new ones

Figure 3. Flow chart of optimization algorithm used to adjust the foci weightings and the heating time.

from figure4(c). This is caused by the lack of focal ability in the y-direction of the system. The total number of foci used, after adjusting the power weightings, can be reduced to 46% of the original one to shorten the scanning period.

3.2. Effect of different aperture sizes

To examine the near-field heating due to the small acoustic window, only 4 strips (5–8) and 8 strips (3–10) of the 1D phased arrays are compared with the 12-strip case. The power weightings of foci in the 4-strip and 8-strip cases were both adjusted, and the averaged power depositions are shown in figure6. When compared with figures4(c) and (d), we find that the power extension towards the skin in the 12-strip case is smaller than for the 4- and 8-strip cases due to a larger geometrical gain.

Heating results produced by 12 1D phased arrays with adjusted power weightings (i.e., both optimizing the power and enlarging the geometrical gain) are shown in figure7. The upper row of figure7shows the temperature distributions at the time when the power is turned off (210 s), and the lower row shows the resulting thermal dose distributions. The required heating times are short due to the absence of cooling intermissions. Moreover, the extension

6 8 10 12 14 X-D (cm) 6 8 10 12 14 Y-D (cm) 6 8 10 12 14 4 2 0 -2 -4 Z-D (cm) Y-D (cm) 6 8 10 12 14 X-D (cm) (a) Z-D (cm) (c) Z-D (cm) (d) Z-D (cm) (b) 4 2 0 -2 -4 4 2 0 -2 -4 4 2 0 -2 -4

Figure 4. Comparison of averaged power depositions: (a), (b) without adjusting the power weightings of foci and (c), (d) with power weighting adjustment. The contour lines are drawn at 10% intervals from 10% to 90%. The PTV is set to 3.2× 3.2 × 3.2 cm3_{with the centre located at}

z= 11 cm. 0 50 100 150 200 1 2 3 4 5 6 7 8 9 10 11 12 Array Driven foci number

without adjusted (1026 foci pairs) adjusted (496 foci pairs)

Figure 5. Comparison of the foci numbers for each 1D phased array between the cases with and without power weighting adjustment in figure 4.

6 8 10 12 14 -4 -2 0 2 4 X-D (cm) 6 8 10 12 14 -4 -2 0 2 4 Z-D (cm) Y-D (cm) 6 8 10 12 14 -4 -2 0 2 4 Y-D (cm) 6 8 10 12 14 -4 -2 0 2 4 X-D (cm) (a) Z-D (cm) (b) Z-D (cm) (c) Z-D (cm) (d)

Figure 6. (a), (b) Averaged power depositions of the 4-strip case (arrays 5–8) and (c), (d) the 8-strip case (arrays 3–10). The contour lines are drawn at 10% intervals from 10% to 90%. The PTV are set to 3.2× 3.2 × 3.2 cm3_{with the centre located at z= 11 cm.}

of TD> 240 min region over the target region towards skin is restricted to within 0.5 cm. This shows that the combination of optimized power and enlarged geometrical gain allows a large tumour to be heated continuously, without near-field heating. In the following study, 12 1D phased arrays with adjusted power weightings are used to investigate other important factors, including the target region size, the blood perfusion rate and the switching frequency.

3.3. Comparison between different target volumes

To compare the heating results of different target volumes, four cubic-shaped PTV of dimensions 2.0, 2.4, 2.8 and 3.2 cm and centred at z= 11 cm are taken as target regions. For these simulations, the blood perfusion rate is set to 1 kg m−3s−1, and the switching frequency between scanned foci is set to 20 Hz. The upper half of table1 summarizes the treatment results. Figure8 shows surface views of the thermal dose at 240 min for each of the four cases. Strong heat conduction effects make it difficult to reach the thermal dose threshold at the corners of the cubes. From table1, the heating times in all of the cases are short (within several minutes), and increase with target region size. This is due to the total number of foci used increasing as the target region becomes larger (see (6)), which forces a lower average power level to achieve the heating result with a longer treatment time.

4 6 8 10 12 14 4 6 8 10 12 14 -4 -2 0 2 4 4 6 8 10 12 14 -4 -2 0 2 4 -4 -2 0 2 4 4 6 8 10 12 14 -4 -2 0 2 4 -4 -2 0 2 4 -4 -2 0 2 4 -4 -2 0 2 4 -4 -2 0 2 4 Z-D (cm) Z-D (cm) Y-D (cm) Z-D (cm) Z-D (cm) Y-D (cm) Y-D (cm) X-D (cm) X-D (cm) Y-D (cm) X-D (cm) X-D (cm) 50 55 45 60 50 55 60 45 45 60 50 50 1 240 60 240 60 1 1 240 1000 60 1000 1000

Figure 7. Temperature distribution (t= 210 s) and respective thermal dose distributions in the

y–z, x–z and x–y planes (at z= 11 cm) for 12 strips of 1D phased arrays. The blood perfusion rate

is set to 1 kg m−3s−1.

Table 1. Comparison of treatment results for different target volumes and blood perfusion rates (ton= heating time, Tz= 8 cm= maximal temperature at z = 8 cm for t = ton, Tmax= maximal

temperature for t= ton, TDz= 8 cm= maximal thermal dose at z = 8 cm, TDmax= maximal thermal

dose). PTV

wb size/vol Foci no ton Tz= 8 cm Tmax TDz= 8 cm TDmax CI EI CI− EI Volume

(kg m−3s−1) (cm/cm3_{) (pairs) (s) (}◦_{C) (}◦_{C) (min) (min) (%) (%) (%) (cm}3₎

1 2.0/8.0 331 88 40.3 56.7 0.1 7960 85.2 7.5 77.7 6.8 2.4/13.8 400 105 41.7 58.6 0.7 25 401 89.5 9.2 80.2 12.3 2.8/22.0 421 120 42.5 61.5 2.2 1.2× 105 _{93.9 10.6 83.3 20.7} 3.2/32.8 483 210 44.8 65.9 12.6 2.8× 106 _{90.7 7.5 83.2 29.7} 5 2.0/8.0 380 133 41.9 57.5 0.2 10 805 87.5 6.8 80.7 7.0 2.4/13.8 443 158 43.8 59.9 1.5 54 900 94.8 15.2 79.6 13.1 2.8/22.0 478 259 43.9 58.8 4.3 37 510 90.2 7.3 82.9 19.8 3.2/32.8 496 450 44.9 60.5 13.2 5.1× 105 _{89.6 6.0 83.0 29.4} 10 2.0/8.0 402 210 41.1 56.4 0.1 10 621 89.4 8.7 80.7 7.2 2.4/13.8 461 350 41.9 55.4 0.6 20 525 88.0 8.8 79.2 12.1

3.4. Comparison between different blood perfusion rates

The effect of different blood perfusion rates is simulated for values of 5 and 10 kg m−3s−1, and compared with the previous results for a rate of 1 kg m−3s−1. Results are summarized in the lower half of table1, and the heating times for the different rates are shown in figure9.

X-D (cm) Y -D (cm) Z-D (cm) Z-D (cm) X-D (cm) X-D (cm) Y -D (cm) Y -D (cm) X-D (cm) Y -D (cm) Z-D (cm) Z-D (cm) (a) (b) (c) (d)

Figure 8. Isosurfaces of TD= 60 and 240 min for different PTV settings with (a) 2.0 × 2.0 × 2.0 cm3_{, (b) 2.4× 2.4 × 2.4 cm}3_{, (c) 2.8× 2.8 × 2.8 cm}3_{and (d) 3.2× 3.2 × 3.2 cm}3_{. The blood}

perfusion rate is set to 1 kg m−3s−1.

It is found that the heating time strongly depends not only on the heating size, but also on the blood perfusion rate. For a blood perfusion rate 10 kg m−3s−1, the treatable PTV is decreased to 2.4× 2.4 × 2.4 cm3_{. This is because the temperature cannot be built up to the required} level due to strong energy dissipation caused by high blood perfusion rate. Figure10shows temperature responses at different blood perfusion rates for the 3.2 cm cube at z= 9.4 cm (front-end of the target region). It can be clearly seen that increasing the blood perfusion rate can decrease the slope of temperature rise, and it fails to reach the required temperature level at 10 kg m−3s−1. The thermal dose conformability is also compared and summarized in table1. The conformal index values in different cases are all above 76% in all treatable cases.

3.5. Comparison between different switching frequencies

The switching frequency for the sequentially scanned foci also plays an important role in the scanning-type heating strategy, especially in large-volume heating (Hynynen et al 1986).

2 2.4 2.8 3.2 0 100 200 300 400 500

Target volume width (cm)

Heating time (s)

10 5 Wb: 1 kg m-3s-1

Figure 9. Heating times for different target volume sizes at blood perfusion rates of 1, 5 and 10 kg m−3s−1. 0 100 200 300 400 500 38 40 42 44 46 48 50 52 54 Time (s) Temperature ( C) 10 5 Wb: 1 kg m-3s-1

Figure 10. Temperature elevations for different blood perfusion rates 1, 5 and 10 kg m−3s−1. The PTV is set to 3.2× 3.2 × 3.2 cm3with the centre located at z= 11 cm.

To investigate the role of switching frequency in determining treatment quality, results for switching frequencies of 10, 20 and 40 Hz are compared. Figure11(a) shows temperature responses under different switching frequencies at z= 9.4 cm and z = 11 cm for the 3.2 × 3.2× 3.2 cm3_{cube. Temperature fluctuation is more apparent at low switching frequencies,} particularly at high temperatures. From figure11(b), it is found that, with the same average temperature level, a larger temperature fluctuation can induce a higher thermal dose value. The phenomenon is caused by the exponential increment of thermal dose under different temperatures (see (6)). This temperature fluctuation caused by low switching frequency also affects the treatment conformability. In figure 12 it is shown that, under different blood perfusion rates, the optimal conformal index values at 20 and 40 Hz are similar, but worsen apparently at 10 Hz, especially when the heating size is smaller than 2.8× 2.8 × 2.8 cm3_.

(a) Temperature (C) 10 10 101 103

Thermal Dose (min)

Time (sec) (b) Time (sec) 0 100 200 300 400 500 40 45 50 55 60 0 100 200 300 400 500 -1 -3 105 240 40 20 SW:10Hz 40 20 SW:10Hz

Figure 11. (a) Temperature and (b) thermal dose responses at z= 9.4 cm (front-end of PTV) and z= 11 cm (centre of PTV) for the 3.2 × 3.2 × 3.2 cm3_{cube under different switching frequencies}

10, 20 and 40 Hz. 2 2.4 2.8 3.2 60 70 80 90 100

Target volume width (cm)

Conformal index (%) 2 2.4 2.8 3.2 60 70 80 90 100 Conformal index (%) (a)

Target volume width (cm) (b) w b = 1 kg m-3s-1 w b = 5 kg m-3s-1 40 20 SW:10Hz 40 20 SW:10Hz

Figure 12. Conformal index values for different target volumes and switching frequencies with the blood perfusion rate set to (a) 1 kg m−3s−1and (b) 5 kg m−3s−1.

The 20 Hz switching frequency is a good approximation for performing sequential heating on different heating sizes for various blood perfusion conditions.

4. Discussion

4.1. Key advantages of proposed system and heating strategy

This study investigates the feasibility of combining power optimization and acoustic window enlargement to continuously heat a large tumour and greatly shorten the treatment time for a 1D ultrasound phased array system. This synthesis exhibits two main advantages. First, tumours up to 3.2× 3.2 × 3.2 cm3can be heated continuously, while maintaining high thermal dose conformability in the designated target region. The heating time is shortened to only several minutes, since no cooling is applied during scanning. Compared to the several-hour

treatment time of existing ultrasound heating strategy, this could be a significant improvement in treatment time reduction (Fan and Hynynen 1996,Daum and Hynynen 1999). The second advantage is that, during scanning, only two 1D phased arrays are driven at one time. Hence, a switching device can switch the 12 1D phased arrays, and an 80-channel power amplifier is enough to actuate the 960 total elements of the 12 1D phased arrays. This largely reduces the complexity of the system and total hardware cost.

4.2. Approximation in the optimization algorithm

While scanning the target region, the power weightings of the scanned foci are adjusted to minimize the incident energy and prevent near-field heating. In the proposed algorithm, the power weightings of scanned foci are not altered directly to optimize the thermal dose, but an averaged power deposition is altered iteratively instead. Then, the power weightings are found by a constrained least-squares optimization method to approach the averaged power deposition. This strategy greatly reduces the number of thermal dose calculations and greatly shortens the total calculation time. However, the actual averaged power deposition of these weighted scanned foci is not identical to the given deposition, and the solution only approximates the optimum. An improved method could regard the algorithm solution as a initial point, which could be used to further fine-tune the weightings.

4.3. Problems in the acoustic window enlargement

Near-field heating can also be reduced by increasing the geometrical gain through enlarging the transducer aperture. In practice, however, the size of the acoustic window may be limited, as a bone may block the incident path of the ultrasonic energy. The ultrasonic energy is then significantly absorbed by the bone, causing severe ultrasound scattering and painful ‘hot spots’ (Fessenden et al 1984). This limits the treatable domain of this system to body portions with large acoustic windows, such as the abdomen or breast. Further, if the target volume is small, fewer strips of 1D arrays are enough to form the thermal lesions. This increases the treatable domain of our heating system.

4.4. Problem in attenuation change during heating

The attenuation coefficient within tissue changes during heating and exhibits a lesion shift effect towards skin when the attenuation is dramatically changed (over 200%) (Kolios et al 1999). Previous study has also shown that, at 65◦C, the ultrasound attenuation coefficient has a 100% maximum increment in muscle and about a 90% increment in liver and kidney (Damianou et al 1997). It has also been found that the attenuation is exposure-time dependent, with a 15% maximum increment within 4 min of exposure at 65◦C. The changes of attenuation are smaller for lower temperatures (Gertener et al 1997,Worthington and Sherar 2001). In the simulation results, only one case approaches 66◦C maximum temperature within 4 min of exposure. All others are under 61◦C. Hence, the constant-attenuation model is still close to the practical attenuation-changed conditions for most cases.

4.5. Considerations of ultrasound cavitation during treatment

The problem of cavitation is one of the most concerning issues in ultrasound thermal therapy. Cavities formed by high energy levels can cause unpredictable high temperatures within the treatment region and induce undesired lesion shape and formation, which should be avoided

during the thermal treatment. Hynynen (1991) proposed that the threshold of the acoustic pressure to generate transient cavitation in dog’s thigh muscle is described by

P = 0.6 + 5.3 × f (7)

under 1 s of exposure time, where P is the pressure magnitude (MPa) and f is the frequency (MHz). The validity of this formula has been verified through many experimental results. For example, Watkin et al (1996) delivered a peak intensity of 2500 W cm−2by using a 1.69 MHz transducer and found no cavitations, where the threshold intensity predicted by Hynynen is 2700 W cm−2. Other than the effects of the frequency and the intensity, previous studies also suggested that the increment of the exposure time can also reduce the cavitation threshold and should be considered (Atchley et al 1988,Sanghvi et al 1996). In our heating scenario, the entire PTV is scanned in a cycling pattern until a pre-determined heating time is reached. The cycling time between two exposures on the same position ranges from 8 to 25 s (calculated by total foci number/switching frequency), and the total exposure time on the same position is less than 1 s for all cases (calculated by total heating time/total assigned foci number) according to table1. For 800 kHz of driving frequency, the predicted threshold of cavitations is about 650 W cm−2with 1 s exposure time. Considering the fact that our maximum power is limited to 500 W cm−2, it is reasonable to assume that cavitations will not be a problem in our case.

4.6. Practical considerations

There are also three practical issues to be discussed: asymmetrical geometry of PTV, PTV larger than 3.2× 3.2 × 3.2 cm3_{, and effect of body movement during treatment. In the} first issue, the proposed power optimization algorithm is intrinsically x- and y-symmetric (the

z-axis can be asymmetric). To show the effect of treating a more realistic tumour under

this structure, one irregular tumour located in the female breast treatment case is considered (Leibel and Phillips 1998). In the x- and y-directions, the PTV is selected to cover the entire tumour region, and in the z-direction, the PTV is set to exactly the same as the tumour shape. The acoustical and thermal parameters are referred to previous studies by Lu (1996) (αbreast= 6 Np m−1, αwater= 0.02 Np m−1, k= 0.5 W m−1 ◦C−1, wb = 1.7 kg m−3 s−1, ct= cb= 4000 J kg−1 ◦C−1). The geometrical relationship between the treatment location and the phased array system as well as the treatment result are shown in figure13.

For the second issue, the maximal generated thermal lesions in this heating scenario are limited to 3.2× 3.2 × 3.2 cm3_{for blood perfusion rates 1 and 5 kg m}−3_s−1_{and 2.4× 2.4 ×} 2.4 cm3_{for 10 kg m}−3_s−1_{. In the case of a PTV larger than 3.2× 3.2 × 3.2 cm}3_{, it cannot be} completely covered by a single lesion with physical geometrical limitations. In this case, the entire PTV could be partitioned into several sub-regions and be treated sequentially. Cooling time is therefore necessary to be introduced between each sub-region’s heating to prevent the near-field heating problem. Although the total treatment time could be extended, the number of cooling periods has been largely reduced due to lesion volume enlargement in one single exposure. As a result, the treatment can be performed in a reasonable amount of time.

The third practical concern is about the PTV shifting during treatment caused by human motions, such as breathing. Such human motions can cause a mistargetting of the PTV and greatly decrease the treatment conformability. The average frequency of human breathing is about 15–20 times/min, i.e., 0.25–0.33 Hz. Normally, a mechanical position-tracking system can have a bandwidth of tens of Hz easily. Theoretically, an ultrasound-guided monitoring system can be built to track the position and compensate the position error, which is also an important future work for our group.

t_on= 140s CI = 93.8% EI = 12.6% Conformal Index = 80.2% PTV = 12.7 cm3 Necrosed vol. = 10.2 cm3 X-D (cm) Y-D (cm) Z-D (cm) X-D (cm) Y-D (cm) Z- D (cm) (a) (b) (c) TD = 240 min. TD = 240 min. Tumour Breast Acoustic window

Y

Z

0

X

Figure 13. One practical treatment case where the PTV is located in female breast: (a) geometrical relation between body portion and the heating system, and (b), (c) the PTVs and isosurfaces of TD= 240 min from different views. Arrows indicate the power-transmitted directions.

5. Conclusions

This study investigates the feasibility of combining power optimization and acoustic window enlargement to achieve a large volume heating in a single exposure without cooling. A multiple 1D phased array system was utilized to examine this, and the resulting heating time for treatments of large thermal lesions was significantly reduced to only several minutes with high treatment conformability. This is a great improvement in shortening the total treatment time of large-volume treatment. Although only preliminary simulation results are presented here, it is shown that the multiple 1D phased array system and the heating strategy have great

potential to treat a large tumour with less treatment time. Practical system construction and simulation verification are currently being investigated in our laboratory.

Appendix

A.1. Constraint least-squares optimization for determining foci weightings

Here we use a constraint least-squares optimization method to determine the power weighting for each focus to obtain a given absorbed power distribution (Golub and Van Loan 1996). First, we define the absorbed power density at the mth focal point induced by the nth generated focus as hmn=αρck 2 4π2 
S u ei ˆθn−(α+ik)(rm−r
) rm− r
ds 2 (A.1) where ˆθnis the phased difference on the driving elements to generate the ultrasonic foci at the predetermined position, m= 1 to M and n = 1 to N. Then, the relation between the driving power for N foci and the absorbed power density for M field positions can be represented as

 q...1 qM   =  h11... ... ...· · · h1N hM1 · · · hMN   ·  u...1 uN   (A.2)

orQ_M×1= H_M×NU_N×1where M N. Here, finding the appropriate power weightings u1to

uNthat obtain a desired absorbed power distribution ˆQ becomes a constrained least-squares

minimization problem:

minimizeH U − ˆQ2 subject to un
0. (A.3)

Since H U − ˆQ2 2=   q1− ˆq... 1 qM− ˆqM   2 2 =M m=1 (qm− ˆqm)2 (A.4)

(A.4) is then equivalent to minimizing the objective function with penalty as follows (Gottfried and Weisman 1973):

E =M m=1 (qm− ˆqm)2+K N  n=1 f (un) f (un) = 0 u_n
0 u2 n un< 0 (A.5) where K is the convergence factor and f (u_n) is the penalty function. Then, the excited weighting adaptation can be simply derived as follows:

un= −_∂u∂E n = − M  m=1 ∂E ∂qm ∂qm ∂un = −2 M  m=1 (qm− ˆqm)hmn+Kf
(un) n = 1 to N. (A.6)

Theoretically, all of the weightings unwill converge to a positive value only when K approaches

these weightings are adjusted to zero in the algorithm. Nevertheless, this solution closely approximates the actual minimum of the energy function (Gottfried and Weisman 1973). In our simulation cases, the omitted foci with their negative magnitude are all smaller than 2% of the maximum; thus, the influence of omitting negative weighting foci can still be small.

References

Atchley A A, Frizzell L A, Apfel R E, Holland C K, Madanshetty S and Roy R A 1988 Thresholds for cavitation produced in water by pulsed ultrasound Ultrasonics 26 280–5

Damianou C A and Hynynen K 1993 Focal spacing and near-field heating during pulsed high temperature ultrasound therapy Ultrasound Med. Biol. 19 777–87

Damianou C A, Hynynen K and Fan X 1995 Evaluation of accuray of a theoretical model for predicting the necrosed tissue volume during focused ultrasound surgery IEEE Trans. Ultrason. Ferroelect. Freq. Contr. 42 182–7 Damianou C A, Sanghvi N T, Fry F J and Maass-Moreno R 1997 Dependence of ultrasonic attenuation and absorption

in dog soft tissues on temperature and thermal dose J. Acoust. Soc. Am. 102 628–34

Daum D R and Hynynen K 1999 A 256-element ultrasonic phased array system for the treatment of larger volumes of deep seated tissue IEEE Trans. Ultrason. Ferroelect. Freq. Contr. 46 1254–68

Dewey W C 1994 Arrhenius relationships from the molecules and cell to the clinic Int. J. Hyperth. 10 457–83 Fan X and Hynynen K 1996 Ultrasound surgery using multiple sonications—treatment time considerations Ultrasound

Med. Biol. 22 471–82

Fessenden P, Lee E R, Anderson T L, Strohbehn J W, Meyer J L, Samulski T V and Marmor J B 1984 Experience with a multitransducer ultrasound system for localized hyperthermia of deep tissues IEEE Trans. Biomed. Eng. 31 126–35

Gertener M R, Wilson B C and Sherar M D 1997 Ultrasound properties of liver tissue during heating Ultrasound Med. Biol. 23 1395–403

Golub G H and Van Loan C F 1996 Matrix Computations (Baltimore, MD: Johns Hopkins University Press) pp 579–90

Gottfried B S and Weisman J 1973 Introduction to Optimization Theory (Englewood Cliffs, NJ: Prentice-Hall) Hutchinson E B and Hynynen K 1996 Intracavitary ultrasound phased arrays for noninvasive prostate surgery IEEE

Trans. Ultrason. Ferroelect. Freq. Contr. 43 1032–42

Hynynen K 1991 The threshold for thermally significant cavitation in dog’s thigh muscle in vivo Ultrasound Med. Biol. 17 157–69

Hynynen K, Roemer R, Moros E, Johnson C and Anhalt D 1986 The effect of scanning speed on temperature and equivalent thermal exposure distributions during ultrasound hyperthermia in vivo IEEE Trans. Ultrason. Ferroelect. Freq. Contr. 34 552–9

Kahn S B, Love R R, Sherman C and Chakravorty R 1983 Concepts in Cancer mediCInes (New York: Grune & Stratton)

Kolios M C, Sharer M D and Hunt J W 1999 Temperature dependent tissue properties and ultrasonic lesion formation Advances in Heat and Mass Transfer in Biotechnology ed E Scott (New York: ASME) pp 113–8

Leibel S A and Phillips T L 1998 Textbook of Radiation Oncology (Philadelphia, PA: Saunders)

Lin W L, Liang T C, Yen J Y, Liu H L and Chen Y Y 2001 Optimization of power deposition and a heating strategy for external ultrasound thermal therapy Med. Phys. 28 2172–81

Lin W L, Yen J Y, Chen Y Y, Jin K W and Shieh M J 1999 Relationship between acoustic aperture size and tumor conditions for external ultrasound hyperthermia Med. Phys. 26 818–24

Lu X-Q, Burdette E C, Bornstein B A, Hansen J L and Svensson G K 1996 Design of an ultrasonic therapy system for breast cancer treatment Int. J. Hyperth. 12 375–99

Malcolm A L and ter Haar G R 1996 Ablation of tissue volumes using high intensity focused ultrasound Ultrasound. Med. Biol. 22 659–69

Nyborg W L 1981 Heat generation by ultrasound in a relaxing medium J. Acoust. Soc. Am. 1 310–2

Ocheltree K B and Frizzell L A 1988 Determination of power deposition patterns for localized hyperthermia: a transient analysis Int. J. Hyperth. 4 281–96

O’Nell H T 1949 Theory of focusing radiators J. Acoust. Soc. Am. 21 516–26

Pearce J and Thomsen S 1995 Rate process analysis of thermal damage Optical–Thermal Response of Laser-Irradiated Tissue ed A J Welch and M J C van Gemert (New York: Plenum) pp 561–606

Pennes H H 1948 Analysis of tissue and arterial blood temperatures in the resting human forearm J. Appl. Phys. 1 93–122

Relkaitis G A, Ravindran A and Ragsdell K M 1983 Engineering Optimization: Methods and Applications (New York: Wiley)

Sanghvi N T, Fry F J, Bihrle R, Foster R S, Phillips M H, Syrus J, Zaitsev A V and Hennige C W 1996 Noninvasive surgery of prostate tissue by high-intensity focused ultrasound IEEE Trans. Ultrason. Ferroelect. Freq. Contr. 43 1099–110

Sapareto S A and Dewey W C 1984 Thermal dose determination in cancer therapy Int. J. Radiat. Oncol. Biol. Phys. 10 787–800

Saw C B and Suntharalingam N 1990 Quantitative assessment of interstitial implants Int. J. Rad. Oncol. Biol. Phys. 20 135–9

Wan H, Aarsvold J, O’Donnell M and Cain C A 1999 Thermal dose optimization for ultrasound tissue ablation IEEE Trans. Ultrason. Ferroelect. Freq. Contr. 46 913–28

Wan H, VanBaren P, Ebbini E S and Cain C A 1996 Ultrasound surgery: comparison of strategies using phased array systems IEEE Trans. Ultrason. Ferroelect. Freq. Contr. 43 1085–98

Watkin N A, ter Haar G R and Rivens I 1996 The intensity dependence of the site of maximal energy deposition in focused ultrasound surgery Ultrasound Med. Biol. 22 483–91

Worthington A E and Sherar M D 2001 Changes in ultrasound properties of porcine kidney tissue during heating Ultrasound Med. Biol. 27 673–82

Wu X and Sherar M 2002 Theoretical evaluation for moderately focused spherical transducers and multi-focus acoustic lens/transducer systems for ultrasound thermal therapy Phys. Med. Biol. 47 1603–21