Transcript Style E 24 by 48

A Fast, Accurate Algorithm Enabling Efficient
Solution of a Drug Delivery Problem
Catherine E. Beni, Oscar P. Bruno
Applied and Computational Mathematics
California Institute of Technology
Introduction
The goal of magnetic drug delivery is to use magnetic fields to direct and
confine magnetically-responsive particles bound to therapeutic agents to
specific regions in a patient’s body-- thus allowing for focused treatment in an
area of interest.
Algorithm
 Graded Mesh
D = .0001, vy = .0001, Ren = .01
Change of variable in the vessel to allow for resolution of the boundary layer
Tissue
To design a method leading to confinement of the magnetically-responsive
particles to a particular region of the body, a predictive capability must be
used to evaluate the effects of external magnetic forces on the convection and
diffusion of magnetic particles through the bloodstream and in membranes
and tissue.
The numerical solution of the Vessel-Membrane-Tissue (VMT) convection
diffusion problem proposed by Grief and Richardson is highly challenging:
Numerical Results
Membrane
COMSOL
VMT Solver
Speed
36 hours
< 5 minutes
Memory Requirements
32 GB
32.7 MB
→ 432 times faster.
→ 1000 times reduction in memory requirements.
Vessel
Alternating Directions Implicit (ADI)
 Discretize in time:
 Split into operators on Cn+1 and Cn
Greatly disparate time-scales
D = .00001, vy = .00001, Ren = .001
Extremely steep boundary layers
Occurrence of very small diffusion coefficients
and factor into differential operators in x and y
COMSOL
VMT Solver
Speed
N/A
< 8 minutes
Memory Requirements
> Available 32Gb
98.3 Mb
Framework
The VMT convection-diffusion problem:
 The parameters D, vx, and vy vary in each layer
 VMT geometry:
 Solve the resulting scheme using Finite Difference methods
 On-and-Off Fluid Freezing Methodology
 Evolve algorithm until concentration in vessel reaches steady state
 “Freeze” concentration in vessel by not applying solver in vessel region
 Iterate twice using large time steps
 Reduce time step, unfreeze vessel, and iterate until concentration in vessel
reaches steady state
 Repeat
 Change of Unknown
Conclusions
 Developed a fast, efficient solver for a drug delivery problem
 432 times faster than commercial package COMSOL Multiphysics
 1000 times reduced memory requirements
 Allows for solution of previously intractable problems
 Use the following change of unknown
Future work
The VMT solver is based on a combination of
Use of a graded mesh to adequately resolve boundary layers
The Alternating Directions Implicit (ADI) method to overcome the
overwhelmingly restrictive CFL condition imposed by the fine spatial
discretization mentioned above
 An on-and-off fluid-freezing methodology that allows for efficient treatment
of the multiple time-scales that coexist in the problem (whose equilibria arise
through a complex balance of fluid-flow, magnetic-pull and diffusion effects)
 A change of unknown that enables evaluation of steady states in tissue and
membrane layers through a highly accelerated time-stepping procedure
to transform the differential operator in y in the membrane and tissue
regions into the Helmholtz equation to make use of a previously known fast
time-stepping method
 Finite Difference methods restrict us to a rectangular geometry
 Room for accuracy improvement
 These two problems will be fixed by solving the ODEs present in the ADI method
with the new Fourier Continuation-Alternating Directions (FC-AD) methodology
 See the talk by O.P. Bruno for more details
References
“The Behaviors of Ferro-Magnetic Nano-Particles in Blood Vessels under Applied
Magnetic Fields”, A. Nacev, C.E. Beni, O.P. Bruno, B. Shapiro (to be submitted)
A Noise-tolerant Fejér-based modified-FBP Reconstruction Algorithm (Fejér-mFBP)
for Positron Emission Tomography
C.E. Beni, O.P. Bruno
Applied and Computational Mathematics
California Institute of Technology
Introduction
 Images can be reconstructed from Positron Emission Tomography (PET) scanners
via two methods: Iterative methods and Direct methods
 Direct methods, such as the well-known Filtered Back Projection (FBP) algorithm
are fast, but reconstruct images that are low resolution.
 Iterative methods, such as ML-EM (Maximum Likelihood-Expectation
Maximization) and OSEM (Ordered Subset Expectation Maximization), are much
slower (each iteration takes the same amount of time as a full reconstruction using a
direct method and approximately 20-30 iterations are required), but provide high
quality reconstructions.
Modified-FBP
 Approximate the Radon transform with its Fourier series:
where ak and bk are the Fourier coefficients
MATLAB’s ‘iradon’
 Combine with the derivative of the Hilbert transform
 Shepp-Logan phantom
Fejér-mFBP
mFBP
Fejér-mFBP
 Reconstructions using 711 values of ½ , 200 values of µ, and 200 Fourier modes
with 18.5% noise present → Realistic noise, unrealistically sensitive device
MATLAB’s ‘iradon’
 Geometry
mFBP
 Approximate CT and ST, the Hilbert transforms of cosine and sine respectively,
as follows:
Framework
 Radon Transform:
MATLAB’s ‘iradon’
 Reconstructions using 100 values of ½ , 200 values of µ, and 200 Fourier modes
→ Unrealistic noise, realistically sensitive device
 Direct methods amplify noise and show a dramatic loss of information in the
reconstructed images
 Goal: to design a fast, accurate reconstruction algorithm that does not
degrade substantially in the presence of noise
 Reconstructions using 711 values of ½ , 200 values of µ, and 200 Fourier modes
→ Unrealistic noise, unrealistically sensitive device
Compute derivative of Hilbert transform:
 Both methods suffer in the presence of noise:
 Iterative algorithms are not guaranteed to converge
Reconstructed Images
mFBP
Fejér-mFBP
and integrate to obtain the modified-Filtered Back Projection (mFBP) algorithm
Fejér-mFBP
 Fejér series of a given function:
 Inverse Radon Transform:
 Reconstructions using 100 values of ½ , 200 values of µ, and 200 Fourier modes
with 18.5% noise present → Realistic noise, realistically sensitive device
MATLAB’s ‘iradon’
 h(½,µ) is the Hilbert transform of the Radon transform
mFBP
Fejér-mFBP
By approximating the Radon transform with a Fejér series instead, we obtain
the Fejér-mFBP algorithm:
Original FBP
 Compute the Hilbert transform and its derivative as follows:
Numerical Results
Conclusions
 Both algorithms were implemented in C++
 Integrate to obtain the inverse Radon transform
 Each reconstruction requires ~3.6 seconds, the same amount of time used by
MATLAB’s built-in ‘iradon’ function
 All reconstructions shown here are of the well-known Shepp-Logan phantom
generated using MATLAB’s built-in ‘phantom’ command
 Developed a new reconstruction algorithm that, in presence of noise, yields
iterative-solver-like quality at FBP computational costs.
References
 “A Noise-Tolerant Fejér-based modified-FBP Reconstruction Algorithm (Fejér-mFBP)
for Positron Emission Tomography”, C.E. Beni, O.P. Bruno (to be submitted)