MRI Based Bayesian Personalization of a Tumor Growth Model

MRI-based digital twins to improve treatment response of breast cancer by optimizing neoadjuvant chemotherapy regimens

We developed a practical framework to construct digital twins for predicting and optimizing triple-negative breast cancer (TNBC) response to neoadjuvant chemotherapy (NAC). This study employed 105 TNBC patients from the ARTEMIS trial (NCT02276443, registered on 10/21/2014) who received Adriamycin/Cytoxan (A/C)-Taxol (T). Digital twins were established by calibrating a biology-based mathematical model to patient-specific MRI data, which accurately predicted pathological complete response (pCR) with an AUC of 0.82. We then used each patient's twin to theoretically optimize outcome by identifying their optimal A/C-T schedule from 128 options. The patient-specifically optimized treatment yielded a significant improvement in pCR rate of 20.95-24.76%. Retrospective validation was conducted by virtually treating the twins with AC-T schedules from historical trials and obtaining identical observations on outcomes: bi-weekly A/C-T outperforms tri-weekly A/C-T, and weekly/bi-weekly T outperforms tri-weekly T. This proof-of-principle study demonstrates that our digital twin framework provides a practical methodology to identify patient-specific TNBC treatment schedules.

Personalized predictions of Glioblastoma infiltration: Mathematical models, Physics-Informed Neural Networks and multimodal scans

Predicting the infiltration of Glioblastoma (GBM) from medical MRI scans is crucial for understanding tumor growth dynamics and designing personalized radiotherapy treatment plans. Mathematical models of GBM growth can complement the data in the prediction of spatial distributions of tumor cells. However, this requires estimating patient-specific parameters of the model from clinical data, which is a challenging inverse problem due to limited temporal data and the limited time between imaging and diagnosis. This work proposes a method that uses Physics-Informed Neural Networks (PINNs) to estimate patient-specific parameters of a reaction-diffusion partial differential equation (PDE) model of GBM growth from a single 3D structural MRI snapshot. PINNs embed both the data and the PDE into a loss function, thus integrating theory and data. Key innovations include the identification and estimation of characteristic non-dimensional parameters, a pre-training step that utilizes the non-dimensional parameters and a fine-tuning step to determine the patient specific parameters. Additionally, the diffuse-domain method is employed to handle the complex brain geometry within the PINN framework. The method is validated on both synthetic and patient datasets, showing promise for personalized GBM treatment through parametric inference within clinically relevant timeframes.

Mathematical modeling of multicellular tumor spheroids quantifies inter-patient and intra-tumor heterogeneity

In the study of brain tumors, patient-derived three-dimensional sphere cultures provide an important tool for studying emerging treatments. The growth of such spheroids depends on the combined effects of proliferation and migration of cells, but it is challenging to make accurate distinctions between increase in cell number versus the radial movement of cells. To address this, we formulate a novel model in the form of a system of two partial differential equations (PDEs) incorporating both migration and growth terms, and show that it more accurately fits our data compared to simpler PDE models. We show that traveling-wave speeds are strongly associated with population heterogeneity. Having fitted the model to our dataset we show that a subset of the cell lines are best described by a "Go-or-Grow"-type model, which constitutes a special case of our model. Finally, we investigate whether our fitted model parameters are correlated with patient age and survival.

Bayesian Inference of Tissue Heterogeneity for Individualized Prediction of Glioma Growth

Reliably predicting the future spread of brain tumors using imaging data and on a subject-specific basis requires quantifying uncertainties in data, biophysical models of tumor growth, and spatial heterogeneity of tumor and host tissue. This work introduces a Bayesian framework to calibrate the two-/three-dimensional spatial distribution of the parameters within a tumor growth model to quantitative magnetic resonance imaging (MRI) data and demonstrates its implementation in a pre-clinical model of glioma. The framework leverages an atlas-based brain segmentation of grey and white matter to establish subject-specific priors and tunable spatial dependencies of the model parameters in each region. Using this framework, the tumor-specific parameters are calibrated from quantitative MRI measurements early in the course of tumor development in four rats and used to predict the spatial development of the tumor at later times. The results suggest that the tumor model, calibrated by animal-specific imaging data at one time point, can accurately predict tumor shapes with a Dice coefficient 0.89. However, the reliability of the predicted volume and shape of tumors strongly relies on the number of earlier imaging time points used for calibrating the model. This study demonstrates, for the first time, the ability to determine the uncertainty in the inferred tissue heterogeneity and the model-predicted tumor shape.

Estimation of Parameter Distributions for Reaction-Diffusion Equations with Competition using Aggregate Spatiotemporal Data

Reaction-diffusion equations have been used to model a wide range of biological phenomenon related to population spread and proliferation from ecology to cancer. It is commonly assumed that individuals in a population have homogeneous diffusion and growth rates; however, this assumption can be inaccurate when the population is intrinsically divided into many distinct subpopulations that compete with each other. In previous work, the task of inferring the degree of phenotypic heterogeneity between subpopulations from total population density has been performed within a framework that combines parameter distribution estimation with reaction-diffusion models. Here, we extend this approach so that it is compatible with reaction-diffusion models that include competition between subpopulations. We use a reaction-diffusion model of glioblastoma multiforme, an aggressive type of brain cancer, to test our approach on simulated data that are similar to measurements that could be collected in practice. We use Prokhorov metric framework and convert the reaction-diffusion model to a random differential equation model to estimate joint distributions of diffusion and growth rates among heterogeneous subpopulations. We then compare the new random differential equation model performance against other partial differential equation models' performance. We find that the random differential equation is more capable at predicting the cell density compared to other models while being more time efficient. Finally, we use k-means clustering to predict the number of subpopulations based on the recovered distributions.

Learn-Morph-Infer: A new way of solving the inverse problem for brain tumor modeling

Current treatment planning of patients diagnosed with a brain tumor, such as glioma, could significantly benefit by accessing the spatial distribution of tumor cell concentration. Existing diagnostic modalities, e.g. magnetic resonance imaging (MRI), contrast sufficiently well areas of high cell density. In gliomas, however, they do not portray areas of low cell concentration, which can often serve as a source for the secondary appearance of the tumor after treatment. To estimate tumor cell densities beyond the visible boundaries of the lesion, numerical simulations of tumor growth could complement imaging information by providing estimates of full spatial distributions of tumor cells. Over recent years a corpus of literature on medical image-based tumor modeling was published. It includes different mathematical formalisms describing the forward tumor growth model. Alongside, various parametric inference schemes were developed to perform an efficient tumor model personalization, i.e. solving the inverse problem. However, the unifying drawback of all existing approaches is the time complexity of the model personalization which prohibits a potential integration of the modeling into clinical settings. In this work, we introduce a deep learning based methodology for inferring the patient-specific spatial distribution of brain tumors from T1Gd and FLAIR MRI medical scans. Coined as Learn-Morph-Infer, the method achieves real-time performance in the order of minutes on widely available hardware and the compute time is stable across tumor models of different complexity, such as reaction-diffusion and reaction-advection-diffusion models. We believe the proposed inverse solution approach not only bridges the way for clinical translation of brain tumor personalization but can also be adopted to other scientific and engineering domains.

Modeling of glioma growth using modified reaction-diffusion equation on brain MR images

BACKGROUND AND OBJECTIVE

Modeling of glioma growth and evolution is of key importance for cancer diagnosis, predicting clinical progression and improving treatment outcomes of neurosurgery. However, existing models are unable to characterize spatial variations of the proliferation and infiltration of tumor cells, making it difficult to achieve accurate prediction of tumor growth.

METHODS

In this paper, a new growth model of brain tumor using a reaction-diffusion equation on brain magnetic resonance images is proposed. Both the heterogeneity of brain tissue and the density of tumor cells are used to estimate the proliferation and diffusion coefficients of brain tumor cells. The diffusion coefficient that characterizes tumor diffusion and infiltration is calculated based on the ratio of tissues (white and gray matter), while the proliferation coefficient is evaluated using the spatial gradient of tumor cells. In addition, a parameter space is constructed using inverse distance weighted interpolation to describe the spatial distribution of proliferation coefficient.

RESULTS

The glioma growth predicted by the proposed model were tested by comparing with the real magnetic resonance images of the patients. Experiments and simulation results show that the proposed method achieves accurate modeling of glioma growth. The interpolation-based growth model has higher average dice score of 0.0647 and 0.0545, and higher average Jaccard index of 0.0673 and 0.0573, respectively, compared to the uniform- and gradient-based growth models.

CONCLUSIONS

The experimental results demonstrate the feasibility of calculating the proliferation and diffusion coefficients of the growth model based on patient-specific anatomy. The parameter space that characterizes spatial distribution of proliferation and diffusion coefficients is established and incorporated into the simulation of glioma growth. It enables to obtain patient-specific models about glioma growth by estimating and calibrating only a few model parameters.

Opportunities for improving brain cancer treatment outcomes through imaging-based mathematical modeling of the delivery of radiotherapy and immunotherapy

Immunotherapy has become a fourth pillar in the treatment of brain tumors and, when combined with radiation therapy, may improve patient outcomes and reduce the neurotoxicity. As with other combination therapies, the identification of a treatment schedule that maximizes the synergistic effect of radiation- and immune-therapy is a fundamental challenge. Mechanism-based mathematical modeling is one promising approach to systematically investigate therapeutic combinations to maximize positive outcomes within a rigorous framework. However, successful clinical translation of model-generated combinations of treatment requires patient-specific data to allow the models to be meaningfully initialized and parameterized. Quantitative imaging techniques have emerged as a promising source of high quality, spatially and temporally resolved data for the development and validation of mathematical models. In this review, we will present approaches to personalize mechanism-based modeling frameworks with patient data, and then discuss how these techniques could be leveraged to improve brain cancer outcomes through patient-specific modeling and optimization of treatment strategies.

Integrating mechanism-based modeling with biomedical imaging to build practical digital twins for clinical oncology

Digital twins employ mathematical and computational models to virtually represent a physical object (e.g., planes and human organs), predict the behavior of the object, and enable decision-making to optimize the future behavior of the object. While digital twins have been widely used in engineering for decades, their applications to oncology are only just emerging. Due to advances in experimental techniques quantitatively characterizing cancer, as well as advances in the mathematical and computational sciences, the notion of building and applying digital twins to understand tumor dynamics and personalize the care of cancer patients has been increasingly appreciated. In this review, we present the opportunities and challenges of applying digital twins in clinical oncology, with a particular focus on integrating medical imaging with mechanism-based, tissue-scale mathematical modeling. Specifically, we first introduce the general digital twin framework and then illustrate existing applications of image-guided digital twins in healthcare. Next, we detail both the imaging and modeling techniques that provide practical opportunities to build patient-specific digital twins for oncology. We then describe the current challenges and limitations in developing image-guided, mechanism-based digital twins for oncology along with potential solutions. We conclude by outlining five fundamental questions that can serve as a roadmap when designing and building a practical digital twin for oncology and attempt to provide answers for a specific application to brain cancer. We hope that this contribution provides motivation for the imaging science, oncology, and computational communities to develop practical digital twin technologies to improve the care of patients battling cancer.

Fast Posterior Estimation of Cardiac Electrophysiological Model Parameters via Bayesian Active Learning

Probabilistic estimation of cardiac electrophysiological model parameters serves an important step toward model personalization and uncertain quantification. The expensive computation associated with these model simulations, however, makes direct Markov Chain Monte Carlo (MCMC) sampling of the posterior probability density function (pdf) of model parameters computationally intensive. Approximated posterior pdfs resulting from replacing the simulation model with a computationally efficient surrogate, on the other hand, have seen limited accuracy. In this study, we present a Bayesian active learning method to directly approximate the posterior pdf function of cardiac model parameters, in which we intelligently select training points to query the simulation model in order to learn the posterior pdf using a small number of samples. We integrate a generative model into Bayesian active learning to allow approximating posterior pdf of high-dimensional model parameters at the resolution of the cardiac mesh. We further introduce new acquisition functions to focus the selection of training points on better approximating the shape rather than the modes of the posterior pdf of interest. We evaluated the presented method in estimating tissue excitability in a 3D cardiac electrophysiological model in a range of synthetic and real-data experiments. We demonstrated its improved accuracy in approximating the posterior pdf compared to Bayesian active learning using regular acquisition functions, and substantially reduced computational cost in comparison to existing standard or accelerated MCMC sampling.

Markov chain Monte Carlo with Gaussian processes for fast parameter estimation and uncertainty quantification in a 1D fluid-dynamics model of the pulmonary circulation

The past few decades have witnessed an explosive synergy between physics and the life sciences. In particular, physical modelling in medicine and physiology is a topical research area. The present work focuses on parameter inference and uncertainty quantification in a 1D fluid-dynamics model for quantitative physiology: the pulmonary blood circulation. The practical challenge is the estimation of the patient-specific biophysical model parameters, which cannot be measured directly. In principle this can be achieved based on a comparison between measured and predicted data. However, predicting data requires solving a system of partial differential equations (PDEs), which usually have no closed-form solution, and repeated numerical integrations as part of an adaptive estimation procedure are computationally expensive. In the present article, we demonstrate how fast parameter estimation combined with sound uncertainty quantification can be achieved by a combination of statistical emulation and Markov chain Monte Carlo (MCMC) sampling. We compare a range of state-of-the-art MCMC algorithms and emulation strategies, and assess their performance in terms of their accuracy and computational efficiency. The long-term goal is to develop a method for reliable disease prognostication in real time, and our work is an important step towards an automatic clinical decision support system.

Evaluating patient-specific neoadjuvant regimens for breast cancer via a mathematical model constrained by quantitative magnetic resonance imaging data

The ability to accurately predict response and then rigorously optimize a therapeutic regimen on a patient-specific basis, would transform oncology. Toward this end, we have developed an experimental-mathematical framework that integrates quantitative magnetic resonance imaging (MRI) data into a biophysical model to predict patient-specific treatment response of locally advanced breast cancer to neoadjuvant therapy. Diffusion-weighted and dynamic contrast-enhanced MRI data is collected prior to therapy, after 1 cycle of therapy, and at the completion of the first therapeutic regimen. The model is initialized and calibrated with the first 2 patient-specific MRI data sets to predict response at the third, which is then compared to patient outcomes (N = 18). The model's predictions for total cellularity, total volume, and the longest axis at the completion of the regimen are significant within expected measurement precision (P< 0.05) and strongly correlated with measured response (P < 0.01). Further, we use the model to investigate, in silico, a range of (practical) alternative treatment plans to achieve the greatest possible tumor control for each individual in a subgroup of patients (N = 13). The model identifies alternative dosing strategies predicted to achieve greater tumor control compared to the standard of care for 12 of 13 patients (P < 0.01). In summary, a predictive, mechanism-based mathematical model has demonstrated the ability to identify alternative treatment regimens that are forecasted to outperform the therapeutic regimens the patients clinically. This has important implications for clinical trial design with the opportunity to alter oncology care in the future.

Quantifying Uncertainty and Robustness in a Biomathematical Model-Based Patient-Specific Response Metric for Glioblastoma

Purpose

Glioblastomas, lethal primary brain tumors, are known for their heterogeneity and invasiveness. A growing body of literature has been developed demonstrating the clinical relevance of a biomathematical model, the proliferation-invasion model, of glioblastoma growth. Of interest here is the development of a treatment response metric, days gained (DG). This metric is based on individual tumor kinetics estimated through segmented volumes of hyperintense regions on T1-weighted gadolinium-enhanced and T2-weighted magnetic resonance images. This metric was shown to be prognostic of time to progression. Furthermore, it was shown to be more prognostic of outcome than standard response metrics. Although promising, the original article did not account for uncertainty in the calculation of the DG metric, leaving the robustness of this cutoff in question.

Methods

We harnessed the Bayesian framework to consider the impact of two sources of uncertainty: (1) image acquisition and (2) interobserver error in image segmentation. We first used synthetic data to characterize what nonerror variants are influencing the final uncertainty in the DG metric. We then considered the original patient cohort to investigate clinical patterns of uncertainty and to determine how robust this metric is for predicting time to progression and overall survival.

Results

Our results indicate that the key clinical variants are the time between pretreatment images and the underlying tumor growth kinetics, matching our observations in the clinical cohort. Finally, we demonstrated that for this cohort, there was a continuous range of cutoffs between 94 and 105 for which the prediction of the time to progression was over 80% reliable.

Conclusion

Although additional validation must be performed, this work represents a key step in ascertaining the clinical utility of this metric.

Embedding high-dimensional Bayesian optimization via generative modeling: Parameter personalization of cardiac electrophysiological models

The estimation of patient-specific tissue properties in the form of model parameters is important for personalized physiological models. Because tissue properties are spatially varying across the underlying geometrical model, it presents a significant challenge of high-dimensional (HD) optimization at the presence of limited measurement data. A common solution to reduce the dimension of the parameter space is to explicitly partition the geometrical mesh. In this paper, we present a novel concept that uses a generative variational auto-encoder (VAE) to embed HD Bayesian optimization into a low-dimensional (LD) latent space that represents the generative code of HD parameters. We further utilize VAE-encoded knowledge about the generative code to guide the exploration of the search space. The presented method is applied to estimating tissue excitability in a cardiac electrophysiological model in a range of synthetic and real-data experiments, through which we demonstrate its improved accuracy and substantially reduced computational cost in comparison to existing methods that rely on geometry-based reduction of the HD parameter space.

Longitudinal brain tumor segmentation prediction in MRI using feature and label fusion

This work proposes a novel framework for brain tumor segmentation prediction in longitudinal multi-modal MRI scans, comprising two methods; feature fusion and joint label fusion (JLF). The first method fuses stochastic multi-resolution texture features with tumor cell density feature to obtain tumor segmentation predictions in follow-up timepoints using data from baseline pre-operative timepoint. The cell density feature is obtained by solving the 3D reaction-diffusion equation for biophysical tumor growth modelling using the Lattice-Boltzmann method. The second method utilizes JLF to combine segmentation labels obtained from (i) the stochastic texture feature-based and Random Forest (RF)-based tumor segmentation method; and (ii) another state-of-the-art tumor growth and segmentation method, known as boosted Glioma Image Segmentation and Registration (GLISTRboost, or GB). We quantitatively evaluate both proposed methods using the Dice Similarity Coefficient (DSC) in longitudinal scans of 9 patients from the public BraTS 2015 multi-institutional dataset. The evaluation results for the feature-based fusion method show improved tumor segmentation prediction for the whole tumor(DSC _WT = 0.314, p = 0.1502), tumor core (DSC _TC = 0.332, p = 0.0002), and enhancing tumor (DSC _ET = 0.448, p = 0.0002) regions. The feature-based fusion shows some improvement on tumor prediction of longitudinal brain tumor tracking, whereas the JLF offers statistically significant improvement on the actual segmentation of WT and ET (DSC _WT = 0.85 ± 0.055, DSC _ET = 0.837 ± 0.074), and also improves the results of GB. The novelty of this work is two-fold: (a) exploit tumor cell density as a feature to predict brain tumor segmentation, using a stochastic multi-resolution RF-based method, and (b) improve the performance of another successful tumor segmentation method, GB, by fusing with the RF-based segmentation labels.

Coupling brain-tumor biophysical models and diffeomorphic image registration

We present SIBIA (Scalable Integrated Biophysics-based Image Analysis), a framework for joint image registration and biophysical inversion and we apply it to analyze MR images of glioblastomas (primary brain tumors). We have two applications in mind. The first one is normal-to-abnormal image registration in the presence of tumor-induced topology differences. The second one is biophysical inversion based on single-time patient data. The underlying optimization problem is highly non-linear and non-convex and has not been solved before with a gradient-based approach. Given the segmentation of a normal brain MRI and the segmentation of a cancer patient MRI, we determine tumor growth parameters and a registration map so that if we "grow a tumor" (using our tumor model) in the normal brain and then register it to the patient image, then the registration mismatch is as small as possible. This "coupled problem" two-way couples the biophysical inversion and the registration problem. In the image registration step we solve a large-deformation diffeomorphic registration problem parameterized by an Eulerian velocity field. In the biophysical inversion step we estimate parameters in a reaction-diffusion tumor growth model that is formulated as a partial differential equation (PDE). In SIBIA, we couple these two sub-components in an iterative manner. We first presented the components of SIBIA in "Gholami et al., Framework for Scalable Biophysics-based Image Analysis, IEEE/ACM Proceedings of the SC2017", in which we derived parallel distributed memory algorithms and software modules for the decoupled registration and biophysical inverse problems. In this paper, our contributions are the introduction of a PDE-constrained optimization formulation of the coupled problem, and the derivation of a Picard iterative solution scheme. We perform extensive tests to experimentally assess the performance of our method on synthetic and clinical datasets. We demonstrate the convergence of the SIBIA optimization solver in different usage scenarios. We demonstrate that using SIBIA, we can accurately solve the coupled problem in three dimensions (256³ resolution) in a few minutes using 11 dual-x86 nodes.

Subject-specific Brain Tumor Growth Modelling via An Efficient Bayesian Inference Framework

An accurate prediction of brain tumor progression is crucial for optimized treatment of the tumors. Gliomas are primarily treated by combining surgery, external beam radiotherapy, and chemotherapy. Among them, radiotherapy is a non-invasive and effective therapy, and an understanding of tumor growth will allow better therapy planning. In particular, estimating parameters associated with tumor growth, such as the diffusion coefficient and proliferation rate, is crucial to accurately characterize physiology of tumor growth and to develop predictive models of tumor infiltration and recurrence. Accurate parameter estimation, however, is a challenging task due to inaccurate tumor boundaries and the approximation of the tumor growth model. Here, we introduce a Bayesian framework for a subject-specific tumor growth model that estimates the tumor parameters effectively. This is achieved by using an improved elliptical slice sampling method based on an adaptive sample region. Experimental results on clinical data demonstrate that the proposed method provides a higher acceptance rate, while preserving the parameter estimation accuracy, compared with other state-of-the-art methods such as Metropolis-Hastings and elliptical slice sampling without any modification. Our approach has the potential to provide a method to individualize therapy, thereby offering an optimized treatment.

Quantifying the uncertainty in model parameters using Gaussian process-based Markov chain Monte Carlo in cardiac electrophysiology

Model personalization requires the estimation of patient-specific tissue properties in the form of model parameters from indirect and sparse measurement data. Moreover, a low-dimensional representation of the parameter space is needed, which often has a limited ability to reveal the underlying tissue heterogeneity. As a result, significant uncertainty can be associated with the estimated values of the model parameters which, if left unquantified, will lead to unknown variability in model outputs that will hinder their reliable clinical adoption. Probabilistic estimation of model parameters, however, remains an unresolved challenge. Direct Markov Chain Monte Carlo (MCMC) sampling of the posterior distribution function (pdf) of the parameters is infeasible because it involves repeated evaluations of the computationally expensive simulation model. To accelerate this inference, one popular approach is to construct a computationally efficient surrogate and sample from this approximation. However, by sampling from an approximation, efficiency is gained at the expense of sampling accuracy. In this paper, we address this issue by integrating surrogate modeling of the posterior pdf into accelerating the Metropolis-Hastings (MH) sampling of the exact posterior pdf. It is achieved by two main components: (1) construction of a Gaussian process (GP) surrogate of the exact posterior pdf by actively selecting training points that allow for a good global approximation accuracy with a focus on the regions of high posterior probability; and (2) use of the GP surrogate to improve the proposal distribution in MH sampling, in order to improve the acceptance rate. The presented framework is evaluated in its estimation of the local tissue excitability of a cardiac electrophysiological model in both synthetic data experiments and real data experiments. In addition, the obtained posterior distributions of model parameters are interpreted in relation to the factors contributing to parameter uncertainty, including different low-dimensional representations of the parameter space, parameter non-identifiability, and parameter correlations.

Low grade glioma growth modeling considering chemotherapy and radiotherapy effects from magnetic resonance images

Studying tumor growth using mathematical models from magnetic resonance (MR) images is an important application that is believed to play a major role in cancer treatment by predicting tumor evolution, quantifying the response to therapy, and treatment planning. Reaction diffusion is the most popular model because of its simplicity and consistency with the biological growth process. However, most of the current growth models focus on presurgical images and likely without treatment. In this paper, we propose a new reaction diffusion model to consider the chemotherapy and radiotherapy effects on the tumor growth modelling for patients with low grade glioma. The proposed model does not consider the tensor information from diffusion tensor imaging. Instead it uses a weighted parameter to promote higher diffusivity in white matter. The radiotherapy and chemotherapy effects are considered as a loss terms in the proposed model. The preliminary results of the proposed model on synthetic and 2 real MR images show that, our model can effectively simulate tumor growth with high accuracies when treatments are administrated to low grade glioma patients.