Improving conformal tumour heating by adaptively removing control points from waveform diversity beamforming calculations: A simulation study

Abstract

Waveform diversity is a phased array beamforming strategy that determines an optimal sequence of excitation signals to maximise power at specified tumour control points while simultaneously minimising power delivered to sensitive normal tissues. Waveform diversity is combined with mode scanning, a deterministic excitation signal synthesis algorithm, and an adaptive control point removal algorithm in an effort to achieve higher, more uniform tumour temperatures. Simulations were evaluated for a 1444 element spherical section ultrasound phased array that delivers therapeutic heat to a 3 cm spherical tumour model located 12 cm from the array. By selectively deleting tumour control points, the tumour volume heated above 42°C increased from 2.28 cm³ to 11.22 cm³. At the expense of a slight increase in the normal tissue volume heated above the target temperature of 42°C, the size of the tumour volume heated above 42°C after tumour points were deleted was almost five times larger than the size of the original heated tumour volume. Several other configurations were also simulated, and the largest heated tumour volumes, subject to a 43°C peak temperature constraint, were achieved when the tumour control points were located along the back edge of the tumour and laterally around the tumour periphery. The simulated power depositions obtained from the results of the adaptive control point removal algorithm, when optimised for waveform diversity combined with mode scanning, consistently increased the penetration depth and the size of the heated tumour volume while increasing the heated normal tissue volume by a small amount.

• phased array beamforming
• waveform diversity
• mode scanning
• spot scanning
• therapeutic ultrasound

Introduction

Therapeutic ultrasound has great potential as a non-invasive modality for breast cancer treatments. Ultrasound provides an effective means of delivering hyperthermia, achieving targeted drug delivery, and thermally ablating cancer cells. Therapeutic ultrasound has a much greater penetration depth than other non-invasive modalities, facilitating non-invasive treatments in much deeper tumours. Therapeutic ultrasound is also advantageous for accessing tumours in the breast due to the absence of bone or air obstructions between the applicator and the tumour target.

To exploit the advantages of therapeutic ultrasound, several different types of phased array and multiple element applicators have been developed for a variety of clinical applications. For example, ultrasound phased array devices have been established for the ablation of uterine fibroids [1], and phased arrays are under development for non-invasive treatments of brain tumours [2]. Other multiple element ultrasound applicators have also been developed, including the SONOTHERM 1000 by Labthermics Technologies [3], the SURLAS device for delivering hyperthermia and radiation therapy simultaneously [4], and sectored interstitial ultrasound applicators that provide spatial control over the power deposition [5]. Phased array ultrasound devices have also been constructed for thermal therapy in the breast [6] and in the prostate [7].

An important limitation of non-invasive ultrasound applicators, particularly electronically focused phased arrays, is the small size of the focal spot relative to the size of the tumour. For example, the typical lateral extent of lesions generated by non-invasive HIFU systems is in the order of 0.1–0.3 cm [8], whereas breast tumours that are candidates for HIFU therapy are between 1.5 and 2.4 cm in diameter [9]. Complete coverage of large tumours can require hundreds of focal spots, and spot scanning approaches cause problems with intervening tissue heating and/or dramatic increases in the time required to complete a treatment. To address this problem, better beamforming strategies are needed for therapeutic applications. At present, the treatment times required for complete tumour coverage with ablative therapies are in the order of several hours [10]. The treatment time can be reduced when a spiral pattern of lesions is generated with single focal spots [11]; however, the steady-state temperature distributions generated by spiral patterns also generate intervening tissue heating. For hyperthermia or heat-modulated drug delivery applications, the objective is to maintain a therapeutic temperature of 42°C or 43°C throughout the entire tumour volume for 60 min, and spot scanning consistently encounters problems with relatively small heated tumour volumes and/or intervening tissue heating. Although the size of the heated volume is increased with split focusing [12], the problems with intervening tissue heating and long treatment times remain unsolved with this approach. Multiple focusing [13], [14] can simultaneously deliver power in several different locations, thereby reducing peak power values and the potential for cavitation. Multiple focusing can reduce intervening tissue heating and treatment times somewhat, but several issues still remain. As demonstrated in [15], waveform diversity has the potential to solve many of these problems by continuously cycling through a small number of focal patterns during a 60-min hyperthermia treatment. However, there is still considerable room for improvement with respect to the implementation of the waveform diversity algorithm.

This paper demonstrates that waveform diversity combined with mode scanning and an approach that adaptively modifies the distribution of control points within the tumour achieves higher, more uniform temperatures in larger volumes for hyperthermia than spot scanning or the waveform diversity approach that was previously evaluated in [15] while simultaneously minimising the effects of intervening tissue heating. Waveform diversity, mode scanning, and analytical bioheat transfer concepts are reviewed, followed by a description of the waveform diversity approach and the deterministic synthesis of excitation signals in [15], a description of the adaptive control point removal algorithm, and a summary of methods for computing the pressure and temperature fields. The results show that the 42°C isothermal surface achieves much better coverage of the 3 cm diameter spherical tumour volume when the deterministic waveform diversity approach is combined with mode scanning and the adaptive control point removal algorithm, where the simulation results are computed for a large spherical ultrasound phased array.

Motivation for improved beamforming

Waveform diversity

Waveform diversity beamforming optimises the covariance matrix R for a set of excitation signals applied to a phased array in an effort to maximise the output power in some locations and minimise the power delivered to others. As opposed to approaches that optimise the individual input signals directly, waveform diversity instead initially determines the covariance matrix R, which specifies the desired properties of a sequence of array excitations. Then, the input signals are extracted from R through an independent algorithm, possibly subject to additional constraints [16]. Traditionally, the generating function utilises a sequence of independent and identically distributed (i.i.d.) random variables, where many samples (in the order of 2¹⁵ to 2¹⁶ for the example in [15]) are required to achieve convergence within 1% of the entries in the covariance matrix R. However, when the covariance matrix is calculated with the deterministic algorithm in [15] that exploits the rank deficiency of R, only a small number of beam patterns (usually less than 6) are required to achieve the same accuracy. After the excitation signals are synthesised with the stochastic approach or the deterministic signal synthesis algorithm, the phased array is excited with continuous wave signals, and then the driving system switches between these signals at a sufficiently high rate to maintain steady-state temperatures at therapeutic levels.

By exploiting diverse excitation signals that produce a sequence of unique beam patterns, waveform diversity achieves several advantages over split focus and standard multiple focus techniques. One advantage is that when the input signals are synthesised with the deterministic algorithm in [15], sensitive normal tissues are spared while tumour coverage is maximised with only a small number of beam patterns. In contrast, traditional split focus and multiple focus approaches reduce the number of beam patterns somewhat relative to spot scanning without offering adequate protection for normal tissues. The multiple focus approach of Ebbini and Cain [13], [14] is effective for a small number of focal spots, but as the number of foci increases, hot spots appear in unpredictable locations. The Ebbini and Cain approach [13], [14] also permits the placement of nulls in sensitive normal tissues; however, only a limited number of these can be inserted without causing hot spots elsewhere. Furthermore, the optimal combination of tumour and normal tissue control points in the context of several unique beam patterns produced by a sequence of excitation signals remains an unsolved problem for the Ebbini and Cain approach [13], [14]. Another advantage of waveform diversity beamforming is that the optimisation of the covariance matrix R is a convex problem that is readily solved with available software such as SeDuMi [17], SDPT3 [18], and CVX [19]. Once the locations of the tumour and normal tissue control points are determined, these software packages calculate the optimal R matrix, and then a sequence of excitation signals are computed via the deterministic signal synthesis algorithm in [15].

Mode scanning

One of the main disadvantages of the standard waveform diversity approach is the large amount of computer memory required for the optimisation of the covariance matrix. This is especially important for waveform diversity calculations with large ultrasound phased arrays that have thousands of elements. For the array and control point geometry evaluated in Zeng et al. [15], if no simplifying assumptions are incorporated into the computer model, the memory required for calculations of the covariance matrix with SeDuMi is about 40 GB. Most modern desktop computers have less than 40 GB of RAM, so the optimisation routines will then spend considerable time swapping variables stored in memory to the hard drive, and the resulting calculations, which should be completed in a few hours, can instead take days or weeks. To address this problem and to exploit other advantages of symmetric beam patterns, waveform diversity is combined with mode scanning. By employing mode scanning, the number of phased array elements and control points in these calculations is reduced by a factor of 4. This in turn decreases the amount of computer memory by roughly a factor of 4³. Thus, a problem that requires tens of GB of RAM is converted into a much smaller problem that instead uses only a few hundred MB of RAM, and the optimisation of the covariance matrix is completed in about an hour.

The symmetric multiple focus patterns produced by mode scanning provide another advantage for tumour heating by cancelling the pressure fields along the two planes of symmetry of the array. This cancellation tends to shift energy away from the central axis, where the two planes of symmetry intersect and where the problems with intervening tissue heating tend to occur. In the context of these waveform diversity calculations, mode scanning is implemented by connecting the excitation signals transmitted to the phased array in symmetric groups of four elements. To achieve cancellation in the pressure fields throughout each plane of symmetry, 180° phase shifts (relative to the baseline signal transmitted to the group) are applied to the elements in quadrants II and IV, and no phase shift is applied to the elements in quadrants I and III. As a result, the pressure fields generated by each group of four elements cancel in the x = 0 and y = 0 planes defined in Figure 1, and the phase of the baseline signal transmitted to each group of four elements is adjusted according to the specifications of the beamforming algorithm. Thus, waveform diversity combined with mode scanning optimises the distribution of the energy deposited outside of the x = 0 and y = 0 planes.

Figure 1. Phased array and spherical tumour models. The 1444 element spherical section array is centred at the origin, and the 3 cm diameter spherical tumour model is located 12 cm from the bottom of the array. A square grid of normal tissue control points in quadrant I is positioned between the array and the tumour model to reduce intervening tissue heating.

Optimal power depositions

Bioheat transfer theory suggests that, in the absence of intervening tissue heating, the optimal analytical power distribution within a spherical tumour volume is given by a continuum of delta functions deposited on the tumour edge and constant power values in the tumour interior [20]. With focused ultrasound, the practical implementation of this ideal mathematical model of the power distribution is given by a set of discrete foci distributed on or near the tumour edge, possibly accompanied by additional low-level foci throughout the tumour interior. When intervening tissue heating is included in the simulation model, the optimal solution obtained from numerical simulations of spot scanning is given by a hemispherical focal point distribution on the distal half of the spherical tumour target [21]. With this approach, intervening tissue heating provides coverage for the tumour interior and the proximal tumour edge. These tumour heating models, which are obtained from theoretical analysis and numerical simulations, suggest that the waveform diversity results presented in Zeng et al. [15] can be further improved if the peak power deposition values are concentrated near the edge of the spherical tumour model. In Zeng et al. [15], tumour control points were distributed throughout the tumour, and approximately equal weightings were assigned to all tumour control points. By removing focal points from the tumour interior, intervening tissue heating problems should be reduced and conformal heating should be improved. Combining the information obtained from the bioheat transfer model with waveform diversity and mode scanning yields several possibilities for improved phased array beamforming in the context of thermal therapy, and one such algorithm that yields effective results is presented below.

Methods

Waveform diversity combined with mode scanning

For an ultrasound phased array subdivided into four symmetric quadrants where the even-numbered quadrants are excited 180° out of phase, the propagation vector defined for mode scanning [22] is given by

In Equation 1, , , , and represent the propagation vectors for the array elements in quadrants I, II, III, and IV, respectively, and the tilde indicates that these vectors are computed in quadrants, incorporating the symmetry required for mode scanning. The entries in these propagation vectors are ordered such that the array elements corresponding to the i^th entry in each vector are symmetric with respect to the x = 0 and y = 0 planes. For example, the first entry in each of these vectors is evaluated for array elements with centres at (x₁, y₁, z₁), (−x₁, y₁, z₁), (−x₁, −y₁, z₁), and (x₁, −y₁, z₁), respectively, where x₁ > 0 and y₁ > 0. Each propagation vector is evaluated at a single observation point r in quadrant I, where the accumulated power depositions are maximised for tumour control points and minimised for normal tissue control points through waveform diversity optimisation. Likewise, the same symmetry defined for , , , and is defined for , , , and such thatdetermines the driving signals for the entire array. The pressures at four symmetric control points are then given bywhich indicates that the pressures at the control points in quadrants II and IV are also 180° out of phase and that these pressures are determined by the product of the condensed propagation vector and the array excitation for elements in quadrant I. In Equation 3, r represents the point (x, y, z) in quadrant I at which the array is focused for mode scanning calculations, and the points at which the pressures are evaluated follow the same symmetric arrangement defined above for the array element coordinates in the propagation vectors.

The total power contributed by N focal patterns at a single point in space r is given bywhere α (Nepers/m) is the absorption coefficient, ρ (kg/m³) is the density of the medium, and c (m/s) is the speed of sound. In Equation 4, n represents the individual focal pattern or the excitation number, * represents the conjugate of a scalar, and *t represents the conjugate transpose of a vector or a matrix. All of the tumour and normal tissue control points are located in quadrant I, and is the covariance matrix for all N excitations when the locations of the control points and the geometry of the array elements satisfy the symmetries required for mode scanning. The propagation vector is a 1 by M row vector, the time-harmonic excitation vector is an M by 1 column vector, and the covariance matrix is an M by M Hermitian matrix, where M is the number of array elements in quadrant I. The array, which is centred at the origin, contains M elements in each quadrant and 4 M elements overall. In each of the equations above, the pressures and normal surface velocities are RMS values determined for individual beam patterns or array excitations, respectively. To generate a desired power deposition, the phased array rapidly switches between N different excitations such that the steady-state 3D output is approximated by the superposition of the N beam patterns computed with waveform diversity optimisation.

To maximise the power delivered to the tumour control points while minimising the power deposited at the normal tissue control points, the following constraints are defined for waveform diversity optimisation:where h̃ in Equation 5 is shorthand notation for h̃_{plane(x, y)} from Equation 1.

In Equation 5, μ represents a single normal tissue control point in quadrant I, ν represents a single tumour control point in quadrant I, and is a representative tumour control point that is arbitrarily selected from the tumour control points in quadrant I. The set of all normal tissue control points is represented by , and the set of all tumour control points is represented by . The first constraint requires that the power at the representative tumour control point is greater than the power deposited at every normal tissue control point (given that t > 0). The second constraint guarantees that the power at every tumour control point is greater than or equal to the power at the representative tumour control point . The third constraint defines as a positive semidefinite Hermitian matrix. The fourth constraint sets the total input power equal to the constant γ [16]. The optimisation problem defined in Equation 5 is then solved with SeDuMi [17].

By sequentially assigning each column of , where is the Hermitian square root of , to a complex array excitation , the deterministic signal synthesis problem is directly solved. Since is rank deficient, the number of distinct excitation signals is reduced with singular value decomposition (SVD). The SVD of is defined by , where and are M × N matrices with orthonormal columns, and is a diagonal N × N matrix that contains the singular values of . The reduced solution is given bywhere is the matrix square root of the diagonal matrix and each column of contains a distinct array excitation . As in Zeng et al. [15], all singular values less than 1/100 of the largest singular value in are discarded, and the number of singular values that remain determine the value of N. Once the N excitations are established in Equation 6, the corresponding array excitations in the other three quadrants are obtained through Equation 2.

Pressure and temperature simulations

Using the fast near-field method [23], the individual pressures generated by each element at each control point are calculated within a pressure input plane and then stored [24]. The total input pressure plane is then obtained from the superposition of the pressure planes generated by the individual elements weighted by the complex excitation signals determined from waveform diversity calculations. Then, the 3-D pressure field is calculated in homogeneous media using the angular spectrum approach [25]. The power depositions for N beam patterns are then superposed to obtain the total power deposition Q.

The temperature response is computed with a finite difference implementation of the steady-state bio-heat transfer equation [20],In Equation 7, K is the thermal conductivity (W/m/°C), T is the tissue temperature increase relative to the baseline defined by the arterial blood temperature (°C), W_b is the blood perfusion rate (kg/m³/s), C_b is the specific heat of blood (J/kg/°C), and Q is the local power deposition (W/m³).

Adaptive tumour control point removal algorithm

Tumour control points are initially distributed uniformly within the spherical tumour volume. Specifically, tumour control points are placed within concentric rings in discrete planes such that all of these are strictly located within the tumour radius. For the 3 cm diameter spherical tumour volume shown in Figure 1, six planes are defined with 0.45 cm spacing between adjacent planes. Within quadrant I, 33 tumour control points and one additional representative tumour control point r₀ are distributed in rings within these six planes. A 12 × 12 square grid of normal tissue control points are also distributed within a single plane in quadrant I at z = 9.45 cm. The tumour control points indicate locations where energy deposition is desired, and the normal tissue control points are placed in a region where ultrasonic energy tends to otherwise accumulate between the tumour and the phased array. After the waveform diversity approach combined with mode scanning determines the optimal , the excitations are synthesised according to the deterministic algorithm described in Zeng et al. [15], the temperature response is computed for the excitations such that the overall peak temperature is equal to 43°C, and the volumes that achieve temperatures of 42°C in the spherical tumour model and the surrounding normal tissues are calculated and stored. The tumour temperature values are evaluated strictly within the spherical tumour model, whereas the normal tissue temperatures are evaluated at all points other than the tumour in the 15 cm × 15 cm × 15.9 cm computational volume.

The adaptive control point removal procedure then deletes the tumour control point at which the temperature is greatest, where the representative tumour control point r₀ is included in the list of candidate points for deletion. Waveform diversity optimisations combined with mode scanning are then repeated, the tumour and normal tissue volumes with 42°C temperatures are stored, and the tumour control point at which the temperature is greatest is again deleted. This process is repeated until a single tumour control point remains in quadrant I. This adaptive algorithm is motivated by the analytical bioheat transfer results described in Section II.C and by the numerical results demonstrating that larger, more uniformly heated tumour volumes are achieved when control points are distributed near the distal edge of the tumour.

Results

The phased array model in Figure 1 has 1444 square elements, where each transducer is 0.24 cm × 0.24 cm. The geometric focus of the array is located at z = 12 cm, and the opening angle in each direction is 60°. The centre of the array is the origin of the coordinate system. The spherical tumour model used for these simulations has a diameter of 3 cm and is centred at (0, 0, 12) cm. The computational grid is a 15 cm × 15 cm ×15.9 cm volume with −7.5 cm ≤ x ≤ 7.5 cm, −7.5 cm ≤ y ≤7.5 cm, and 4.05 cm ≤ z ≤ 19.95 cm. The acoustic and thermal properties in the following simulations are: ρ = 1000 kg/m³, c = 1500 m/s, α = 0.5 dB/cm/MHz, W_b = 5 kg/m³/s, K = 0.55 W/m/°C, and C_b = 4000 J/kg/°C. These parameter values are the same as those used in [15]. The parameter values are modelled as constants throughout the entire simulated hyperthermia treatment. The array is driven by a 1 MHz excitation signal, so the wavelength is λ = 0.15 cm. This is the same array structure that was modelled in [15], which was retained to facilitate comparisons with the results presented herein.

The results of the waveform diversity simulations before any tumour control points are removed are shown in Figures 2 and 3. For these simulation results, one representative tumour control point r₀ and 33 additional tumour control points are distributed within the spherical tumour volume in quadrant I. In addition, 144 normal tissue control points were uniformly spaced on a regular grid in quadrant I between the phased array and the tumour target. The results in Figures 2 and 3 are very similar to those shown in Figures 5b and 6 in Zeng et al. [15]. The main differences are 1) slightly different constraints are defined for waveform diversity optimisation, and 2) only the 42°C isothermal contours are highlighted here, whereas the 41°C and 42°C isothermal contours are highlighted in Zeng et al. [15]. The similarity between the waveform diversity results in Figures 2 and 3 and the corresponding results in Zeng et al. [15] suggest that setting a lower bound for the power and minimising t in Equation 5 achieves approximately the same result as setting the upper and lower limits for the power and minimising t in the corresponding expression from Zeng et al. [15].

Figure 2. Simulated temperature (°C) in the y = 0 plane calculated in response to the power deposition obtained with waveform diversity and mode scanning before tumour control points are removed. The centre of the tumour model is located at z = 12 cm, and the 33 initial tumour control points are uniformly distributed throughout the 3 cm diameter tumour in quadrant I. The overall peak temperature in all bioheat transfer simulations is 43°C.

Figure 3. The 3 cm diameter spherical tumour volume and the 42°C isothermal surface calculated in response to the power deposition obtained with waveform diversity and mode scanning, where the initial 33 tumour control points are distributed uniformly throughout the spherical tumour volume in quadrant I. The 42°C isothermal surface, which covers only 16% of the tumour volume, is completely encompassed by the 3 cm diameter spherical tumour model, and no normal tissues are heated above 42°C.

Figure 5. The volume of the normal tissue that exceeds 42°C plotted with respect to the number of tumour control points. When tumour control points are initially removed, the heated normal tissue volume remains small. As additional tumour control points are removed, the heated normal tissue volume increases until the maximum value is reached at 12 tumour control points.

Figure 6. Simulated temperature (°C) in the y = 0 plane calculated from the power deposition obtained with waveform diversity and mode scanning using the 13 optimal tumour control points in quadrant I. This result corresponds to the largest value of the tumour volume heated above 42°C in Figure 4.

Figure 2 shows the temperature response in the y = 0 plane, and Figure 3 shows the 42°C isothermal surface. In Figures 2 and 3, the percentage of the tumour volume heated to 42°C is 16% (2.28 cm³ out of 14.14 cm³), and none of the normal tissue volume reaches 42°C. The temperature response in Figure 2, which shows isothermal contours in the y = 0 plane evaluated at 1°C increments while indicating the external tumour boundary with a thicker solid line, indicates that the waveform diversity result for 33 tumour control points plus one representative control point r₀ avoids problems with intervening tissue heating (with respect to the 42°C threshold) and that tumour coverage is incomplete. The 42°C isothermal contour in Figure 3, which is encapsulated by the spherical tumour volume, confirms this result. The isothermal contour in Figure 3 is also important for this analysis because intervening tissue heating, if present, is not always evident in individual temperature cross-sections, especially when mode scanning cancels pressures throughout the y = 0 plane and x = 0 plane. For the results shown in Figures 2 and 3, the deterministic algorithm for synthesising the excitation signals outlined in Zeng et al. [15] identified three singular values in so N = 3 focal patterns were defined for these calculations, which is the same number of focal patterns computed for the result shown in Zeng et al. [15].

Figures 4 and 5 show how the heated spherical tumour and normal tissue volumes, respectively, change as the adaptive algorithm deletes tumour control points. Figure 4 indicates the volume of the tumour model that exceeds 42°C, and Figure 5 describes the normal tissue volume that exceeds 42°C, where these volumes are plotted with respect to the number of control points remaining. The vertical axis in both plots ranges from 0 to 14 cm³, where the upper limit is the approximate volume of the spherical tumour model. In Figures 4 and 5, the values corresponding to the results shown in Figures 2 and 3 for 33 tumour control points plus the representative tumour control point r₀ are located on the far right of each plot.

Figure 4. The volume of the tumour that exceeds 42°C plotted with respect to the number of tumour control points. The heated tumour volume obtained with the initial tumour control point distribution is indicated on the far right. As tumour control points are removed, the heated tumour volume increases until the maximum size of the heated tumour volume is achieved with 13 tumour control points.

As control points are removed one by one and the waveform diversity simulations are repeated, the horizontal axes of Figures 4 and 5 are traversed from right to left. In Figure 4, the tumour volume heated above 42°C increases almost monotonically until the maximum value occurs at 13 tumour control points. In the vicinity of this maximum tumour volume, the curve in Figure 4 is relatively flat, and the tumour volume exceeds 10 cm³ (71% of the tumour volume) when the number of tumour control points is in the 12–20 range. When fewer than eight tumour control points are defined, the tumour volume heated above 42°C tends to decrease. Finally, the minimum value for the heated tumour volume is achieved when a single tumour control point is defined.

Figure 5 shows the change in the normal tissue volume heated above 42°C as tumour control points are removed according to the adaptive algorithm described above. Initially, there is very little change in the heated normal tissue volume as tumour control points are removed. Then, the heated normal tissue volume increases somewhat as the number of tumour control points drops below 25. A local minimum in the heated normal tissue volume is evident at 16 tumour control points, where the heated normal tissue volume is 2 cm³. The largest overall values for the volume of normal tissues heated above 42°C occur when 13 to 19 tumour control points are defined. The normal tissue volume heated above 42°C is generally smaller outside this range.

Based on the results shown in Figures 4 and 5, several different heating strategies can be defined. One such strategy emphasises the need to limit normal tissue heating over all other considerations. In this case, the number of tumour control points is restricted to the 25–33 range. Over this range, the volume of heated intervening normal tissues remains small. As the number of tumour control points is reduced to 25, the tumour coverage more than doubles (relative to the results shown in Figures 2 and 3 for 33 tumour control points), and only a small increase in normal tissue heating is observed. Throughout this range, the heated normal tissue volume remains below 0.522 cm³. Another strategy maximises the heated tumour volume at the expense of increased normal tissue heating. This approach, when applied to the results shown in Figures 4 and 5, restricts the number of tumour control points to the very largest heated tumour volumes, which correspond to the 13–18 range. If normal tissue heating is considered, then the result with 16 tumour control points is selected. If the size of the heated tumour volume is the overriding consideration, then the result obtained with 13 tumour control points is chosen.

The results obtained with 13 tumour control points in quadrant I (in addition to the representative tumour point at r₀) are shown in Figures 6 and 7. Figure 6 shows the temperature contour plot in the y = 0 plane, and Figure 7 contains the 42°C isothermal surface. These figures demonstrate the maximum heated tumour volume achieved with the adaptive control point removal algorithm combined with waveform diversity and mode scanning for the array and tumour geometry depicted in Figure 1. The tumour volume heated above 42°C in Figures 6 and 7 is 11.22 cm³ (79.35% of the 14.14 cm³ tumour volume), whereas the corresponding normal tissue volume heated above 42°C is 3.203 cm³. This result is achieved with N = 3 beam patterns. A comparison between the results shown in Figures 2 and 3 and in Figures 6 and 7 confirms that the tumour coverage achieved with the 13 optimised tumour control points is much greater than the tumour coverage with the 33 initial tumour control points. Although all of these results are computed with waveform diversity combined with mode scanning, adaptively removing control points increases the size of the heated tumour volume by 8.94 cm³ or 392% relative to the result obtained with the initial 33 control points. This improvement is achieved at the expense of an increase in normal tissue heating from 0 cm³ to 3.203 cm³. If the result for 16 tumour control points is instead selected, the normal tissue heating would then be reduced by about 38% (from 3.203 cm³ to 2 cm³) with only a small reduction (of about 4%) in the heated tumour volume relative to the results obtained with 13 tumour control points.

Figure 7. The 3 cm diameter tumour volume and the 42°C isothermal surface calculated in response to the power deposition obtained with waveform diversity and mode scanning, where the 13 optimal tumour control points are distributed on the back edge of the tumour and laterally about the tumour periphery in quadrant I. The 42°C isothermal surface covers 79% of the 3 cm diameter tumour, and normal tissue heating is observed along the near and far edges of the tumour relative to the centre of the array.

The 42°C isothermal contour obtained with 13 tumour control points in quadrant I is roughly conformal to the outline of the spherical tumour volume in Figure 6, and the 42°C isothermal surface covers most of the spherical tumour volume in Figure 7. In Figure 6, some normal tissue heating occurs on-axis beyond the tumour, and the under-heated region within the tumour model is near the tumour periphery at z = 12 cm. In Figure 7, some additional intervening tissue heating is observed in regions that are slightly off-axis and adjacent to the tumour. The near edge of the tumour is closely aligned with the 42°C isothermal contour in the cross-section depicted in Figure 6. Elsewhere, the 42°C isothermal contour is representative of the temperature distribution indicated by the 42°C isothermal contour in Figure 7. The result in Figures 6 and 7, when compared to the result in Figures 2 and 3, show that the tumour penetration depth increases when the optimal tumor control points are used in waveform diversity calculations.

Figure 8 shows the distribution of tumour control points within quadrant I. The solid square markers (▪) indicate the locations of the 13 optimal tumour control points that produced the result in Figures 6 and 7, where the tumour volume heated above 42°C reached the maximum value shown in Figure 4. The solid circle marker (•) represents the location of the representative tumour control point r₀. The solid triangle markers (▾) indicate the tumour control points that were removed from the initial distribution of 33 tumour control points. In Figure 8, the tumour control points that are closer to the top of the figure are near the back edge of the tumour, and the tumour control points that are closer to the bottom of the figure are near the front face of the tumour. These control point locations are evaluated relative to the location of the ultrasound phased array in Figure 1. The tumour control points are arranged along lines with equal x and y coordinates, where the z coordinate varies along each line with a constant increment. Figure 8 indicates that all of the tumour control points closest to the array axis (i.e. x = y = 0) are deleted by the adaptive control point removal algorithm. In the next two lines of tumour control points that are closest to the central axis of the phased array, only the points that are close to the far edge of the tumour are retained and the remaining points are deleted. In the other four lines of control points, which are closest to the tumour periphery, only one control point in each line is deleted, and the other three points are retained. Thus, all of the tumour control points that are retained are either located near the back edge of the tumour or laterally near the tumour periphery, and all other points are deleted.

Figure 8. Distribution of tumour control points in quadrant I within the 3 cm diameter spherical tumour model. The points that were removed from the initial control point distribution are indicated by a solid triangle (▾), the 13 optimal points that were retained are represented by a solid square (▪), and the representative point r₀ is indicated by a solid circle (•). To maximise 42°C tumour coverage, the tumour control points are located on the far edge of the tumour and laterally around the periphery of the spherical tumour model.

Discussion

Phased array model

Large ultrasound phased arrays are presently under development for external therapeutic applications, as demonstrated by the 1372 element phased array described by Song and Hynynen [2]. These large arrays are motivated by the need to generate higher intensities at the focus and to spread the energy generated by the phased array across a wider aperture. With these large array designs, undesirable side-effects are reduced in intervening tissues. Similar goals motivated the design of the array aperture shown in Figure 1, which is much larger than the 3 cm diameter spherical tumour model. If a much smaller aperture is used, efforts to continuously deliver approximately uniform temperatures for hyperthermia throughout a large, deep-seated volume are guaranteed to cause significant intervening tissue heating with or without optimal beamforming. As observed in Figures 2, 3, 6, and 7, the thermal contributions from the grating lobes generated by this dense 38 × 38 element array are negligible, which indicates that this design sufficiently samples the array aperture.

Comparison to previous results

In Zeng et al. [15], simulation results were evaluated for the same array geometry, the same tumour shape and size, the same tissue control point distribution, and the same initial tumour control point distribution. The main differences are the modified constraints in Equation 5, the temperature threshold for tumour and normal tissue evaluations, and the adaptive tumour control point removal algorithm. In Zeng et al. [15], the results were evaluated for a 4°C temperature increase (relative to a core temperature of 37°C), whereas here, a 5°C temperature increase was evaluated. In each case, the results are normalised such that the peak temperature is 43°C. The results in [15] showed that 56.5% of the tumour volume achieved temperatures of 41°C or higher, and the normal tissue volume heated above 41°C was 2 cm³. For the same volume of normal tissue heated above the threshold temperature, the adaptive control point removal algorithm achieves a higher temperature (42°C versus 41°C) in a larger tumour volume (10.61 cm³ or 75% tumour coverage with 16 tumour control points versus 7.99 cm³ with 33 tumour control points). If the input powers for the simulations evaluated herein are scaled by 4/5 (the ratio of the temperature increases evaluated with these two waveform diversity approaches) and the temperature target is then set to 41°C for both approaches, the heated tumour volume also increases from 7.99 cm³ to 10.61 cm³ (i.e. by 32.8%) for the same normal tissue volume due to the linearity of the bioheat transfer equation. In addition, the results shown here confirm that higher, more uniform temperatures are achieved when tumour control points are concentrated near the back edge of the tumour as in [21].

Comparison to spot scanning

The adaptive control point removal algorithm was also applied to single focus spot scanning (not shown). The initial focal point distribution in quadrant I was the same as that shown in Figure 8, and the focal points in the other three quadrants were symmetrically distributed with respect to the x = 0 and y = 0 planes. A total of 136 focal points were therefore distributed within the 3 cm diameter tumour model with 34 foci in each quadrant. All focal points were weighted equally. The total power deposition was calculated, and then the temperature response was computed. The focal point in quadrant I at which the temperature was highest was identified and deleted, along with the corresponding focal points in the other three quadrants that demonstrate symmetry with respect to the x = 0 and y = 0 planes. The power deposition and temperature response were recalculated for the updated focal point distribution, the focal point in quadrant I with the highest temperature was identified, and then the corresponding group of four symmetric focal points was deleted. This procedure was repeated until a single group of four symmetric focal points remained. After each iteration, the tumour volume and the normal tissue volume exceeding 42°C was evaluated and stored. The maximum tumour volume that exceeded 42°C was obtained with 20 focal points in quadrant I, so a total of 80 focal points were symmetrically distributed throughout the four quadrants. The tumour volume that exceeded 42°C with this focal point distribution was 4.576 cm³ or 32.367%, and the normal tissue volume that exceeded 42°C was 2.9084 cm³. More than three quarters of the tumour control point distributions computed for waveform diversity combined with mode scanning that were evaluated in Figure 4 produced tumour volumes above 42°C that are larger than this result, which is the largest tumour volume heated above 42°C with this set of focal locations using single focus spot scanning. In addition, only four of the tumour control point distributions computed for waveform diversity combined with mode scanning and evaluated in Figure 5 generated normal tissue volumes above 42°C that are larger than this result. Furthermore, more than half of the volumes shown in Figures 4 and 5 correspond to results that simultaneously achieve larger tumour volumes and smaller normal tissue volumes that are heated above 42°C, where several of these waveform diversity results double the heated volume while maintaining normal tissue volumes that are less than or equal to those produced by spot scanning. Thus, the results obtained with waveform diversity combined with mode scanning when tumour control points are selectively deleted are significantly better than the best steady-state temperature distributions obtained with spot scanning.

Heating a 2 cm diameter spherical tumour model

The adaptive control point removal algorithm, when applied to waveform diversity combined with mode scanning, is also effective for heating other tumour sizes. For example, a 2 cm diameter spherical tumour model (not shown) was populated with 11 tumour control points in quadrant I (plus one additional representative tumour point r₀), the power deposition was optimised for waveform diversity combined with mode scanning, and the temperature response was computed. The tumour control points were deleted one by one according to the adaptive removal algorithm, then the tumour and normal tissue volumes heated above 42°C were recorded. The optimal distribution associated with the largest tumour volume heated above 42°C consisted of seven tumour control points. The initial distribution with 11 tumour control points heated 0.964 cm³ of the 4.1888 cm³ (23%) tumour volume above 42°C, and the corresponding normal tissue volume heated above 42°C was 0.07425 cm³ with N = 2 beam patterns. The optimal distribution with seven tumour control points delivered temperatures of 42°C or greater to 2.496 cm³ or 59.6% of the tumour volume with N = 1 beam pattern. The normal tissue volume heated above 42°C with seven tumour control points was 0.38644 cm³. Similar to the optimal tumour control point distribution shown in Figure 8 for a 3 cm diameter tumour, all of the tumour control points closest to the central axis of the array were deleted, and the seven optimal tumour control points were all distributed laterally about the tumour periphery. Comparisons with the spot scanning results (also not shown) were similar to the comparisons described above for a 3 cm diameter tumour.

The heated tumour and normal tissue volumes were also evaluated for a 2 cm diameter tumour volume that was more densely populated with tumour control points. In quadrant I, 41 tumour control points were uniformly distributed throughout the 2 cm tumour volume, and the adaptive tumour control point algorithm was evaluated for results obtained with waveform diversity beamforming and mode scanning. The tumour volume heated above 42°C with 41 tumour control points generated by N = 2 beam patterns was 0.69 cm³, which corresponds to 16.4% tumour coverage. The normal tissue volume heated above 42°C with 41 tumour control points was 0 cm³. The maximum tumour coverage was reached with six tumour control points generated by N = 2 beam patterns, which heated 2.414 cm³ or 57.6% of the tumour above 42°C. The normal tissue volume heated above 42°C was 0.1485 cm³. Thus, with a higher initial density of tumour control points, the tumour volume and the normal tissue volume heated above 42°C with waveform diversity and mode scanning decrease slightly relative to the quantities obtained with the less dense tumour control point distribution. However, the absolute change in the heated tumour and normal tissue volumes is very small. One interpretation of this result suggests that, within a certain range of values, the beamforming approach evaluated here is relatively insensitive to the tumour control point spacings. Especially large tumour control point spacings are expected to underheat large portions of the tumour volume, and very small tumour control point spacings are not expected to achieve any significant benefit considering the excessive amount of computer time required to evaluate the adaptive control point removal algorithm with a large number of tumour control points.

Computation times

The total time required for a single iteration of the adaptive control point removal algorithm is approximately 65 min. Thus, calculating the results shown in Figures 4 and 5 takes just over 35 h on a computer with a quad core processor and 8 GB of RAM running the Linux operating system. Since the FNM computations are performed before the waveform diversity simulations, each iteration takes 63 min 42s for the waveform diversity calculation, 2.7s to superpose the input pressure planes [24] that were previously calculated with the fast nearfield method [23], 40.6s for the angular spectrum simulation, and 31.9s to calculate the bioheat transfer response. Clearly, the most time consuming procedure is the waveform diversity calculation. The pressure calculations take much less time due to the time-efficient and numerically accurate combination of the fast nearfield method and the angular spectrum approach.

An approach that is expected to reduce the total computation time instead selectively places all of the initial tumour control points near the outer edge of the tumour volume. This approach is suggested by previous numerical results [21] and by the optimal tumour control point distributions evaluated herein, which consistently delete all of the interior tumour control points and the majority of the tumour control points closest to the array. After the initial tumour control point distribution is defined, the adaptive control point removal algorithm would then determine the combination of these that heats the largest tumour volume above 42°C (or obtains the solution to another objective function based on the tumour and normal tissue volumes heated above the target temperature). This approach will reduce the number of initial tumour control points, thereby decreasing the total computation time by reducing the number of iterations required.

Better tumour coverage is also expected as more tumour control points are added in underheated tumour regions. The potential downside of adding more tumour control points in these locations is that normal tissue heating will most likely increase. This suggests that methods for determining the initial distribution of tumour control points are suboptimal. The adaptive control point removal algorithm obtains the optimal result for a given initial distribution, but even better results are anticipated with optimised initial control point distributions.

Reducing the number of elements in the ultrasound phased array will also greatly reduce the computation time. For example, waveform diversity is also ideal for beamforming with 1D phased arrays with 64 or 128 elements. Tumour and normal tissue control points would then be placed in a single plane, and the adaptive control point removal algorithm would determine the combination of tumour control points that heats the largest tumour volume to the target temperature (possibly subject to restrictions on the heated normal tissue volume that are incorporated into the objective function).

Non-monotonic variations in the heated tumour volume

In Figures 4 and 5, the change in the heated tumour volume is non-monotonic over a certain range of values away from the peak. Analysis of the individual power depositions and temperature responses suggests that this behaviour occurs when certain tumour control points are deleted. In other words, certain tumour control points contribute more strongly to non-uniform tumour heating, and when these tumour control points are removed, a sudden change in the heated tumour volume is observed. Likewise, the normal tissue volume also changes rapidly when certain tumour control points are removed. This suggests that when certain tumour control points are removed, normal tissue heating either increases or decreases depending on whether the deleted control point hinders or facilitates, respectively, the reduction of power deposited in regions where normal tissue heating tends to occur.

Temperature evaluations

The target temperature of 42°C selected for these simulations is intended to satisfy the treatment goal of at least 41–41.5°C in greater than 90% of the measured tumour volume during 60-min hyperthermia treatments in patients with locally advanced breast cancer [26]. This target temperature also achieves targeted drug release with temperature-sensitive liposomes designed for mild hyperthermia [27]. Thus, the 42°C target temperature in the tumour volume is motivated by clinical hyperthermia temperatures achieved in breast cancer patients, and this value is also above the phase transition temperature of thermally-sensitive liposomes designed for mild hyperthermia. The peak overall temperature of 43°C enforces relatively uniform tumour temperatures, which is one of the objectives of the waveform diversity beamforming algorithm evaluated here. Maintaining a constant peak temperature across all simulations also normalises the simulation results so that larger heated tumour volumes are not achieved through excessive normal tissue heating. Finally, the normal tissue threshold of 42°C enforces conformal temperature delivery within the tumour target. Smaller temperature thresholds for evaluations of normal tissue heating generate excessively large numbers for the calculated volume, even when the 42°C isothermal surface closely conforms to the tumour boundary. If the temperature threshold for evaluating normal tissue heating is much larger than the tumour target temperature, then the 42°C isothermal surface generated by the optimal beamforming result can conform poorly to the outline of the target volume. Since conformal heat delivery is preferable for these waveform diversity beamforming simulations, 42°C was also selected as the temperature threshold in normal tissues.

The threshold values determined for these comparisons of scanning methods are steady-state temperatures. To achieve steady-state during a 60-min hyperthermia treatment, both single focus spot scanning and waveform diversity methods rapidly switch between patterns, thereby reducing transient effects [28]. Evaluating the transient temperatures produced by waveform diversity beam patterns and incorporating these values into thermal dose calculations is a topic for future research. Evaluating the effect of changing blood perfusion during a mild hyperthermia treatment and determining how this impacts the results of waveform diversity beamforming is also a topic for future research. If the tumour control points are distributed near the edge of the tumour and if perfusion changes uniformly throughout the tumour volume, feedback compensation based on measured temperature values is expected to achieve acceptable results by adaptively scaling the power levels associated with each multiple focus pattern, where similar approaches are also applicable to single focus spot scanning. Better results are expected when waveform diversity beamforming results are optimised for a set of predicted perfusion distributions during treatment planning. The treatment planning results would then be incorporated into the feedback loop, and different sets of multiple focus patterns will be selected by the feedback controller when blood perfusion increases or decreases as indicated by changes in measured temperatures. If blood perfusion changes by different amounts in different regions of the tumour, then much faster optimisation procedures will be needed for waveform diversity beamforming with unequal tumour control point weightings, more complicated treatment planning algorithms will be required, or methods that combine waveform diversity beamforming results with single focus spot scanning will be needed. As the objective functions for waveform diversity evolve into simpler expressions, the time required for optimisation of the excitation signals for waveform diversity beamforming is expected to decrease by several orders of magnitude. Ideally, these calculations will eventually be completed in real-time. If a combined approach is required, the waveform diversity beamforming result will establish the baseline power deposition, and then additional focal spots will be generated wherever additional power is needed. This hybrid scanning approach will therefore employ waveform diversity beamforming to reduce intervening tissue heating and spot scanning to adapt to changing blood perfusion values, where the optimal combination of waveform diversity and spot scanning for mild hyperthermia is expected to simultaneously achieve better conformal heating and lower temperatures in intervening tissues than methods that strictly employ spot scanning.

The effect of inhomogeneities

In a clinical setting, inhomogeneities change the power depositions generated by invasive and non-invasive ultrasound devices by (1) spatially shifting the location of the focal peak(s), by (2) diminishing the peak acoustic power output, and by (3) increasing the sidelobe levels. These are problems for both fixed focus and electronically steered ultrasound devices. Fixed focus devices solve the first problem by mechanically shifting the transducer until the focal spot is directed at the desired target, the second problem is addressed by increasing the electrical input power until the required acoustic power output is achieved, and the third problem cannot be resolved with a fixed focus system. Likewise, phased array devices solve the first problem by electronically and/or mechanically shifting the focal spot, the second problem is addressed by phase aberration correction and/or increased electrical input powers, and the third problem is addressed with phase aberration correction. All of these solutions point to the need for non-invasive technologies that enable compensating feedback during clinical treatments. Two such non-invasive technologies include MRI thermometry and ultrasound-based methods for phase aberration correction. MRI thermometry measures the temperature distribution, and from this information, the location of the focus and the acoustic power output is extracted. Phase aberration correction compensates for speed of sound variations along the propagation path by determining the additional time delay or phase offset required for each element.

These solutions are equally applicable to single focus spot scanning and multiple focus beamforming with waveform diversity. If a multiple focus waveform diversity pattern is shifted relative to the desired target, then the applicator will be mechanically scanned and/or the appropriate linear phase terms will be applied to the excitation signals to compensate for this spatial shift. Likewise, if the acoustic power outputs are insufficient, then the electrical input powers will be scaled accordingly and/or phase aberration correction will compensate for diminished powers due to inhomogeneities. If sidelobe levels are increased or if the multiple focus pattern is distorted, then phase aberration correction will compensate for the phase errors introduced by the inhomogeneities. Each solution indicated here is also applicable to single focus spot scanning, and similar results are anticipated when these compensations are applied to both scanning methods.

Other applications

The adaptive tumour control point removal algorithm combined with waveform diversity and mode scanning, unlike most other beamforming approaches defined for thermal therapy with large ultrasound phased arrays, considers both power and temperature during optimisation. The excitation signals obtained from waveform diversity specify the optimal power deposition subject to the constraints in Equation 5, and the optimal temperature distribution is determined by the adaptive tumour control point removal algorithm. The desired result is selected from one of several possible objective functions defined relative to the tumour and normal tissue volumes that exceed the target temperature of 42°C in Figures 4 and 5, respectively. Thus, the power deposition and the temperature response are optimised in normal tissue and in the tumour volume, which facilitates more conformal heat delivery within the tumour target.

The waveform diversity beamforming approach described herein is also applicable to other tumour shapes, tumour sites, and phased array geometries. For example, ellipsoidal geometries and other four-fold symmetric shapes that are aligned with the coordinate axes are amenable to waveform diversity beamforming combined with mode scanning. Although the results shown here emphasise symmetric phased array and tumour geometries, the symmetry requirements imposed by mode scanning were primarily needed so that the problem considered could be solved with the available computer resources (namely, a quad core computer with 8 GB of RAM). Since four symmetric foci were removed at each iteration of the adaptive control point removal algorithm, mode scanning also reduced the overall simulation time. Other preliminary simulation results (not shown) suggest that waveform diversity and the adaptive control point removal algorithm are the main sources of improved temperature localisation and that the main advantages of mode scanning for these calculations are therefore the reductions in computer memory and computation time. Furthermore, if waveform diversity beamforming calculations are evaluated on a desktop computer with at least 40 GB of RAM, the symmetry imposed on the phased array and the tumour volume will no longer be necessary, thus facilitating waveform diversity beamforming without mode scanning. Although the computation time will increase accordingly, removing the symmetry constraint imposed by mode scanning will enable waveform diversity calculations for asymmetric phased array geometries and for asymmetric tumour volumes. Eliminating the symmetry constraint will also enable waveform diversity calculations in laterally shifted tumour volumes. Waveform diversity beamforming for sparse and random arrays can also be optimised with this approach. For symmetric and asymmetric tumour models, the initial tumour control point distribution will either be specified by uniformly depositing control points throughout the tumour or by preferentially placing tumour control points near the edge of the tumour. Normal tissue control points will be placed in sensitive normal tissues and in locations where excessive normal tissue heating is anticipated.

These beamforming algorithms are also applicable to ablative therapies with ultrasound phased arrays. Instead of ablating a tumour volume with a spot scanning approach with a waiting time between individual focal spots, the power, temperature, and thermal dose would be optimised for a collection of focal points. Larger regions containing multiple tumour control points would then be ablated simultaneously. The waveform diversity beam-forming approach, as applied to ablative therapies, is also expected to alleviate problems with intervening tissue heating, thereby enabling a potentially significant reduction in the total treatment time.

Conclusion

Waveform diversity and mode scanning were combined with an adaptive control point removal algorithm, and the simulation results showed that after the tumour control points identified by the algorithm were removed, much larger tumour volumes were heated above 42°C within a 3 cm diameter tumour volume. This was achieved at the expense of increased normal tissue heating. However, these results were significantly better than the results obtained using the initial tumour control point distribution for waveform diversity combined with mode scanning or the results obtained using the adaptive tumour control point removal algorithm with spot scanning. The optimal tumour control points extracted from the initial distribution defined for waveform diversity beamforming were consistently distributed along the back edge of the tumour and laterally about the tumour periphery. Interior tumour control points tend to overheat the near edge of the tumour, contributing to the non-uniformity of the tumour temperature distribution and to the heat delivered within intervening tissues. Simulation results also showed that the penetration depth, as demonstrated by the extent of the 42°C isothermal surface, was increased when the power deposition was computed with waveform diversity combined with mode scanning for the optimal tumour control points determined by the adaptive algorithm. These effects were also observed in simulation results obtained for a smaller 2 cm diameter tumour volume and a more densely sampled 2 cm diameter tumour volume. The deterministic excitation signal synthesis algorithm [15] determined that the optimal value for N was between one and three, so only a small number of multiple focus patterns were needed to heat each tumour. Approaches to reduce the overall computation time were suggested, and other applications, including ablative therapies, that are expected to benefit from calculations that employ waveform diversity and the adaptive control point removal algorithm, were identified. In each of these applications, as demonstrated here for simulations of hyperthermia delivered by a 1444 element ultrasound phased array, power depositions computed with waveform diversity for the optimal distribution of tumour control points obtained from the adaptive tumour control point removal algorithm are expected to achieve higher, more uniform temperatures throughout the tumour volume.

Declaration of interest: This work was supported in part by NSF Grant 0634786. The authors report no conflicts of interest. The authors alone are responsible for the content and writing of the paper.