Search

In this article, we present a multispecies reaction-advection-diffusion partial differential equation coupled with linear elasticity for modeling tumor growth. The model aims to capture the phenomenological features of glioblastoma multiforme observed in magnetic resonance imaging (MRI) scans. These include enhancing and necrotic tumor structures, brain edema and the so-called "mass effect", a term-of-art that refers to the deformation of brain tissue due to the presence of the tumor. The multispecies model accounts for proliferating, invasive and necrotic tumor cells as well as a simple model for nutrition consumption and tumor-induced brain edema. The coupling of the model with linear elasticity equations with variable coefficients allows us to capture the mechanical deformations due to the tumor growth on surrounding tissues. We present the overall formulation along with a novel operator-splitting scheme with components that include linearly-implicit preconditioned elliptic solvers, and a semi-Lagrangian method for advection. We also present results showing simulated MRI images which highlight the capability of our method to capture the overall structure of glioblastomas in MRIs.

We present a generative approach for simultaneously registering a probabilistic atlas of a healthy population to brain magnetic resonance (MR) scans showing glioma and segmenting the scans into tumor as well as healthy tissue labels. The proposed method is based on the expectation maximization (EM) algorithm that incorporates a glioma growth model for atlas seeding, a process which modifies the original atlas into one with tumor and edema adapted to best match a given set of patient's images. The modified atlas is registered into the patient space and utilized for estimating the posterior probabilities of various tissue labels. EM iteratively refines the estimates of the posterior probabilities of tissue labels, the deformation field and the tumor growth model parameters. Hence, in addition to segmentation, the proposed method results in atlas registration and a low-dimensional description of the patient scans through estimation of tumor model parameters. We validate the method by automatically segmenting 10 MR scans and comparing the results to those produced by clinical experts and two state-of-the-art methods. The resulting segmentations of tumor and edema outperform the results of the reference methods, and achieve a similar accuracy from a second human rater. We additionally apply the method to 122 patients scans and report the estimated tumor model parameters and their relations with segmentation and registration results. Based on the results from this patient population, we construct a statistical atlas of the glioma by inverting the estimated deformation fields to warp the tumor segmentations of patients scans into a common space.

This paper investigates the problem of atlas registration of brain images with gliomas. Multiparametric imaging modalities (T1, T1-CE, T2, and FLAIR) are first utilized for segmentations of different tissues, and to compute the posterior probability map (PBM) of membership to each tissue class, using supervised learning. Similar maps are generated in the initially normal atlas, by modeling the tumor growth, using reaction-diffusion equation. Deformable registration using a demons-like algorithm is used to register the patient images with the tumor bearing atlas. Joint estimation of the simulated tumor parameters (e.g., location, mass effect and degree of infiltration), and the spatial transformation is achieved by maximization of the log-likelihood of observation. An expectation-maximization algorithm is used in registration process to estimate the spatial transformation and other parameters related to tumor simulation are optimized through asynchronous parallel pattern search (APPSPACK). The proposed method has been evaluated on five simulated data sets created by statistically simulated deformations (SSD), and fifteen real multichannel glioma data sets. The performance has been evaluated both quantitatively and qualitatively, and the results have been compared to ORBIT, an alternative method solving a similar problem. The results show that our method outperforms ORBIT, and the warped templates have better similarity to patient images.

Intravascular ultrasound (IVUS) is a catheter-based medical imaging technique that produces cross-sectional images of blood vessels and is particularly useful for studying atherosclerosis. In this paper, we present a novel method for segmentation of the luminal border on IVUS images using the radio frequency (RF) raw signal based on a scattering model and an inversion scheme. The scattering model is based on a random distribution of point scatterers in the vessel. The per-scatterer signal uses a differential backscatter cross-section coefficient (DBC) that depends on the tissue type. Segmentation requires two inversions: a calibration inversion and a reconstruction inversion. In the calibration step, we use a single manually segmented frame and then solve an inverse problem to recover the DBC for the lumen and vessel wall (kappa(l) and kappa(w), respectively) and the width of the impulse signal theta. In the reconstruction step, we use the parameters from the calibration step to solve a new inverse problem: for each angle theta(i) of the IVUS data, we reconstruct the lumen-vessel wall interface. We evaluated our method using three 40MHz IVUS sequences by comparing with manual segmentations. Our preliminary results indicate that it is possible to segment the luminal border by solving an inverse problem using the IVUS RF raw signal with the scatterer model.

We present a framework for black-box and flexible simulation of soft tissue deformation for medical imaging and surgical planning applications. Our main motivation in the present work is to develop robust algorithms that allow batch processing for registration of brains with tumors to statistical atlases of normal brains and construction of brain tumor atlases. We describe a fully Eulerian formulation able to handle large deformations effortlessly, with a level-set-based approach for evolving fronts. We use a regular grid-fictitious domain method approach, in which we approximate coefficient discontinuities, distributed forces and boundary conditions. This approach circumvents the need for unstructured mesh generation, which is often a bottleneck in the modeling and simulation pipeline. Our framework employs penalty approaches to impose boundary conditions and uses a matrix-free implementation coupled with a multigrid-accelerated Krylov solver. The overall scheme results in a scalable method with minimal storage requirements and optimal algorithmic complexity. We illustrate the potential of our framework to simulate realistic brain tumor mass effects at reduced computational cost, for aiding the registration process towards the construction of brain tumor atlases.

In this article, we propose a framework for modeling glioma growth and the subsequent mechanical impact on the surrounding brain tissue (mass-effect) in a medical imaging context. Glioma growth is modeled via nonlinear reaction-advection-diffusion, with a two-way coupling with the underlying tissue elastic deformation. Tumor bulk and infiltration and subsequent mass-effects are not regarded separately, but captured by the model itself in the course of its evolution. Our formulation is fully Eulerian and naturally allows for updating the tumor diffusion coefficient following structural displacements caused by tumor growth/infiltration. We show that model parameters can be estimated via optimization based on imaging data, using efficient solution algorithms on regular grids. We test the model and the automatic optimization framework on real brain tumor data sets, achieving significant improvement in landmark prediction compared to a simplified purely mechanical approach.

BackgroundMRI characteristics of brain gliomas have been used to predict clinical outcome and molecular tumor characteristics. However, previously reported imaging biomarkers have not been sufficiently accurate or reproducible to enter routine clinical practice and often rely on relatively simple MRI measures. The current study leverages advanced image analysis and machine learning algorithms to identify complex and reproducible imaging patterns predictive of overall survival and molecular subtype in glioblastoma (GB).MethodsOne hundred five patients with GB were first used to extract approximately 60 diverse features from preoperative multiparametric MRIs. These imaging features were used by a machine learning algorithm to derive imaging predictors of patient survival and molecular subtype. Cross-validation ensured generalizability of these predictors to new patients. Subsequently, the predictors were evaluated in a prospective cohort of 29 new patients.ResultsSurvival curves yielded a hazard ratio of 10.64 for predicted long versus short survivors. The overall, 3-way (long/medium/short survival) accuracy in the prospective cohort approached 80%. Classification of patients into the 4 molecular subtypes of GB achieved 76% accuracy.ConclusionsBy employing machine learning techniques, we were able to demonstrate that imaging patterns are highly predictive of patient survival. Additionally, we found that GB subtypes have distinctive imaging phenotypes. These results reveal that when imaging markers related to infiltration, cell density, microvascularity, and blood-brain barrier compromise are integrated via advanced pattern analysis methods, they form very accurate predictive biomarkers. These predictive markers used solely preoperative images, hence they can significantly augment diagnosis and treatment of GB patients.

With this work we release CLAIRE, a distributed-memory implementation of an effective solver for constrained large deformation diffeomorphic image registration problems in three dimensions. We consider an optimal control formulation. We invert for a stationary velocity field that parameterizes the deformation map. Our solver is based on a globalized, preconditioned, inexact reduced space Gauss-Newton-Krylov scheme. We exploit state-of-the-art techniques in scientific computing to develop an effective solver that scales to thousands of distributed memory nodes on high-end clusters. We present the formulation, discuss algorithmic features, describe the software package, and introduce an improved preconditioner for the reduced space Hessian to speed up the convergence of our solver. We test registration performance on synthetic and real data. We demonstrate registration accuracy on several neuroimaging datasets. We compare the performance of our scheme against different flavors of the Demons algorithm for diffeomorphic image registration. We study convergence of our preconditioner and our overall algorithm. We report scalability results on state-of-the-art supercomputing platforms. We demonstrate that we can solve registration problems for clinically relevant data sizes in two to four minutes on a standard compute node with 20 cores, attaining excellent data fidelity. With the present work we achieve a speedup of (on average) 5 with a peak performance of up to 17 compared to our former work. [ABSTRACT FROM AUTHOR]

Kernel-based methods are a powerful tool in a variety of machine learning and computational statistics methods. A key bottleneck in these methods is computations involving the kernel matrix, which scales quadratically with the problem size. Previously, we introduced ASKIT (Approximate Skeletonization Kernel Independent Treecode), an efficient, scalable, kernel-independent method for approximately evaluating kernel matrix-vector products. ASKIT is based on a novel, randomized method for efficiently factoring off-diagonal blocks of the kernel matrix using approximate nearest neighbor information. In this context, ASKIT can be viewed as an algebraic fast multipole method for arbitrary dimensions. In this paper, we introduce our open-source implementation of ASKIT. Features of our ASKIT library include linear dependence on the input dimension of the data, the ability to approximate kernel functions with no prior information on the kernel, and scalability to tens of thousands of compute cores and data with billions of points or hundreds of dimensions. We also introduce some new extensions and improvements of ASKIT, included in our library. We introduce a new method for adaptively selecting approximation ranks and correctly partition the nearest neighbor information, both of which improve the performance of ASKIT over our previous implementation. We describe the ASKIT algorithm in detail, and collect and summarize our previous theoretical complexity and error bounds in one place. We present a brief selection of experimental results illustrating the accuracy and scalability of ASKIT. We then provide some details and guidance for users of ASKIT. [ABSTRACT FROM AUTHOR]

The solution of a constant-coefficient elliptic Partial Differential Equation (PDE) can be computed using an integral transform: A convolution with the fundamental solution of the PDE, also known as a volume potential. We present a Fast Multipole Method (FMM) for computing volume potentials and use them to construct spatially adaptive solvers for the Poisson, Stokes, and low-frequency Helmholtz problems. Conventional N-body methods apply to discrete particle interactions. With volume potentials, one replaces the sums with volume integrals. Particle N-body methods can be used to accelerate such integrals. but it is more efficient to develop a special FMM. In this article, we discuss the efficient implementation of such an FMM. We use high-order piecewise Chebyshev polynomials and an octree data structure to represent the input and output fields and enable spectrally accurate approximation of the near-field and the Kernel Independent FMM (KIFMM) for the far-field approximation. For distributed-memory parallelism, we use space-filling curves, locally essential trees, and a hypercube-like communication scheme developed previously in our group. We present new near and far interaction traversals that optimize cache usage. Also, unlike particle N-body codes, we need a 2:1 balanced tree to allow for precomputations. We present a fast scheme for 2:1 balancing. Finally, we use vectorization, including the AVX instruction set on the Intel Sandy Bridge architecture to get better than 50% of peak floating-point performance. We use task parallelism to employ the Xeon Phi on the Stampede platform at the Texas Advanced Computing Center (TACC). We achieve about 600gflop/s of double-precision performance on a single node. Our largest run on Stampede took 3.5s on 16K cores for a problem with 18e+9 unknowns for a highly nonuniform particle distribution (corresponding to an effective resolution exceeding 3e+23 unknowns since we used 23 levels in our octree). [ABSTRACT FROM AUTHOR]

We present a comparison of different multigrid approaches for the solution of systems arising from high-order continuous finite element discretizations of elliptic partial differential equations on complex geometries. We consider the pointwise Jacobi, the Chebyshev-accelerated Jacobi, and the symmetric successive over-relaxation smoothers, as well as elementwise block Jacobi smoothing. Three approaches for the multigrid hierarchy are compared: (1) high-order h-multigrid, which uses high-order interpolation and restriction between geometrically coarsened meshes; (2) p-multigrid, in which the polynomial order is reduced while the mesh remains unchanged, and the interpolation and restriction incorporate the different-order basis functions; and (3) a first-order approximation multigrid preconditioner constructed using the nodes of the high-order discretization. This latter approach is often combined with algebraic multigrid for the low-order operator and is attractive for high-order discretizations on unstructured meshes, where geometric coarsening is difficult. Based on a simple performance model, we compare the computational cost of the different approaches. Using scalar test problems in two and three dimensions with constant and varying coefficients, we compare the performance of the different multigrid approaches for polynomial orders up to 16. Overall, both h-multigrid and p-multigrid work well; the first-order approximation is less efficient. For constant coefficients, all smoothers work well. For variable coefficients, Chebyshev and symmetric successive over-relaxation smoothing outperform Jacobi smoothing. While all of the tested methods converge in a mesh-independent number of iterations, none of them behaves completely independent of the polynomial order. When multigrid is used as a preconditioner in a Krylov method, the iteration number decreases significantly compared with using multigrid as a solver. Copyright © 2015 John Wiley & Sons, Ltd. [ABSTRACT FROM AUTHOR]

Abstract: We propose an algorithm to compute an approximate singular value decomposition (SVD) of least-squares operators related to linearized inverse medium problems with multiple events. Such factorizations can be used to accelerate matrix-vector multiplications and to precondition iterative solvers. We describe the algorithm in the context of an inverse scattering problem for the low-frequency time-harmonic wave equation with broadband and multi-point illumination. This model finds many applications in science and engineering (e.g., seismic imaging, subsurface imaging, impedance tomography, non-destructive evaluation, and diffuse optical tomography). We consider small perturbations of the background medium and, by invoking the Born approximation, we obtain a linear least-squares problem. The scheme we describe in this paper constructs an approximate SVD of the Born operator (the operator in the linearized least-squares problem). The main feature of the method is that it can accelerate the application of the Born operator to a vector. If N ω is the number of illumination frequencies, N s the number of illumination locations, N d the number of detectors, and N the discretization size of the medium perturbation, a dense singular value decomposition of the Born operator requires operations. The application of the Born operator to a vector requires work, where μ(N) is the cost of solving a forward scattering problem. We propose an approximate SVD method that, under certain conditions, reduces these work estimates significantly. For example, the asymptotic cost of factorizing and applying the Born operator becomes . We provide numerical results that demonstrate the scalability of the method. [Copyright &y& Elsevier]

We present algorithms for the solution of a class of source identification problems for systems governed by elliptic partial differential equations (PDEs) on two-dimensional regular geometries. The state is the solution of the PDE, which is driven by an unknown source field. Given observations of the state, we seek to reconstruct the source field. We consider the cases of full domain observations and boundary observations. The problem is formulated as a least-squares PDE-constrained optimization problem. We use a reduced space approach in which we "invert" the associated Hessian using a preconditioned conjugate gradient (PCG) algorithm. Using standard Fourier analysis, we derive analytical solutions for the case in which the governing PDE has constant coefficients. Based on these solutions, we construct preconditioners that accelerate the convergence of PCG in the case of variable-coefficient elliptic PDE constraints. We performed numerical experiments to show the effectiveness of the pre-conditioners for variable coefficients with different contrasts and smoothness properties. We observed mesh-independent and β-independent convergence for different cases of the variable coefficients. The computational complexity of solving the source identification problem is O(N log N). The construction of the preconditioner costs O(N3/2), where N is the discretization size for the source. [ABSTRACT FROM AUTHOR]

We present a multigrid algorithm for the solution of source identification inverse problems constrained by variable-coefficient linear parabolic partial differential equations. We consider problems in which the inversion variable is a function of space only. We consider the case of L2 Tikhonov regularization. The convergence rate of our algorithm is mesh-independent—even in the case of no regularization. This feature makes the method algorithmically robust to the value of the regularization parameter, and thus useful for the cases in which we seek high-fidelity reconstructions. The inverse problem is formulated as a PDE-constrained optimization. We use a reduced-space approach in which we eliminate the state and adjoint variables, and we iterate in the inversion parameter space using conjugate gradients. We precondition the Hessian with a V-cycle multigrid scheme. The multigrid smoother is a two-step stationary iterative solver that inexactly inverts an approximate Hessian by iterating exclusively in the high-frequency subspace (using a high-pass filter). We analyze the performance of the scheme for the constant coefficient case with full observations; we analytically calculate the spectrum of the reduced Hessian and the smoothing factor for the multigrid scheme. The forward and adjoint problems are discretized using a backward-Euler finite-difference scheme. The overall complexity of our inversion algorithm is O(NtN + N log2 N), where N is the number of grid points in space and Nt is the number of time steps. We provide numerical experiments that demonstrate the effectiveness of the method for different diffusion coefficients and values of the regularization parameter. We also provide heuristics, and we conduct numerical experiments for the case with variable coefficients and partial observations. We observe the same complexity as in the constant-coefficient case. Finally, we examine the effectiveness of using the reduced-space solver as a preconditioner for a full-space solver. [ABSTRACT FROM AUTHOR]

State-of-the art algorithms for deformable image registration are based on the minimization of an image similarity functional that is regularized by adding a penalty term on the deformation map. The penalty function typically represents a smoothness regularization. In this article, we use a constrained optimization formulation in which the image similarity functional is coupled to a biophysical model. This formulation is pertinent when the data have been generated by imaging tissue that undergoes deformations due to an actual biophysical phenomenon. Such is the case of coregistering tumor-bearing brain images from the same individual. We present an approximate model that couples tumor growth with the mechanical deformations of the surrounding brain tissue. We consider primary brain tumors-in particular, gliomas. Glioma growth is modeled by a reaction-advection-diffusion PDE, with a two-way coupling with the underlying tissue elastic deformation. Tumor bulk, infiltration, and subsequent mass effects are not regarded separately but are captured by the model itself in the course of its evolution. Our formulation allows for updating the tumor diffusion coefficient following structural displacements caused by tumor growth/infiltration. Our forward problem implementation builds on the PETSc library of Argonne National Laboratory. Our reformulation results in a very small parameter space, and we use the derivative-free optimization library APPSPACK of Sandia National Laboratories. We test the forward model and the optimization framework by using landmark-based similarity functions and by applying it to brain tumor data from clinical and animal studies. State-of-the-art registration algorithms fail in such problems due to excessive deformations. We compare our results with previous work in our group, and we observed up to 50% improvement in landmark deformation prediction. We present preliminary validation results in which we were able to reconstruct deformation fields using four degrees of freedom. Our study demonstrates the validity of our formulation and points to the need for richer datasets and fast optimization algorithms. [ABSTRACT FROM AUTHOR]

In this article, we propose new parallel algorithms for the construction and 2:1 balance refinement of large linear octrees on distributed memory machines. Such octrees are used in many problems in computational science and engineering, e.g., object representation, image analysis, unstructured meshing, finite elements, adaptive mesh refinement, and N-body simulations. Fixed-size scalability and isogranular analysis of the algorithms using an MPI-based parallel implementation was performed on a variety of input data and demonstrated good scalability for different processor counts (1 to 1024 processors) on the Pittsburgh Supercomputing Center's TCS-1 AlphaServer. The results are consistent for different data distributions. Octrees with over a billion octants were constructed and balanced in less than a minute on 1024 processors. Like other existing algorithms for constructing and balancing octrees, our algorithms have O(N logN) work and O(N) storage complexity. Under reasonable assumptions on the distribution of octants and the work per octant, the parallel time complexity is O( N/np log( N/np ) + np log np), where N is the size of the final linear octree and np is the number of processors. [ABSTRACT FROM AUTHOR]

Large-scale optimization of systems governed by partial differential equations (PDEs) is a frontier problem in scientific computation. Reduced quasi-Newton sequential quadratic programming (SQP) methods are state-of-the-art approaches for such problems. These methods take full advantage of existing PDE solver technology and parallelize well. However, their algorithmic scalability is questionable; for certain problem classes they can be very slow to converge. In this two-part article we propose a new method for steady-state PDE-constrained optimization, based on the idea of using a full space Newton solver combined with an approximate reduced space quasi-Newton SQP preconditioner. The basic components of the method are Newton solution of the first-order optimality conditions that characterize stationarity of the Lagrangian function; Krylov solution of the Karush--Kuhn--Tucker (KKT) linear systems arising at each Newton iteration using a symmetric quasi-minimum residual method; preconditioning of the KKT system using an approximate state/decision variable decomposition that replaces the forward PDE Jacobians by their own preconditioners, and the decision space Schur complement (the reduced Hessian) by a BFGS approximation initialized by a two-step stationary method. Accordingly, we term the new method {\it Lagrange--Newton--Krylov--Schur} (LNKS). It is fully parallelizable, exploits the structure of available parallel algorithms for the PDE forward problem, and is locally quadratically convergent. In part I of this two-part article, we investigate the effectiveness of the KKT linear system solver. We test our method on two optimal control problems in which the state constraints are described by the steady-state Stokes equations. The objective is to minimize dissipation or the deviation from a given velocity field; the control variables are the boundary velocities. Numerical experiments on up to 256 Cray T3E processors and on an SGI Origin 2000 include scalability and performance assessment of the LNKS algorithm and comparisons with reduced SQP for up to $1,000,000$ state and 50,000 decision variables. In part II of the article, we address globalization and inexactness issues, and apply LNKS to the optimal control of the steady incompressible Navier--Stokes equations. [ABSTRACT FROM AUTHOR]

We present a new fast multipole method for particle simulations. The main feature of our algorithm is that it does not require the implementation of multipole expansions of the underlying kernel, and it is based only on kernel evaluations. Instead of using analytic expansions to represent the potential generated by sources inside a box of the hierarchical FMM tree, we use a continuous distribution of an equivalent density on a surface enclosing the box. To find this equivalent density, we match its potential to the potential of the original sources at a surface, in the far field, by solving local Dirichlet-type boundary value problems. The far-field evaluations are sparsified with singular value decomposition in 2D or fast Fourier transforms in 3D. We have tested the new method on the single and double layer operators for the Laplacian, the modified Laplacian, the Stokes, the modified Stokes, the Navier, and the modified Navier operators in two and three dimensions. Our numerical results indicate that our method compares very well with the best known implementations of the analytic FMM method for both the Laplacian and modified Laplacian kernels. Its advantage is the (relative) simplicity of the implementation and its immediate extension to more general kernels. [Copyright &y& Elsevier]

Abstract: Vesicles are locally-inextensible fluid membranes that can sustain bending. In this paper, we extend the study of Veerapaneni et al. [S.K. Veerapaneni, D. Gueyffier, G. Biros, D. Zorin, A numerical method for simulating the dynamics of 3D axisymmetric vesicles suspended in viscous flows, Journal of Computational Physics 228 (19) (2009) 7233–7249] to general non-axisymmetric vesicle flows in three dimensions. Although the main components of the algorithm are similar in spirit to the axisymmetric case (spectral approximation in space, semi-implicit time-stepping scheme), important new elements need to be introduced for a full 3D method. In particular, spatial quantities are discretized using spherical harmonics, and quadrature rules for singular surface integrals need to be adapted to this case; an algorithm for surface reparameterization is needed to ensure stability of the time-stepping scheme, and spectral filtering is introduced to maintain reasonable accuracy while minimizing computational costs. To characterize the stability of the scheme and to construct preconditioners for the iterative linear system solvers used in the semi-implicit time-stepping scheme, we perform a spectral analysis of the evolution operator on the unit sphere. By introducing these algorithmic components, we obtain a time-stepping scheme that circumvents the stability constraint on the time-step and achieves spectral accuracy in space. We present results to analyze the cost and convergence rates of the overall scheme. To illustrate the applicability of the new method, we consider a few vesicle-flow interaction problems: a single vesicle in relaxation, sedimentation, shear flows, and many-vesicle flows. [Copyright &y& Elsevier]