Graded Mesh Change of variable in the vessel to allow for resolution of the boundary layer

1. A Fast, Accurate Algorithm Enabling EfficientSolution of a Drug Delivery ProblemCatherine E. Beni, Oscar P. BrunoApplied and Computational MathematicsCalifornia Institute of Technology Introduction Algorithm Numerical Results • 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. • 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: • Greatly disparate time-scales • Extremely steep boundary layers • Occurrence of very small diffusion coefficients • Graded Mesh • Change of variable in the vessel to allow for resolution of the boundary layer • Alternating Directions Implicit (ADI) • Discretize in time: • Split into operators on Cn+1 and Cn • and factor into differential operators in x and y • 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 • Use the following change of unknown • 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 Tissue Membrane Vessel Framework • The VMT convection-diffusion problem: • The parameters D, vx, and vy vary in each layer • VMT geometry: • 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 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 Future work • 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)

2. A Noise-tolerant Fejér-based modified-FBP Reconstruction Algorithm (Fejér-mFBP) for Positron Emission Tomography C.E. Beni, O.P. BrunoApplied and Computational MathematicsCalifornia Institute of Technology Introduction Modified-FBP Reconstructed Images • 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. • Both methods suffer in the presence of noise: • Direct methods amplify noise and show a dramatic loss of information in the reconstructed images • Iterative algorithms are not guaranteed to converge • Goal: to design a fast, accurate reconstruction algorithm that does not degrade substantially in the presence of noise • Approximate the Radon transform with its Fourier series: • where ak and bk are the Fourier coefficients • Compute derivative of Hilbert transform: • Approximate CT and ST, the Hilbert transforms of cosine and sine respectively, as follows: • Combine with the derivative of the Hilbert transform • and integrate to obtain the modified-Filtered Back Projection (mFBP) algorithm • Reconstructions using 711 values of ½ , 200 values of µ, and 200 Fourier modes → Unrealistic noise, unrealistically sensitive device • Reconstructions using 100 values of ½ , 200 values of µ, and 200 Fourier modes → Unrealistic noise, realistically sensitive device Framework • Reconstructions using 711 values of ½ , 200 values of µ, and 200 Fourier modes with 18.5% noise present → Realistic noise, unrealistically sensitive device • Radon Transform: • Geometry • Inverse Radon Transform: • h(½,µ)is the Hilbert transform of the Radon transform • Shepp-Logan phantom Fejér-mFBP • Fejér series of a given function: • Reconstructions using 100 values of ½ , 200 values of µ, and 200 Fourier modes with 18.5% noise present → Realistic noise, realistically sensitive device • By approximating the Radon transform with a Fejér series instead, we obtain the Fejér-mFBP algorithm: Original FBP Numerical Results • Compute the Hilbert transform and its derivative as follows: • Integrate to obtain the inverse Radon transform Conclusions • Both algorithms were implemented in C++ • 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)