Numerical Modelling of High-Intensity Focused Ultrasound (HIFU) Phased Arrays

1
Paul Godden, Gail ter Haar, Ian Rivens
Therapeutic Ultrasound Group Joint Department of
Physics The Institute of Cancer Research The
Royal Marsden NHS Trust
2
Outline of the Presentation

• Part I Physics Clinical Applications
• Overview of HIFU
• HIFU as a treatment modality
• Tumour ablation
• Acoustic cavitation

• Part II Mathematical Modelling - Phased Array
  Transducers Acoustic Beam Propagation
• Model equations and the parabolic
  approximation
• The extended angular spectrum approach
• Quadrature for highly-oscillatory integrals
• Finite volume methods for hyperbolic
  conservation laws

3
Part I High-Intensity Focused Ultrasound
Physics and Clinical Applications
4
Ultrasound - what is it? A pressure wave,
requires a dense medium
gt16 kHz (audible frequency)
5
Ultrasound as a Function of Acoustic Intensity
6
Ultrasound Propagation in Tissue
Attenuation Absorption Scatter
7
HIFU Treatment of Liver Tumours I
Lesion at focus
(12x3mm)
Skin
Target organ (e.g. liver)
Transducer
Tumour
Undamaged tissue in front of focus
8
HIFU Treatment of Liver Tumours II

• Volume ablation animation and MR image courtesy
  of Feng Wu

Kennedy JE (2005) Nature Reviews Cancer 18 March
321-327
9
Focussing of Ultrasound for Tumour Ablation

• An example one of the ICR transducers with 10
  equal
• area therapy strips

Phased arrays allow the focus to be
dynamically moved in a three-dimensional region
Constructed in a collaboration with Imasonic SA
(France)
10
Examples of Transducer Designs
Integrated therapy MR imaging device for
treatment of prostate tumours
Single-element 1.7 MHz extra-corporeal transducer
Single-element HIFU transducer with central
aperture for imaging transducer

• Frequencies from 0.5 to 4.0 MHz, with spatial
  peak intensities up to 20 kWcm-2
• Focal lengths from 3 to 15 cm
• Many different geometries
• Phased arrays allow multiple simultaneous foci
  and electronic steering of the focal region.

11
Why HIFU?

• Achieves the necessary temperature rise in the
  timescale required to ablate tumours

• Margin of ablation is accurate to the order of a
  few cell diameters

• Ablated tissue is re-absorbed by the body

• Non-invasive (extracorporeal transducers)

• No blood loss

• No therapeutic impasse and repeatable

unlike radiotherapy!
12
Trans-Rectal Combined Imaging Therapy Probe
Treatment of localised prostate cancer (image
courtesy of EDAP)
13
Cavitation
Associated with various fluid flow regimes
Increased energy deposition can lead to
mechanical (as opposed to thermal) damage
Hypothesised enhanced thermal transfer and,
consequently, tissue damage
However

• Acoustic cavitation inception physics currently
  unknown
• Even an order-of-magnitude estimate of bubble
  population not possible
• Modelling is therefore difficult!

14
Dynamics of Acoustic Cavitation in an Applied
Ultrasonic Field

• Cavitation bubbles act as both wave scatterers
  (useful for imaging) and sources (due to
  nonlinear oscillation)
• Individual bubbles may coalesce or collapse, as
  well as oscillate
• Previous modelling work has used fast multipole
  methods (for small populations of relatively
  large bubbles), some attempts at DNS-type
  analyses (simple ones!)
• Monte Carlo methods have also been attempted
  (with limited success)

15
Increased Thermal Effects due to Cavitation
Bubbles immediately
Bubbles after 0.5s
No bubbles
Passive cavitation (bubble) detection
16
Part 2 Mathematical Modelling of HIFU Propagation
17
Numerical Methods Simple Domain Decomposition
De-gassed water
Transducer (multi-element phased array)
Liver
Transducer
Linear region, close to the transducer E-ASA
method, FE methods
Linear region, close to the transducer E-ASA
method, FE methods
Region of significant nonlinear wave propagation
spectral, finite volume, finite difference methods
Region of significant nonlinear wave propagation
finite volume methods
18
The Angular Spectrum Approach Overview

• The Angular Spectrum Approach (ASA) is a useful
  method for computing the linear approximation to
  the acoustic field from a multi-element
  transducer array
• The ASA method was originally derived for planar
  transducers, but has been modified to accommodate
  focused arrays.
• In particular, the method is efficient for
  arrays with large numbers of transducer elements
  (some arrays have 1000 individual elements)

5 ring approximation
Array of 9 multi-width strips
Array of 18 strips
256 element array template
19
The ASA Method
Hydrophone positioning system
Transducer
Hydrophone
Oscilloscope
PC (via data acquisition interface)
Array drive circuitry
20
The ASA Method

• The basic ASA algorithm is
• - sample the pressure field in a 2D plane at a
  distance z from the transducer
• - take the DFT of the pressure field
• - multiply each pressure sample in the plane by
  a propagator term to advance the
• the sampled field to a distance z1 from the
  transducer
• - take the inverse DFT to acquire the pressure
  field in the plane z1

Assume harmonic pressure variation
Substitution into wave equation gives Helmholtz
form
Assume angular spectrum in FT form
Linear relation for pressure
Propagator term
21
The E-ASA Method

• Consideration of the contribution from each
  elemental curved surface in the array leads to a
  frequency-domain expression for the angular
  spectrum from a focused transducer array (Wu
  Stepinski, 1999)

22
Levin-type Quadrature for Highly-Oscillatory
Integrals

• A method for efficient and accurate solution of
  highly-oscillatory multidimensional integrals
• Original formulation (Levin, 1982) for problems
  involving a rectangular domain subject to a
  nonresonance condition and an absence of critical
  points in the domain of integration

No critical points in the domain
23
Levin-type Quadrature for Highly-Oscillatory
Integrals
24
Error Analysis for Univariate Levin-type
Quadrature
(Levin, 1997)
The second lemma deals with the error in the
collocation approximation to the non-oscillatory
solution.
25
E-ASA Comments

• Preliminary results indicate that 3rd-order
  expansion in bivariate monomials gives acceptable
  accuracy
• Consequently, for each transducer element a 9x9
  (dense) linear system needs to be solved
• Hence, for large transducer arrays (hundreds or
  thousands of elements), computational
  requirements are smaller than for Fourier
  transform methods
• Nonresonance condition holds since transducer
  has an aperture in the centre, thus eliminating
  the turning point from the computational domain

26
Approximate Methods for Ultrasound Propagation
The Khokhlov-Zabolotskaya-Kuznetsov Equation

• Assumes a parabolic approximation for the
  acoustic wave propagation.
• First consequence expression only valid up to
  30o off-axis.
• Second consequence the approximation is
  paraxial.

is the parabolic term,
Where
and
is the Khokhlov number a dimensionless ratio of
characteristic nonlinearity scale to diffraction
scale.
27
Numerical Solution of the KZK Equationthe
Standard Method

• All of the theoretical work in HIFU has been
  formulated according to engineering principles
• Models are simplified and feature assumed
  approximations
• Many aspects of the models, such as attenuation
  in tissue, are not empirically quantifiable at
  present!

• Standard current methods based upon Christopher
  Lee (1991) formulation
• - Models refraction and reflection (but not
  multiple reflections) in a region of multiple
  parallel fluid layers.
• - Limited to axisymmetric single-element
  sources.
• - Not limited by parabolic approximation, but
  not full-wave either.

28
The Method of Christopher Parker Overview
ALGORITHM (extended form of Christopher Lee,
1991)

• Take discrete Hankel transform (DHT) of planar
  pressure field

• Propagate over small distance, account for
  attenuation and diffraction of the linear field
  contribution

• Acquire a frequency domain solution to Burgers
  equation to account for nonlinear accretion
  depletion of harmonics

• Perform inverse DHT and update new pressure field

However, note that (c.f Korpel, 1980)
Solution to Burgers equation uses weakly
nonlinear theory and a Fourier series expansion
i.e. assume Burgers equation in the form
where
29
3D Numerical Solution for the KZK Equation
Implemented by Yang Cleveland (2005)

• Operator splitting fractional steps on each KZK
  term

- Diffraction
- Nonlinearity
- Thermoviscous absorption
- A finite number of (fluid) relaxation
processes
30
The Method of Christopher Parker Solutions for
an Axisymmetric Piston-Type Source
UNFORTUNATELY, the open-source Texas KZK code
generated some errors. Compiling under g77 14
errors, 9 warnings compiling under gcc gt50
errors compiling the C (??) code with g
too many to count
31
LeVeques Formulation of Finite Volume Methods

• Apply Godunovs method by solving the Riemann
  problem for each grid face per time step (coarse
  grid cells depicted for clarity)
• High-resolution corrections via Lax-Wendroff

• The Riemann problem is solved by Roes method
  (parametric integration via the parameter
  )

(Godden et. al., 2006)
32
The Riemann Problem

• The nonlinear wave propagation problem at each
  cell interface is posed as a Riemann problem a
  piecewise-constant jump discontinuity

Cell i
Cell i-1
33
LeVeques Formulation of Finite Volume Methods

• Formulation
• - cell-centred (Riemann problem formulation)
• - Godunovs method
• 
  
• - Lax-Wendroff for 2nd-order
• - Limiters for flux correction, e.g.
• - 3D via flux differencing methods
• (rather than fractional steps)
  LeVeque (1997)

34
Time Discretisation

• Dosimetry
• - Need to compute total thermal energy delivered
  to different tissue regions
• - Over continuous beam exposures, symplectic
  integrators necessary in order to correctly model
  conservation properties? (Assuming the Euler
  equations as governing equations)

35
Application 1 HIFU for Brachytherapy Treatment
Failure

• Treatment modality for localised prostate cancer
  insertion of radioactive seeds (e.g.
  iodine-125) into the prostate to produce
  continual localised low-level radiotherapeutic
  dose for approximately a year (brachytherapy)
• Cure rates gt80 at 10 years but in a well
  selected patient group with low-volume and
  non-aggressive disease
• Men who fail this treatment have very limited
  treatment options the only widely accepted
  treatment is hormone therapy which is
  non-curative - controls the disease for typically
  3-4years

• Important questions
• - is it possible to use HIFU as a salvage
  treatment for patients in which
  brachytherapy was not succesful?
• - if so, how do the seeds affect the
  HIFU field?

36
Application 1 HIFU for Brachytherapy Treatment
Failure

• Model the orientation of the seeds relative to
  the US imaging field
• Model the scatter, attenuation and reflection
  due to the presence of the seeds
• Possible to obtain solutions indicating the
  effect of the seeds upon steering of the focus?
• What is the magnitude of the change in heat flux
  due to the thermal conductivity of the seeds?
• Do the seeds act as preferential cavitation
  seeding regions? (Difficult to model!)

37
Application 2 Modelling Random Tissue
Inhomogenities

• Fogarty LeVeque (1999) modelled linear
  acoustic propagation through a randomly-varying
  inhomgeneous medium

• Model 1 Periodic Medium
• - material parameters are piecewise constant

38
Application 2 Modelling Random Tissue
Inhomogenities

• Possible to model nonlinear HIFU propagation
  through tissue using this method?
• A new limiter was derived to overcome numerical
  instabilities when standard limiters were used

39
Application 3 Modelling Large-Scale Tissue
Inhomogenities

• Large-scale inhomogeneities, such as blood
  vessels, may affect both beam propagation and
  thermal conductivity
• Accurate imaging of blood vessels is necessary
  to avoid haemorrhage during treatment
• Simple approximation define the vessel using
  flux function method in the finite volume mesh
• Fluid dynamics significant? (Lynch et. al. 1996)
• Simple conservation or advection models
  sufficent?

Angiogram of rat circulatory system
40
Conclusions Acknowledgements

• Professor David Levin (Department of
  Mathematical Sciences, University of Tel-Aviv)
• Dr Ping Wu (Department of Signals Systems,
  University of Uppsala)
• Professor Tadeusz Stepinski (Department of
  Signals Systems, University of Uppsala)
• Professor Victor Humphrey (ISVR, University of
  Southampton)
• Professor Ron Roy (Department of Aerospace and
  Mechanical Engineering, University of Boston)
• Dr Greg Clement (Brigham Womens Hospital,
  Harvard University)
• EPSRC (grant no EP/C511298)
• And..
• Hugh Morris, Alex Chapman, Jim McLaughlin,
  Chaturika Jayadewa, Yixin Ma, John
  Civale (Therapeutic Ultrasound Team, Institute of
  Cancer Research)

Thank
you!
41
Error Analysis for Univariate Levin-type
Quadrature
Introduce an m-vector of linearly independent
highly-oscillatory functions (oscillations may be
irregular)
Perform collocation separately on each
subinterval using monomial basis functions, and
sum the individual approximations
42
Error Analysis for Univariate Levin-type
Quadrature
Outline of the analysis, following Levin (1997).
We assume that the ODE may be written in the form
Assuming
These assumptions, together with Theorem 1, give
43
Error Analysis for Univariate Levin-type
Quadrature Lemma 1
Proof
Re-write the ODE as
Define a sequence of successive approximations to
the solution
44
Differentiating the previous expression (using
same assumptions as per Theorem 1)
From the expression
it follows that
45
which satisfies
with initial data set to zero and where
By standard ODE theory (Picards Existence
Theorem) the solution of
46
We have
and hence
with
Thus, we have
and
47
(No Transcript)
48
Error Analysis for Univariate Levin-type
Quadrature Lemma 2
Proof
49
this is a system for the vector of unknowns
we can re-write this in the form
where
50
(No Transcript)
51
From approximation theory we have the result