Where did the tumor start? An inverse solver with sparse localization for tumor growth models

TL;DR

This paper introduces a novel inverse solver for tumor growth models that uses sparse localization and regularization to accurately identify the tumor's origin, infiltration, and proliferation from MRI data, outperforming existing methods.

Contribution

The paper presents a new inverse solver with sparsity-based regularization for tumor localization, improving initial condition recovery in reaction-diffusion models from MRI scans.

Findings

Successfully localizes tumor origin in synthetic and clinical images.

Outperforms existing solvers in initial condition localization.

Provides qualitatively different solutions on clinical data.

Abstract

We present a numerical scheme for solving an inverse problem for parameter estimation in tumor growth models for glioblastomas, a form of aggressive primary brain tumor. The growth model is a reaction-diffusion partial differential equation (PDE) for the tumor concentration. We use a PDE-constrained optimization formulation for the inverse problem. The unknown parameters are the reaction coefficient (proliferation), the diffusion coefficient (infiltration), and the initial condition field for the tumor PDE. Segmentation of Magnetic Resonance Imaging (MRI) scans from a single time snapshot drive the inverse problem where segmented tumor regions serve as partial observations of the tumor concentration. The precise time relative to tumor initiation is unknown, which poses an additional difficulty for inversion. We perform a frozen-coefficient spectral analysis and show that the inverse problem is severely ill-posed. We introduce a biophysically motivated regularization on the tumor initial condition. In particular, we assume that the tumor starts at a few locations (enforced with a sparsity constraint) and that the initial condition magnitude in the maximum norm equals one. We solve the resulting optimization problem using an inexact quasi-Newton method combined with a compressive sampling algorithm for the sparsity constraint. Our implementation uses PETSc and AccFFT libraries. We conduct numerical experiments on synthetic and clinical images to highlight the improved performance of our solver over an existing solver that uses a two-norm regularization for the calibration parameters. The existing solver is unable to localize the initial condition. Our new solver can localize the initial condition and recover infiltration and proliferation. In clinical datasets (for which the ground truth is unknown), our solver results in qualitatively different solutions compared to the existing solver.

Tables

Table 1. Table 1: Notation for the inverse tumor problem formulation.

Description	Notation
spatial coordinates	$𝒙=(x,y,z)$
time	$t$
healthy tissue probability maps (see Eq. (1))	$𝒎\left(𝒙\right)$
patient tumor from imaging (see Eq. (9a))	$d\left(𝒙\right)$
initial tumor concentration (see Eq. (9c))	$c(𝒙,0)$
simulated tumor concentration at time $t$ (see Eq. (9b))	$c(𝒙,t)$
reaction operator (see Eq. (7))	$ℛ$
scalar reaction coefficient (see Eq. (8))	$\rho$
diffusion operator (see Eq. (5))	$𝒟$
scalar diffusion coefficient (see Eq. (6))	$\kappa$
parameterization basis functions (see Eq. (10))	$\varphi \left(𝒙\right)={\left\{{\varphi }_i\left(𝒙\right)\right\}}_{i=1}^m$
parameterization coefficients in ${𝐑}^m$ (see Eq. (10))	$𝒑$
observation operator (see Eq. (11))	$𝑶\left(𝒙\right)$
tumor concentration detection threshold (see Eq. (11))	$c_d$
regularization parameter (see Eq. (9a))	$\beta$

Table 2. Table 2: Numerical parameters and their values used for tumor inversion.

Parameter	Value
spatial discretization for synthetic tumors	${128}^3$
volume fraction for Gaussian basis selection, ${\tau }_V$	$0.99$
standard deviation of Gaussian basis functions, $\sigma$	$2\pi /64$
relative gradient tolerance, ${\tau }_g$	$1\text{E}-5$
relative objective tolerance, ${\tau }_J$	$1\text{E}-5$
maximum quasi-Newton iterations	$50$
maximum compressive sampling iterations	$2$
sparsity level, $s$	10
regularization parameter, $\beta$	$1\text{E}-4$
initial guess for reaction coefficient, ${\rho }^0$	$15$

Equations

$${\partial }_{t}c(\mathbf{x},t)=\kappa \mathcal{D}c(\mathbf{x},t)+\rho \mathcal{R}c(\mathbf{x},t),t\in [0,T],\text{mbox}withinitialconditionc(\mathbf{x},0),$$

$$\boldsymbol{m}(\boldsymbol{x})=\left(\begin{array}[]{c}m_{g}(\boldsymbol{x})\\ m_{w}(\boldsymbol{x})\\ m_{f}(\boldsymbol{x})\end{array}\right)\in\mathbf{R}^{3},$$

$${\partial }_{\tilde{t}}c-T\kappa \mathcal{D}c-T\rho \mathcal{R}c=0.$$

$${\partial }_{t}c-\kappa \mathcal{D}c-\rho \mathcal{R}c$$

$$c(0)$$

$$\mathcal{D}c=\text{mbox}div{k}_{\mathbf{m}}\nabla c,$$

$$k_{\mathbf{m}}(\mathbf{x})={\kappa }_g{m}_g(\mathbf{x})+{\kappa }_w{m}_w(\mathbf{x}),$$

$$\mathcal{R}c={\rho }_{\mathbf{m}}c(1-c),$$

$${\rho }_{\mathbf{m}}(\mathbf{x})={\rho }_g{m}_g(\mathbf{x})+{\rho }_w{m}_w(\mathbf{x}),$$

$$\underset{\mathbf{p},\kappa ,\rho }{\min }\mathcal{J}(\mathbf{p},\kappa ,\rho ):=\frac{1}{2}\Vert \mathbf{O}c(1)-d{\Vert }_{L_2(\Omega )}^2+\frac{\beta }{2}\Vert {\mathbf{\varphi }}^T\mathbf{p}{\Vert }_{L_2(\Omega )}^2$$

$${\partial }_{t}c-\kappa \mathcal{D}c-\rho \mathcal{R}c$$

$$c(0)$$

$$c(\mathbf{x},0)={\mathbf{\varphi }}^T(\mathbf{x})\mathbf{p}=\sum _{i=1}^m{\varphi }_i(\mathbf{x})p_i,\mathbf{p}\in {\mathbf{R}}^m.$$

$$\mathbf{O}(\mathbf{x})=\{\begin{array}{ll}1,& \text{if }d(\mathbf{x})\ge c_d\\ 0,& \text{otherwise.}\end{array}$$

$${\partial }_{t}c-\rho c-\kappa {\partial }_x^{2}c$$

$${{\partial }_{x}c|}_{x=0,\pi }$$

$$c(x,0)-g(x)$$

$$c(x,t)=\sum _{n=0}^{\infty }{g}_n{e}^{(\rho -n^2\kappa )t}\cos (nx),$$

$$\underset{g_n,\rho ,\kappa }{\min }\mathcal{J}:=\frac{1}{2}\sum _{n=0}^N{\left(g_n{e}^{(\rho -n^2\kappa )}-d_n\right)}^2.$$

$${\partial }_{g_n}\mathcal{J}$$

$${\partial }_{\rho }\mathcal{J}$$

$${\partial }_{\kappa }\mathcal{J}$$

$$g_n{e}^{(\rho -n^2\kappa )}-d_n=0.$$

$$g(x)=\sum _{i=0}^m{\varphi }_i(x)p_i={\mathbf{\varphi }}^T(x)\mathbf{p},$$

$$c(x,t)=\sum _{n=0}^N{\mathbf{\varphi }}_n^T\mathbf{p}{e}^{(\rho -n^2\kappa )}\cos (nx),$$

$$\underset{\mathbf{p},\rho ,\kappa }{\min }\mathcal{J}:=\frac{1}{2}\sum _{n=0}^N{\left({\mathbf{\varphi }}_n^T\mathbf{p}{e}^{(\rho -n^2\kappa )}-d_n\right)}^2.$$

$$\mathbf{p}$$

$$\kappa$$

$$\rho$$

$$\max ({\mathbf{\varphi }}^T\mathbf{p})$$

$$\text{mbox}forward{\partial }_{t}c-\kappa \mathcal{D}c-\rho \mathcal{R}c$$

$$c(0)-{\mathbf{\varphi }}^T\mathbf{p}$$

$$\text{mbox}adjoint-{\partial }_t\alpha -\kappa \mathcal{D}\alpha -\rho {\partial }_{\mathbf{c}}(\mathcal{R}c)\alpha$$

$${\mathbf{O}}^T(\mathbf{O}c(1)-d)+\alpha (1)$$

$$\text{mbox}gradient{\mathbf{g}}_{\mathbf{p}}=\beta \left({\int }_{\Omega }\mathbf{\varphi }{\mathbf{\varphi }}^{T}d\mathbf{x}\right)\mathbf{p}-{\int }_{\Omega }\mathbf{\varphi }\alpha (0)d\mathbf{x}$$

$$g_{\kappa }={\int }_0^1{\int }_{\Omega }{k}_{\mathbf{m}}\left((\nabla c{)}^T\nabla \alpha \right)d\mathbf{x}dt$$

$$g_{\rho }={\int }_0^1{\int }_{\Omega }{\rho }_{\mathbf{m}}\left(-\alpha c(1-c)\right)d\mathbf{x}dt$$

$${\varphi }_i(\mathbf{x})=\exp \left(-\frac{1}{2{\sigma }^2}\Vert \mathbf{x}-{\mathbf{x}}_i{\Vert }_2^2\right)1(m_{CSF}(\mathbf{x})=0),$$

Full text

(Preprint submitted to Journal of Inverse Problems)

1 Introduction

The integration of biophysical macroscopic brain tumor growth models with Magnetic Resonance Imaging (MRI) scans is an important tool in clinical cancer research. Phenomenological growth models capture phenomena like tumor cell invasion and proliferation and can be used to compute patient-specific biomarkers. Such biomarkers can potentially facilitate diagnosis (e.g., tumor grading), prognosis (e.g., predicting recurrence and survival), and treatment (e.g., pre-operative planning and radiotherapy). Additionally, biophysical models can advance our understanding of the disease by using imaging data to test model-driven hypotheses on disease progression and treatment. For example, one such hypothesis is that the location of tumor initiation relates to distinct tumor mutations, thereby paving the way to molecularly-informed models of tumor growth [50]. The main problem we focus here is whether we can localize the initiation of the tumor while simultaneously estimating physics-based parameters that characterize infiltration (spreading of the tumor) and proliferation (growth).

We use a single-species, reaction-diffusion model of tumor growth (described in detail in §2.1.1),

[TABLE]

where $c(\boldsymbol{x},t)$ is the tumor volume fraction (or concentration, assuming constant density), $\kappa$ is a scalar diffusion parameter, $\rho$ is a scalar reaction (or growth) parameter, ${\mathcal{D}}$ is a diffusion operator (with free-flux boundary conditions), and ${\mathcal{R}}$ is an non-linear logistic growth operator. This is one of the most common tumor growth models owing to its simplicity (it has just two parameters and the initial condition) and good performance in capturing the basic characteristics of brain tumor growth observed in MRI scans [58, 53]. The parameters of this model are the reaction coefficient $\rho$ (modeling proliferation), the diffusion coefficient $\kappa$ (modeling invasion), and the initial condition $c(\boldsymbol{x},0)$. For the inverse problem (described in §2.1.2), we assume that we are given partial observations of $c(\boldsymbol{x},T)$ at some, unknown relative to the tumor initiation, time $T$ and we seek to reconstruct $\kappa$, $\rho$, and $c(\boldsymbol{x},0)$, which are unknown. The time horizon $T$ is also unknown but a simple calculation (see §2.1) shows that among $\kappa$, $\rho$, and $T$, only two are independent. Hence, we set $T=1$ and invert for the remaining parameters. Notice that $c(\boldsymbol{x},0)$ is a field. In our formulation, we parameterize it as a finite sum of radial functions. This re-parameterization is not only for regularization but mainly to allow for $\ell_{1}$ sparsity constraints and can be viewed as a modeling assumption.

The data used for the inverse problem are partial observations of the tumor concentration at a single time point. This scenario corresponds to the clinical setting in which we seek to estimate the tumor grade and aggressiveness from a single-time-point scan (as opposed to the case in which we have a series of scans taken at different points in time). Indeed, when a patient comes with symptoms, the tumor is typically quite large and treatment (surgery, chemotherapy, or radiation) starts immediately. Hence, it is not practical to assume multiple time snapshots of the tumor growth because otherwise we would have to model treatment (which comes with its own modeling challenges and introduces additional unknown parameters). Given this setup, the inverse problem is formulated as a PDE-constrained optimization problem. We seek to estimate the model parameters that minimize the mismatch between the model-generated tumor and the partially observed tumor concentration. This problem is highly ill-posed. Notice that just inverting for the initial condition problem is exponentially ill-posed (see §2.2). The inherent non-linearity in the tumor model and the strong ill-posedness of the associated inverse problem make this effort a computationally challenging task. In this paper, we discuss this ill-posedeness, and we propose a set of algorithms for its numerical solution.

1.1 Contributions

The main innovation in the paper is the formulation of an $\ell_{1}$-constrained inverse problem for the initial condition, and two model parameters, the reaction $\rho$ and the diffusion $\kappa$. Our contributions are summarized below.

• •

We formulate the inverse problem as a PDE-constrained optimization problem (see §2.1). We perform an analysis of the linearized operator to demonstrate the ill-posedness and the need for a specific regularization (see §2.2). Although an $\ell_{2}$ regularization (for all parameters) removes the ill-posedness, it does not result in good localization of the initial condition and tumor parameters. For this reason, we introduce an $\ell_{1}$ sparsity constraint to enable sparse111We denote an initial condition parameterization with only few non-zero coefficients as sparse. reconstruction for the initial condition, combined with a maximum norm condition for $c(\boldsymbol{x},0)$. This allows us to reconstruct the model parameters $\kappa$ and $\rho$. This regularization is biophysically motivated.

• •

We introduce a novel numerical scheme to simultaneously invert for all model parameters while preserving the sparse content of the initial tumor condition. We employ a quasi-Newton optimization solver along with a reduced-space parameterization of the tumor initial condition (see §3). The scheme alternates between an $\ell_{2}$ regularized solve and an $\ell_{1}$ projection. The underlying linear and non-linear solves are based on preconditioned inexact Krylov solves accelerated with PETSc and AccFFT libraries that support distributed-memory parallelism and enable fast solves for the challenging derivative-based 3D problem.

• •

We test our scheme on 3D clinical and synthetic datasets. The synthetic datasets resemble the complexity in tumor patterns observed in MRI scans and are used to verify our algorithm since we know the ground truth. With the clinical images, we examine the prediction error and compare with a classical $\ell_{2}$ regularization scheme. We report results that highlight the performance improvement of our solver over existing inversion strategies in estimating $c(\boldsymbol{x},0)$, $\kappa$, and $\rho$. We further demonstrate the robustness of our solver by testing it on noisy data and partial observations. This inversion setup is related to the clinical setting because the extent of tumor invasion is not fully observable in MRI scans.

1.2 Limitations

While the reaction-diffusion model can predict the whole tumor region, it does not distinguish between the different tumor phenotypes observed in imaging data (such as enhancing and necrotic tumor regions). Further, it does not account for any mechanical deformations on the ventricles and surrounding healthy tissue due to tumor growth. We have developed more complex models that do so [57], but their integration with the inverse problem is beyond the scope of this paper. We do not have direct observations of the tumor concentration, we only have the segmentation of the MRI scan into healthy tissue, edema, enhancing tumor, non-enhancing tumor, and necrotic tumor. In this paper, we ignore these glioma-specific details, and assign a volume fraction one to regions of enhancing, non-enhancing, or necrotic tumor. Our formulation requires the healthy patient brain scan, which is typically not available. The standard approach in medical imaging is to resort to methods such as image registration to construct a healthy atlas that resembles the observed healthy structures of the patient scans. This can be inaccurate for very large tumors. We are currently working on incorporating the ideas presented here with our previous work [56], in which we coupled diffeomorphic image registration with tumor inversion for this purpose.

1.3 Related work

Although there has been a lot of work on forward problems for tumor growth, there has been less work on inverse problems. The latter has different aspects. The first is the underlying biophysical model. The second is the inverse problem setup, observation operators and the existence of scans at multiple points, the noise models, inversion parameters and constraints. And the third one is the solution algorithm.

Regarding the underlying model, like us, most researchers focus on parameter calibration of a handful of model parameters using single-species reaction-diffusion equations [7, 22, 24, 28, 33, 39, 58, 59, 51]. While more complex models describing processes like mass effect, angiogenesis and chemotaxis [54, 48, 60, 23, 26] exist, they have not been considered for calibration due to theoretical and computational challenges. However, several groups, including ours, are working to address these challenges.

Regarding the inverse problem setup, in most studies the authors use scans from two or more time points and assume that the tumor concentration is fully observed. In that case the inverse problem is much more tractable since we need to invert for just $\kappa$ and $\rho$. Some authors use manual calibration, which is possible for a handful of parameters [7, 59]. But most research groups resort to derivative-free optimization algorithms including Bayesian methods [6, 32, 37, 38, 42, 64, 21, 35, 34]. In [33], the authors use the traveling wave approximation to the reaction-diffusion model to estimate the diffusivity by fixing the reaction coefficient and the source seed. They explain that the non-uniqueness of the solution renders the inversion of all parameters infeasible for their model. Derivative-based optimization methods have been employed in [9, 25, 15, 31, 37, 39, 49, 56].

Now we discuss, in some detail, attempts to reconstruct the initial condition of the tumor using a single scan. In our previous works [15, 56], we employed an $L_{2}$ regularization scheme for $c(\boldsymbol{x},0)$ and we inverted for $\kappa$ assuming known $\rho$. We used reduced-space parallel algorithms allowing for the calibration of highly complex tumor shapes in reasonable time. There have been a few studies exploring the identifiability of all parameters and their correlation with the tumor source. The traveling wave approximation (using an Eikonal equation) from [33], is extended in [52] to estimate a localized initial condition and the diffusion coefficient. To our knowledge, this is the only attempt to localize the initial condition of the tumor. Its advantage is that it is simple to implement and quite robust since it circumvents non-convexity. But it is limited to the specific Eikonal equation for the propagation of the tumor front, which is then time-reversed to localize the tumor. Unlike our scheme, this inversion scheme cannot be used with more realistic biophysical models. In [27], the authors use the 3D non-linear model and attempt to identify $c(\boldsymbol{x},0)$ by fixing the growth rate and the diffusion coefficient to set values. The scheme is essentially an unregularized gradient descent and assumes noise-free, full observations of $c(\boldsymbol{x},1)$. The authors consider only synthetic experiments and assume $\kappa$ and $\rho$ are the ground truth. In this idealistic setup everything works well and one can localize $c(\boldsymbol{x},0)$. As we will see, in the presence of noise, partial observations, and uncertainty in $\kappa$ and $\rho$, such a formulation fails.

Let us also mention that, in general (unrelated to tumor modeling), there is little work on sparsity constraints for initial condition inverse problems for reaction-diffusion PDEs. While some groups have used sparsity promoting approaches to estimate sources in the heat equation in a 1-D and 2-D setting [36, 13], there is not much work on sparse localization for the non-linear reaction-diffusion equations. In conclusion, there is no prior work on a generalizable scheme for inverting for a sparse initial condition for the tumor in addition to model parameters.

Regarding numerical methods for $\ell_{1}$ regularized inverse problems, there is a vast amount of work in literature. Most algorithms are based on compressive sampling pursuit [62, 11, 46, 45, 3], $\ell_{1}$ convex optimization approaches using methods like interior points [4, 30], projected gradient algorithms [12], or iterative shrinkage algorithms [29, 2, 18]. We base our algorithm on the enhanced greedy pursuit strategy described in [45], which yields optimal theoretical guarantees on convergence rates for quadratic mismatch functions. However, our overall problem is ill-posed and non-convex. To our knowledge, fast algorithms with theoretical guarantees for problems like ours do not exist.

1.4 Outline

In §2, we present the forward problem and the mathematical formulation of the parameter estimation inverse problem. We explore the invertibility of the unknown parameters of the growth model in a simplified 1D setting and introduce new assumptions to solve the parameter estimation problem. In §3, we provide details regarding the algorithm and numerical schemes used. Finally, we conduct numerical experiments to assess the performance of our solver in §4.

2 Methods

We introduce some notation and subsequently present the forward problem and the mathematical formulation of the inverse problem below.

2.1 Formulation

We have briefly discussed the forward model for the tumor concentration $c(\boldsymbol{x},t)$. It involves two operators, the diffusion operator ${\mathcal{D}}$ and the reaction operator ${\mathcal{R}}$. In order to define these operators, we need to take into account the material heterogeneity of the brain tissue. For example, tumors do not grow inside the ventricles, so the reaction and diffusion coefficients should be set to zero. Also, we need to account for the fact that we do not have the precise tissue type, but only probability maps related to the MRI scans. These probability maps are obtained by statistical analysis of manual or machine-learning based segmentations [19]. So, in our setting, a brain is described by separate probability maps for each brain tissue type, namely gray matter, white matter, and cerebrospinal fluid (CSF; which includes ventricles). We represent such probability maps as a vector field:

[TABLE]

where $\boldsymbol{x}=(x,y,z)$ are the spatial coordinates and $m_{g}(\boldsymbol{x})$, $m_{w}(\boldsymbol{x})$ and $m_{f}(\boldsymbol{x})$ denote the gray matter, white matter and CSF probability maps, respectively. These probability maps are obtained from binary segmentation labels for each tissue type of a healthy brain MRI scan.222In a pre-processing step, the binary labels are smoothed and re-scaled to probabilities. Our data $d(\boldsymbol{x})$ is derived from actual patient MRI tumor segmentations.333We assign $d(\boldsymbol{x})=1$ to the voxels that have been identified as tumor. More complex scenarios are possible, for example see [5], which uses apparent diffusion coefficient (ADC) values to identify tumor tissue. Figure 1 shows an example of a healthy brain MRI scan (and its associated segmentation) and a patient brain MRI scan (and its associated tumor segmentation). As we discussed the healthy patient brain scan is rarely available in practice. There have been many approaches to address this, such as using a statistical average of many individual brains or image registration methods to approximate the healthy patient brain (see [56] for details). For the most part of this paper, we assume that the healthy patient brain segmentation is given. For our test on an actual clinical scan, we use an atlas-based segmentation in §4.7.

2.1.1 Forward problem

The forward problem is a single-species, reaction-diffusion partial differential Eq. [22, 24, 28, 33, 39, 58, 44]. In this model we group all the MRI-observed tumor phenotypes into a single one. (For gliomas, these phenotypes are typically enhancing, non-enhancing, and necrotic. We do not include edema in the tumor although it is quite typical to have significant tumor infiltration in the edema region.) As mentioned in the introduction, the time horizon (the time since tumor initiation at which we observe the tumor concentration) is unknown. However, we can non-dimensionalize the PDE and use a normalized time horizon equal to one for tumor growth so that $t\in[0,1]$. We discuss this further below. Let $\tilde{t}=t/T$ be the normalized time. Then, the PDE can be written as:

[TABLE]

Denoting $\tilde{\kappa}=\kappa T$ and $\tilde{\rho}=\rho T$, we observe that the normalization of time is equivalent to a re-scaling of the reaction and diffusion coefficients. In the remaining of the paper, we use $t,\kappa,\rho$ instead of $\tilde{t},\tilde{\kappa},\tilde{\rho}$ to represent the re-scaled parameters. This normalization is because we cannot possibly invert for the three parameters $\rho$, $\kappa$, and $T$, since only two are independent. So, we set $T=1$ and we consider the reconstruction of $\rho$ and $\kappa$ (in addition to the initial condition). The brain scan is re-scaled to be in $\Omega=[0,2\pi]^{3}$. The non-dimensionalized forward problem is given as:

[TABLE]

where $c=c(\boldsymbol{x},t)$ is the tumor volume fraction (or mass concentration assuming constant density), $c(0)=c(\boldsymbol{x},0)$ is the initial tumor concentration set to $g=g(\boldsymbol{x})$, $\mathcal{D}$ is the diffusion operator, $\kappa$ is a scalar diffusion coefficient, $\rho$ is a scalar growth coefficient, and $\mathcal{R}$ is the reaction operator. We use periodic boundary conditions on $\partial\Omega$. In reality, we have no-flux boundary conditions on the boundary of the skull and CSF boundaries. These no-flux conditions are approximated using a penalty method, the details are discussed in [15]. The approximation amounts to having very small diffusivity in the ventricles and the image space surrounding the brain. We define $\mathcal{D}$ to be an inhomogeneous isotropic diffusion operator given by

[TABLE]

where $k_{\boldsymbol{m}}(\boldsymbol{x})$ is defined as:

[TABLE]

where $\kappa_{g}$ and $\kappa_{w}$ are scalar parameters which enforce inhomogeneity of the diffusion rate in the gray matter and white matter, respectively. Similarly, we define the reaction operator $\mathcal{R}$ as:

[TABLE]

with $\rho_{\boldsymbol{m}}(\boldsymbol{x})$ defined as:

[TABLE]

where $\rho_{g}$ and $\rho_{w}$ are scalar parameters which control the spatial variability of the growth rate in the gray matter and white matter, respectively. The tumor does not grow or diffuse in the cerebrospinal fluid-filled ventricles (CSF).

2.1.2 Inverse problem

Our formulation of the inverse problem requires the definition of the observation operator, the mismatch function, and the initial condition parameterization. Existing MRI technologies (the predominant imaging modality used in neuro-oncology) do not provide the ability to directly observe the extent of tumor concentration everywhere in the parenchyma. For this reason, we introduce an observation operator $\boldsymbol{O}(\boldsymbol{x})$ that is determined by the clearly observable tumor margin. We also assume that in that region the tumor concentration is equal to one. We now present the formulation for the inverse problem and explain its components below.

The inverse tumor growth problem is given by the following minimization problem:

[TABLE]

subject to the reaction-diffusion forward problem $\mathcal{F}(\boldsymbol{p},\kappa,\rho)$:

[TABLE]

Here, we have introduced a parameterization for the tumor initial condition $c(0)=c(\boldsymbol{x},0)$ as $\boldsymbol{\phi}^{T}(\boldsymbol{x})\boldsymbol{p}$, where $\boldsymbol{p}$ is an $m$-dimensional vector, $\phi_{i}(\boldsymbol{x})=\phi_{i}(\boldsymbol{x}-\boldsymbol{x_{i}})$ is a radial function centered at point $\boldsymbol{x}_{i}$ (one of the grid points after spatial discretization), and $\boldsymbol{\phi}(\boldsymbol{x})=\{\phi_{i}(\boldsymbol{x})\}_{i=1}^{m}$, i.e.,

[TABLE]

This parameterization is mainly a modeling assumption. We assume that the spatial distribution of the tumor at its very early stages is given by $\phi_{i}$ for some $i$. There is no clear understanding on how tumors start, but the general consensus is that it starts at a single cell due to some mutation. At this resolution the continuum hypothesis breaks down and it is unclear how to couple such a single cell mutation to the PDE model. Instead, we use the PDE model after the tumor has grown sufficiently large while remaining localized. We label this time as zero and assume that the tumor volume fraction can be approximated well by $\phi(\boldsymbol{x}-\boldsymbol{x}_{i})$, where $\boldsymbol{x}_{i}$ is the center of mass of this tumor. However, we allow a small number of $\boldsymbol{x}_{i}$ points to be activated. This is to admit more complex shapes for the initial condition and also to compensate for the modeling errors in our tumor growth PDE.

The objective function minimizes the mismatch between the simulated tumor at $t=1$ and data on all observation points using the standard $L_{2}$ distance measure on square integrable functions. We balance this mismatch with an $L_{2}$ regularization on the inverted initial condition with regularization parameter $\beta$.

Similar to the works of [15], we observe all tumorous data points with concentration above a detection threshold ($c_{d}$). We define our observation operator $\boldsymbol{O}$ as follows:

[TABLE]

In our synthetic examples, we know the ground truth $c$. In clinical images, we use the segmentation probability maps as a proxy for $c$. Tab. 1 summarizes the notations. Next, we discuss the difficulty of inverting for $\boldsymbol{p}$, $\kappa$ and $\rho$ simultaneously. To do that, we use a simplified analytical model in one dimension.

2.2 Parameter inversion analysis for a simple linear analytical model

There have been several studies that outline the infeasibility of inverting for all parameters of a reaction-diffusion system. In particular, the authors in [33] note the non-uniqueness of the solution of the forward model with respect to different parameter configurations thereby rendering the inverse problem ill-posed. Here, we bolster this intuitive explanation through an analytical investigation of a simplified model which underpins our strategy in circumventing this issue.

We consider a 1D linear reaction-diffusion model with constant coefficients on the domain $\omega=[0,\pi]$ with periodic boundary conditions as our growth model. The forward model is given by:

[TABLE]

where $c(x,t)=c$ is the tumor volume concentration, $\rho$ is the reaction coefficient, and $\kappa$ is the diffusion coefficient. Our goal is to invert for $g(x)$, $\rho$ and $\kappa$. Notice that, for a moment here, we omit the parameterization we introduced in the previous section for $c(x,0)$ and we just try to invert for a field $g(x)$. We can solve Eq. (12) analytically using separation of variables and we obtain

[TABLE]

where $g_{n}$ are the spectral cosine coefficients of the initial condition $g(x)$. We truncate this infinite sum to $N$ terms, where we assume that the solution is fully resolved for sufficiently large $N$. For our inverse problem, we assume band-limited data given by: $d(x)=\sum_{n=0}^{N}d_{n}\cos(nx)$. Now, we can formulate the parameter inversion as the following optimization problem for $\{g_{n}\}_{n=0}^{N}$, $\rho$ and $\kappa$ as:

[TABLE]

The first order optimality conditions to Eq. (14) are:

[TABLE]

We observe that all equations are trivially satisfied when:

[TABLE]

Eq. (16) is informative in two ways. First, we see that even if $\rho$ and $\kappa$ are known, inverting for $g_{n}$ is exponentially ill-posed. Indeed, $g_{n}=e^{-\rho}e^{n^{2}\kappa}d_{n}$. Any perturbation in $d_{n}$ is amplified exponentially in the frequency and thus, in practice, $g_{n}$ cannot be identified for large $n$. But in our case, even small $n$ is problematic when $\rho$ and $\kappa$ are unknown. We cannot invert for $g_{n}$, $\rho$ and $k$ simultaneously since changing $\rho$ or $\kappa$ amounts to a scalar scaling of the initial condition.

This ill-posedness persists even if we re-parameterize the initial condition as before using $m$ degrees of freedom:

[TABLE]

where $\phi_{i}$ are some pre-selected basis functions and $\boldsymbol{p}\in\mathbf{R}^{m}$ is the (low dimensional) coefficient vector. Taking the Fourier decomposition of $\phi_{i}$, we have: $\phi_{i}(x)=\sum_{n=0}^{N}\phi_{in}\cos(nx)$, with spectral cosine coefficients $\phi_{in}$. Now we can write $g(x)=\sum_{n=0}^{N}(\boldsymbol{\phi}_{n}^{T}\boldsymbol{p})\cos(nx)$, with $\boldsymbol{\phi}_{n}=\{\phi_{in}\}_{i=1}^{m}$. Hence, the solution to the forward model is given by:

[TABLE]

with the optimization problem as:

[TABLE]

Similar to the non-parametric case, it is not possible to invert for the initial condition and $\rho$ simultaneously as any perturbation in $\rho$ simply scales the initial condition. In contrast to $\rho$, a perturbation in $\kappa$ cannot be absorbed into the initial condition anymore due to the frequency dependency of the factor $e^{-n^{2}\kappa}$.

In order to recover the reaction coefficient, we make the following modeling assumption: $\max_{x}(g(x))=1$. Notice that, since $c(x,t)$ and $g(x)$ are volume fractions, their values are between zero and one. But here we assert something stronger. We assume that $t=0$ corresponds to the first time the tumor volume concentration reaches 100% in some locations. This fixes the scale of the tumor initial condition and allows us to determine a reaction coefficient. Our second modeling assumption is to use a Gaussian function for $\phi_{i}(x)$ (see §3 for more details). Also, notice that the equation for the diffusion is better defined.

Hence, in summary, we can invert for $\kappa$ and the initial condition parameterization simultaneously, and can recover $\rho$ with an additional modeling assumption.

2.3 Optimality conditions for the full 3D problem

For the optimization, we additionally impose the bound constraints $\rho\geq 0$, $\kappa\geq 0$ and $\boldsymbol{p}\geq 0$ by definition of the forward problem. Hence, our formulation in Eq. (9) is now constrained by:

[TABLE]

Eq. (23) is difficult to impose formally. So we will omit it from the optimality conditions and we will use it when we discuss the numerical scheme. We introduce the Lagrangian and take variations with respect to $(c(\boldsymbol{x},t),\alpha(\boldsymbol{x},t),\boldsymbol{p})$, where $\alpha(\boldsymbol{x},t)$ is the adjoint variable of $c(\boldsymbol{x},t)$. This results in the following strong form of the first order optimality conditions.

[TABLE]

The basic solution strategy is to iterate in $(\boldsymbol{p},\kappa,\rho)$ variables using a bound constrained optimization algorithm in which at each iteration Eq. (24a, 24b) (forward) and Eq. (24c, 24d) (adjoint) are eliminated. Next, we describe the details of the numerical schemes for the forward and inverse problems.

3 Numerical Algorithms

First, we summarize the numerical scheme for solving for forward and adjoint variables given $\boldsymbol{p}$, $\rho$ and $\kappa$. We also give the precise definition of how we choose the number and centers of the Gaussian basis functions $\{\phi_{i}(\boldsymbol{x})\}_{i=1}^{m}$, which we use to parameterize $c(\boldsymbol{x},0)$. Then, we discuss the optimization algorithm and how we impose the sparsity and “$\max$” constraints on $\boldsymbol{\phi}^{T}\boldsymbol{p}$.

3.1 Forward solver numerics

We use a pseudo-spectral Fourier method on a regular mesh for spatial discretization and employ 3D fast Fourier transforms to compute all spatial differential operators. We enforce our boundary conditions on the surface of the brain and CSF boundaries using a penalty approach (see [15, 25]). We use periodic boundary conditions on $\Omega=[0,2\pi]^{3}$. We use an unconditionally stable, second order Strang operator-splitting method to numerically solve the forward and adjoint reaction-diffusion PDEs in Eq. (24). The forward equation is split into reaction and diffusion terms (for more details, see [15, 23]). We use the Crank-Nicolson [10] method for the diffusion term. This is an implicit method and the linear system is solved with a matrix-free preconditioned Conjugate Gradient method. We precondition the resulting implicit system of equations by solving them with constant coefficients computed as the average diffusion coefficient in the domain. The reaction term is solved analytically. The main cost of the forward and adjoint solves are fast Fourier transforms to apply derivatives and preconditioners.

3.2 Selection of basis functions for the tumor initial condition

Given a tumor segmentation from an MRI image, we proceed to define the Gaussians, i.e., their number $m$, their centers $\{\boldsymbol{x}_{i}\}_{i=1}^{m}$, and the standard deviation $\sigma$ of each Gaussian. To avoid spurious growth and diffusion of cancer cells in the CSF, we set the Gaussians to zero in this region. In other words, we have:

[TABLE]

where $\mathbb{1}$ is the indicator function. We construct a candidate set of Gaussians whose locations are placed at a subset of a regular grid with spacing $\delta=2\sigma$ so that we can control the conditioning of $M:=\int_{\Omega}\boldsymbol{\phi}\boldsymbol{\phi}^{T}d\boldsymbol{x}$.444We empirically determine the ratio $\delta/\sigma$ by observing the condition number of $M$. We set $\sigma=\lambda h_{\mathrm{image}}$, where $h_{\mathrm{image}}$ is the resolution of the input MRI image and $\lambda$ is an integer constant. Notice that $\sigma$ is a modeling parameter. If we refine the image resolution or the underlying discretization, $\sigma$ should remain unchanged.

From this candidate set, we select the Gaussians based on the amount of tumor cells (volume fraction) within a radius of one standard deviation around the candidate location $\boldsymbol{x}_{i}$. We choose a candidate Gaussian if the volume of tumor around it exceeds a user-defined tolerance $\tau_{V}$, i.e.,

[TABLE]

where $B^{\sigma}(\boldsymbol{x}_{i})$ is a ball of radius $\sigma$ around $\boldsymbol{x}_{i}$ and $d(\boldsymbol{x})$ is the tumor data. We show an illustration of this method in Fig. 2.

This adaptive approach is fully automatic and is effective, for example, in identifying multifocal tumorous regions without any manual seeding.

3.3 The $\ell_{1}$ regularized inverse problem

While an $L_{2}$ penalty is effective in regularizing the inverse problem for $\rho$, $\kappa$ and $\boldsymbol{p}$, it results in an initial condition that has an unrealistically large support and raises the question of how much we can trust the reconstructed $\rho$ and $\kappa$ coefficients. In fact, the reconstructed initial condition $c(\boldsymbol{x},0)$ closely resembles the data $d(\boldsymbol{x})$ up to a scaling, leading to $\kappa$ and $\rho$ being underestimated (see §4.3 for numerical experiments highlighting this issue). As discussed in the introduction, to truly restrict $c(\boldsymbol{x},0)$, we need another modeling assumption – that the tumor originated in a small number of locations. We call this property the sparsity of the initial condition.

**Enforcing sparsity. **The most common technique for imposing sparsity is to use the $\ell_{1}$ norm (see [63] for a comprehensive review of such techniques). Here, we assume that only a small number of Gaussians are activated and impose an $\ell_{1}$ penalty on the solution vector $\boldsymbol{p}$. We pursue an iterative recovery strategy similar to the algorithm described in [45]. The main difficulty in using this algorithm directly is that our formulation involves a non-linear and non-convex inverse problem unlike the operators for the mismatch term in [45]. This also means we lose any theoretical guarantees of finding a solution. In particular, the crucial step of identifying the important components of the signal relies on certain properties of the sampling matrix, such as the restricted isometry property, which may no longer hold. For our formulation, we use the gradient of the objective calculated using adjoints to identify the important components (similar to the original algorithm described in [45]). Before we outline our algorithm to reconstruct all the variables ($\boldsymbol{p},\kappa$ and $\rho$), we describe our strategy to enforce the non-convex and non-differentiable max constraint in Eq. (23) which allows us to identify $\rho$.

**Enforcing the max constraint. ** While there are methods such as relaxing the max constraint to inequalities and using slack variables to handle them or mixed integer programming formulations, we do not pursue them here, in interest of algorithmic simplicity. Instead, we use the following strategy:

(i) We start with an initial guess555We discuss our choice for this initial guess in §4.1. for $\rho$ and invert for the parameters $(\boldsymbol{p},\kappa)$. The resulting initial condition does not necessarily satisfy Eq. (23).

(ii) Then we re-scale the reconstructed initial condition appropriately to satisfy Eq. (23). This fixes the scale of the reaction coefficient, as discussed in §2.2.

(iii) Finally, we invert for the variables $(\rho,\kappa)$ keeping the (re-scaled) initial condition fixed.

Before describing the complete algorithm, we introduce some notation. We write

[TABLE]

for the support of $\boldsymbol{p}$, and with

[TABLE]

we denote the $s$* restriction of $\boldsymbol{p}$*, where $S$ is the set of indices of the $s$ largest-magnitude components of $\boldsymbol{p}$.

In a first phase, the algorithm simultaneously inverts for $\boldsymbol{p}$ and $\kappa$ following an iterative strategy:

1. (i)

*Identify significant Gaussians: *Compute the gradient using the current approximation for $\boldsymbol{p}$ and $\kappa$. The largest components of the gradient of $\boldsymbol{p}$ indicate the location of important Gaussians and form the new support of $\boldsymbol{p}$ (cf. Alg. 1, line 7). 2. (ii)

Merge supports: Merge the new support derived from the gradient with the support of the current approximation of $\boldsymbol{p}$ (cf. Alg. 1, line 8). 3. (iii)

Signal estimation: Solve an $L_{2}$ regularized inverse problem on the restricted subspace (defined by the merged support). This results in a new approximation for $\boldsymbol{p}$ and $\kappa$, (cf. Alg. 1, line 9). 4. (iv)

Signal restriction: Restrict the new approximation to $\boldsymbol{p}$ to the required sparsity using Eq. (28), (cf. Alg. 1, line 10). 5. (v)

Iterate until stopping criterion is true.

We terminate the compressive sampling when a maximum number of iterations has been reached or the relative change in objective function value is below a user-defined threshold $\tau_{J}$. This stopping criterion is typically used in $\ell_{1}$ minimization problems. Other criteria such as monitoring the monotonic decrease of the mismatch two-norm have been used in subspace pursuit algorithms but are valid only in the linear case with restricted isometry properties. Hence, we do not employ these rules.

In a second phase, the signal estimation step is performed once more (cf. Alg. 1, line 12) in the subspace defined by the support of $\boldsymbol{p}$ to get the final approximation (rather than just keeping the largest entries). This step is important as the reconstructed initial condition is used in the subsequent inversion for $\rho$: We re-scale the initial condition and solve the tumor inversion problem for parameters $\rho$ and $\kappa$ keeping the initial condition fixed. Next, we describe the details of the signal estimation step.

**The $L_{2}$ regularized inverse problem. **An algorithm for solving the $L_{2}$ regularized inverse problem (Eq. (9)) for $\boldsymbol{p}$ and $\kappa$ was presented in [15, 56]. Here, we summarize the main components for solving the system of optimality conditions in Eq. (24). We employ a quasi-Newton method, L-BFGS [47, p. 135ff]) with matrix-free approximations to the inverse Hessian and the Moré-Thuente criterion [43] for line-search. We enforce the bound constraints on our parameters $(\boldsymbol{p},\kappa,\rho)$ using projected gradients to approximate the quasi-Newton Hessian. We terminate our solver when the relative change in the norm of the gradient falls below a user-defined threshold $\tau_{g}$ or when a maximal number of quasi-Newton iterations is reached. We discuss the choice for these parameters in §4.1. We compute our reference gradient with a zero initial guess for $(\boldsymbol{p}^{0},\kappa^{0})$.

We present a pseudocode in Alg. 1.

4 Results

We test the performance of our solver on synthetically generated tumors and clinical datasets.

4.1 General setup

All 3D numerical experiments were executed on the Stampede2 system at the Texas Advanced Computing Center (TACC) at The University of Texas at Austin on its dual socket Intel Xeon Platinum 8160 (“Skylake”) nodes. Our solver is written in C++. We use MPI for parallelism compiled with Intel 17 compiler and 64 MPI tasks on 3 nodes for all our simulations. We use PETSc for all linear algebra operations, PETSc-TAO for the bound constrained L-BFGS optimization solver used in the $L_{2}$ regularized inverse problem, AccFFT for the 3D parallel FFTs and PnetCDF for I/O. We use a 2D pencil decomposition of our data for 3D distributed FFTs (see [14] for the parallel implementation details). For 1D test-cases, we use Matlab R2018b.

Performance measures. The primary goal of our numerical experiments is to demonstrate the ability of the solver to reconstruct the reaction and diffusion coefficients, and to reconstruct the initial condition with good localization accuracy. To evaluate the errors, we define the following error metrics to quantify our reconstruction:

1. (i)

relative error in $\kappa$ and $\rho$ reconstruction:

[TABLE]

where $\kappa^{\star}$ and $\rho^{\star}$ are the true coefficient values and $\kappa^{\text{rec}}$ and $\rho^{\text{rec}}$ are the reconstructed values, 2. (ii)

relative error in initial tumor reconstruction:

[TABLE]

where $c^{\text{rec}}(0)$ is the reconstructed initial condition and $c^{\star}(0)$ is the ground truth initial condition, 3. (iii)

relative error in final tumor reconstruction:

[TABLE]

where $c^{\text{rec}}(1)$ is the reconstructed tumor final concentration and $c^{\star}(1)$ is the final tumor concentration grown with the ground truth parameters $(\boldsymbol{p}^{\star},\kappa^{\star},\rho^{\star})$, 4. (iv)

relative $\ell_{1}$ error in initial tumor parameterization:

[TABLE]

with reconstructed parameterization vector $\boldsymbol{p}^{\text{rec}}$, 5. (v)

relative change in the norm of the gradient to help assess the convergence of our numerical scheme:

[TABLE]

where $g$ is the gradient at the final iteration and $g^{0}$ is the gradient with respect to the initial solution guess.

We evaluate the performance of two solvers, namely the $\mathbold{L_{2}}$ solver and the CS solver. The former refers to the $L_{2}$ (two-norm) regularized solver introduced in our previous work [15, 56] that does not use any sparsity constraints, while the latter is the compressive sampling solver described in Alg. 1.

Numerical parameters: We list all numerical parameters in Tab. 2 and outline any parameter change specific to a test-case in its respective section. We now provide a justification for our parameter choices:

1. (i)

For the selection of Gaussians, we choose a volume fraction $\tau_{V}$ of $0.99$ to ensure that only Gaussians occupied completely by tumorous cells within a neighbourhood of one standard deviation $\sigma$ get selected. We consider Gaussians with smaller volume fractions as unimportant candidates for the support of the tumor initial condition and hence we neglect them. 2. (ii)

We refer to [56] for the $L_{2}$ tumor regularization parameter $\beta$, which was chosen based on the data reconstruction quality and an L-curve analysis. For the compressive sampling (CS) solver, we set $\beta$ to zero as the inversion is performed for typically 10–20 parameters and we have experimentally found that no regularization is necessary. 3. (iii)

We pick the relative gradient convergence tolerance $\tau_{g}$ for all $L_{2}$ solves (Eq. (9)) to be $110-5$ in order to get sufficient accuracy in the reconstruction.666We note that our formulation is based on an optimize-then-discretize approach and hence our gradient is only accurate up to discretization errors. 4. (iv)

We set a maximal number of quasi-Newton iterations as 50 to prevent prohibitively large times-to-solution. But in most of our experiments, the target tolerance is reached within the given iteration budget. In the rare occasion this fails to happen (for example, in test-cases with extremely large tumors), we observe that additional quasi-Newton iterations do not lead to significant accuracy gains, but they do significantly increase the computational cost. 5. (v)

Based on numerical experimentation with different values of reaction and diffusion coefficients, we find $\rho=15$ to be a good representative reaction coefficient to capture quite “large” tumors ($\scriptstyle\sim$30% of white matter volume). Hence, we set this to be our initial guess for the reaction coefficient ($\rho^{0}$). 6. (vi)

The required sparsity level is a question that needs to be addressed empirically using a large study using a clinical dataset. Here, we hypothesize 5–10 activated sites for tumor initiation and verify numerically that, using this sparsity, we capture important activations in all our synthetic test-cases. For clinical studies, we plan to use an image-analysis based tool to estimate the focality of the tumor and then assign a number of sites per focal region. This is beyond the scope of this paper. 7. (vii)

In all our experiments, we observe inconsequential changes in sparsity after two compressive sampling iterations and hence set it as an upper bound to the number of iterations for the $\ell_{1}$ sparsity.

We examine different types of test-cases to illustrate the performance of our model. We describe the setup of each experiment and the corresponding results in the following sub-sections.

4.2 1D constant coefficients test-case

We start with a simple 1D example to demonstrate that the CS solver outperforms the $\mathbold{L_{2}}$ solver.

*Setup. *The data $d(\boldsymbol{x})$ is generated using the non-linear forward model with constant reaction and diffusion coefficients. We use a spatial resolution of $256$ grid points in $\Omega=[0,2\pi]$ and a temporal resolution of 100 time steps within a time horizon of 1. We consider two variations for our test-cases:

[TABLE]

We consider more variations for our 3D test-cases to capture the heterogeneity in tumor shapes generally observed in clinical data. We report our performance measures, relative gradients and solver timings in Tab. 3 and present qualitative results of the reconstruction in Fig. 3.

*Observations. *While both solvers are able to fit the data satisfactorily (reconstruction errors less than $5$%), we observe that the compressive sampling (CS) solver outperforms its $\mathbold{L_{2}}$ counterpart in biophysical parameter estimation. This is evident in Fig. 3a and Fig. 3c, where the reconstructed initial conditions resemble grown tumors compared to the sparse ground truth. This can also be observed in the $\ell_{1}$ relative error norms of the initial condition parameterization ($\mu_{\boldsymbol{p},\ell_{1}}$) which show about an order of magnitude difference between the two solvers. Here, we would like to emphasize that it is because of the unknown reaction coefficient that the $\mathbold{L_{2}}$ solver reconstructs a rich initial condition, and consequently a low (and inaccurate) reaction coefficient. Hence, the tumor does not grow much, and instead diffuses to fit the data – this is seen in the reconstructed values for the reaction and diffusion coefficients. As we shall see, this issue is mitigated a little in 3D, but still persists. We also note that the initial conditions might not be perfectly recoverable due to the ill-conditioned nature of the inverse problem. For instance, the multi-focal test-case (1D-C2) finds activations in neighbouring Gaussians to the ground truth. This localization problem is unavoidable due to ill-conditioning. If the time horizon is so large that the maximum tumor concentration is everywhere in the spatial domain, it is impossible to recover any information about the initial location of the tumor. So, the larger the tumor (or equivalently, the larger the time horizon), the harder it is to localize the initial condition.

4.3 Artificial tumor (AT) test-case

In this test-case, the data $d(\boldsymbol{x})$ is generated using the non-linear forward model with inhomogeneous coefficients to capture the heterogeneity of tumors. These experiments are intended as a proof-of-concept of the performance of our solver on 3D synthetic datasets.

*Setup. *We grow synthetic tumors in the segmentation of a statistical brain atlas [8] MRI scan (an atlas obtain by averaging several MR images). The segmentation labels are gray matter, white matter, and cerebrospinal fluid-filled ventricles (CSF). We use sparse initial conditions to grow the tumor to a significant size resembling clinical observations and invert for all biophysical parameters using the grown tumor at $t=1$ as input data to the solver. For all our test cases, we assume that the tumor grows and diffuses only in white matter. For this test-case, we observe the data everywhere, i.e., $c_{d}=0$. We consider the following variations:

[TABLE]

We report our quantitative results (performance measures, convergence, and solver timings) in Tab. 4 and visualize the reconstruction to qualitatively assess the performance of our solvers in Fig. 4 - Fig. 7. Each figure shows an isometric view of the patient brain with 3D volumes of the grown tumor data, ground truth, and reconstructed initial condition. We additionally visualize the ventricles and sections of healthy tissue to reveal the heterogeneity of the brain (and hence the tumor).

*Observations. *We observe that for the medium sized tumor (AT-C1) the reconstruction is nearly perfect for the CS solver while the $\mathbold{L_{2}}$ solver produces a lower performance, particularly in the inversion of the reaction coefficient. This is primarily due to the rich initial conditions produced by the $\mathbold{L_{2}}$ solver which impedes its ability to predict the correct reaction scaling using our method. We note that if the reaction scaling is known beforehand, the $\mathbold{L_{2}}$ solver can potentially have better performance (see [56, 16] for similar synthetic experiments). This problem is magnified for larger tumors (AT-C2 and AT-C3) where the predicted reaction coefficient shows about $66\%$ and $54\%$ relative error in the $\mathbold{L_{2}}$ solver for the two test-cases respectively (as compared to around $15\%$ and $1\%$ error with sparsity constraints). The reconstruction in diffusivity also suffers similarly without enforcing sparsity. The final test-case shows a multi-focal tumor with two far apart proliferation sites. While both the solvers produce disjoint initial conditions for the tumor, the CS solver does better in reconstructing the model parameters: 2% versus 37% relative error in diffusivity and 9% versus 50% relative error in reaction coefficient. Both solvers show similar performance in their ability to fit the data (as measured by the relative error in final tumor reconstruction). Hence, our synthetic observations show consistent improvement in target parameter estimation for all test-cases with sparse initial conditions, without sacrificing the final tumor reconstruction quality.

4.4 Artificial tumor with observation operator test-case

As mentioned before, an observation operator is necessary since the full extent of tumor invasion is not detectable using existing imaging technologies. The purpose of this experiment is to investigate the effect of the observation operator $\boldsymbol{O}$ on the model parameter reconstruction abilities of the solver. In particular, we would like to see if partial observations significantly deteriorate the reconstruction quality.

*Setup. *We assess the performance and sensitivity of our solver by choosing different values for the detection threshold $c_{d}\in\{0,0.1,0.2,0.3,0.4,0.5\}$, which defines the observation points according to Eq. (11). We define the parameter $n_{d}$ to quantify the number of voxels observed with each threshold value. We choose the synthetic test-case AT-C2 as a suitable representation of a large and significantly invasive (aggressive) tumor (large parameter values: $\rho^{\star}=12$, $\kappa^{\star}=$510-2$$). We report all our performance metrics for the different observation thresholds in Tab. 5. We emphasize that the final tumor reconstruction error is with respect to the fully grown tumor with the ground truth model parameters and not the data restricted to the observation points. This definition is more useful in helping us understand if the model is reconstructing the full extent of the tumor, including the areas which are typically invisible in imaging data.

*Observations. *The most important observation is that the solver is able to approximately establish the extent of tumor invasion up to a detection threshold of $c_{d}=0.3$. This is reflected in the performance metrics that do not register significant drops over these threshold values. The final tumor reconstruction error shows perturbations of up to 4% till this detection threshold. The values for the inverted model parameters also do not exhibit appreciable degradation and roughly remain around the ground truth coefficients – $9.8$ to $11.1$ for the reaction coefficient and $0.0389$ to $0.06$ for the diffusion coefficient. Larger detection thresholds (from $c_{d}=0.4$ onwards) lead to worsening reconstruction errors and significant deterioration in target parameter estimation. This is expected since there is insufficient data for any reliable inversion, especially since this tumor has substantial diffusion resulting in considerable loss of data when choosing higher values for the detection threshold. Finally, while there is no consensus regarding the value of $c_{d}$ in literature (generally values around 0.2 - 0.4 have been used [33, 61]), this value can easily be modified in our formulation once more specific information becomes available.

4.5 Artificial tumor with varying spatial resolution test-case

The ground truth for parameter estimation should ideally be independent of the numerical discretization used. Since we use a pseudo-spectral discretization scheme which demands smooth functions (to avoid aliasing errors), a mesh-dependent numerical diffusion is introduced in our formulation. The purpose of this test-case is to assess the degree of information loss due to spatial discretization and to understand the extent of reliability of our inverted parameters with mesh size.

*Setup: *We consider the synthetic test-case AT-C2 as a suitable representation of a large and invasive tumor. Since there is no analytical solution to the forward problem, we grow a synthetic tumor in a high spatial resolution of $512^{3}$ and assume this to be the “true” tumor. We sub-sample this tumor to resolutions $128^{3}$ and $256^{3}$ and use the resulting tumors as data for inversion in the respective spatial discretizations. The initial radial basis $\sigma$ is kept constant across all grids (set to the spacing of the $128^{3}$ grid). We consider the following test cases:

[TABLE]

We use different orders of magnitude for the diffusion coefficient to gauge the extent of numerical diffusion, relative to actual physical diffusion. We report our performance measures in Tab. 6.

*Observations: *We observe that we cannot recover diffusion coefficients smaller than around $110-3$ using a spatial resolution of $128^{3}$. This can be seen in the values for $\kappa^{\text{rec}}$ (diffusivity reconstruction) which fail to recover any smaller values (test-case C1; C2). As a result, the values for $\rho^{\text{rec}}$ (reconstructed reaction magnitude) also improve for test-cases with higher ground truth diffusivity values as the effect of mesh-dependent diffusion is reduced. The higher resolution mesh ($256^{3}$) is able to recover diffusion coefficient values up to an order of magnitude of $110-5$. Consequently, we observe larger differences for the $\rho^{\text{rec}}$ values between the two resolutions as the ground truth diffusivity drops. The higher resolution inversion does a better job in estimating the diffusivity (due to lower mesh-diffusion errors) and hence can estimate the reaction coefficient reasonably where the lower resolution inversion performs poorly. Hence, our estimated parameters are well founded only if they are not in the range of mesh-induced numerical diffusion. One way to circumvent this error is to perform grid continuation where we invert for parameters in successively higher mesh resolutions using the sparsity and inverted parameters from the coarser meshes as an initial guess and settle on a resolution as the results of the inversion stagnate.

4.6 Artificial tumor with noise test-case

The purpose of this test-case is to examine the robustness of our inversion against noise in the data. Instead of using high frequency uncorrelated noise (white noise), we are interested in observing the effects of a more realistic low-frequency noise model.

*Setup. *We approximate the source of noise by modifying the low-frequency spectrum of our data. We add a uniform random noise scaled by the inverse squared frequency to the respective frequency components of the data . Then, we take the inverse Fourier transform of this noisy spectrum to obtain our noise-corrupted data. We use AccFFT to compute our frequency spectrum and IFFTs in a distributed manner. We choose our ground truth parameters from the artificial tumor test-case AT-C2 since this test-case features the largest invasive synthetic tumor making the inversion most ill-conditioned. We add various levels of noise (as measured by the relative error in the true data and noise-corrupted data) and report our findings in Tab. 7.

*Observations. *We observe that our model does not exhibit any significant sensitivity for this noise model. We consider the small changes in performance and convergence metrics between the different noise levels inconsequential and attribute them to the inherent randomness in generating the data.

4.7 Real tumor (RT) test-case

The purpose of this test-case is to evaluate our solver on real patient brains diagnosed with glioblastoma. While we have no knowledge of the ground truth model parameters for real data, this experiment attempts to demonstrate the similarity in reconstruction trends observed in the artificial tumor test-cases (as designed above) and a clinical dataset.

*Setup: *Since we do not have access to the healthy patient brain, we cannot resort to direct inversion of model parameters. There have been many approaches in literature (see [16, 56] for a comprehensive review) to approximate the healthy brain through methods such as diffeomorphic image registration and registration-tumor coupled optimization formulations. While such techniques are necessary to capture the healthy patient brain, we do not pursue them since this test-case is intended to primarily serve as an illustration. Instead, we simply replace all tumor segmented voxels in the patient image with white matter and invert only for the diffusivity and reaction in the white matter. We note that by doing so we have lost any information concerning the inhomogeneity of our model coefficients in white and gray matter. We use MRI scans from the BraTS dataset [40]. These images are annotated with voxel-wise segmentation of tumor acquired through manual delineation performed by expert radiologists. For the healthy tissue segmentation, we refer to [17] where the brain is further segmented into gray matter, white matter and cerebrospinal fluid-filled ventricles through diffeomorphic image registration with healthy atlases. Since our model incorporates only a single species of tumor cells, we combine the enhancing and necrotic tumor sub-classes into one tumor class (also known as the tumor core) and use this as our tumorous species. We report the inverted values for our model coefficients, data reconstruction error, and convergence metrics in Tab. 8. We visualize our reconstruction in Fig. 8.

*Observations: *We observe sparse initial condition reconstruction with the CS solver and consequently larger diffusion and reaction coefficients for the tumor model. This is in accordance with our observations in the synthetic tumor test-cases. Both solvers produce roughly equivalent reconstruction errors for the final tumor condition which is around 16-19%. We note that if our primary goal is to reduce this mismatch we might fare better with more complex initial conditions which resemble fully grown tumors. At the extreme scenario, pure interpolation of the data with radial basis, i.e., setting all model coefficients to zero, would lead to a significantly better mismatch. However, this incurs the severe penalty of losing all information about the biophysical model at hand. Since the focus of this work is on the latter, we accept higher reconstruction errors as long as we can identify the target model parameters reliably. This is ultimately more useful for tumor growth understanding and prediction. We emphasize that approaches like a coupled image registration and tumor inversion formulation are necessary for access to a more informative healthy patient (see [56] for more details). Our plan is to conduct a systematic study of the solver considering hundreds of cases. However, such a study is beyond the scope of this paper and requires a separate treatment.

5 Conclusions

We presented an inverse problem formulation for the simultaneous estimation of all parameters in a commonly used reaction-diffusion model, namely the reaction coefficient, diffusion coefficient, and the spatial distribution of the tumor initial condition. The data that drive the inversion come from MRI scans, in particular the segmentation of tumor regions. We motivated the need for the sparse localization of the initial tumor in order to reliably reconstruct the other parameters. Additionally, this constraint is biophysically motivated and hence can be of clinical relevance in understanding the origin of the tumor and its role in the subsequent evolution of the disease. We presented results demonstrating the performance of our solver in reconstructing the ground truth model parameters, even for large and complex synthetic tumors. The results also highlight significant improvement in reconstruction quality with the use of sparsity constraints. We observe similar reconstruction trends in real patient data, thus providing a proof-of-concept of applicability in clinical practice. Finally, we employed novel reduced-space optimization algorithms and compressive sampling to solve our non-linear inverse problem and have tested its robustness against partial data observations and noise. Our software is implemented in a scalable manner enabling realistic time-to-solutions, a significant factor for applications in a diagnostic setting. While limitations like the necessity of more complex growth models exist, we envision this work as an important stepping stone towards the personalization of tumor growth models.

Our next step is to conduct a comprehensive evaluation of the model using a large number of clinical images. Our goal with such a study will be to correlate the model parameters with clinical outcome, for example survival. There are several ways that this work can be extended. Even if multiple time-snapshot scans exist, the presence of partial observations still poses difficulties. The same framework we discussed here can be used, but with multiple scans the problems of inverting for $\rho$ and $\kappa$ and the tumor initiation can be possibly decoupled. A more significant challenge is incorporating more complex biophysical models, while still having a single scan. All these directions are ongoing work in our group.

Acknowledgments

This material is based upon work supported by NIH award 5R01NS042645-14; by NSF awards CCF-1817048 and CCF-1725743; by the U.S. Department of Energy, Office of Science, Office of Advanced Scientific Computing Research, Applied Mathematics program under Award Number DE-SC0019393; and by the U.S. Air Force Office of Scientific Research award FA9550-17-1-0190. Any opinions, findings, and conclusions or recommendations expressed herein are those of the authors and do not necessarily reflect the views of the AFOSR, DOE, NIH, and NSF. Computing time on the Texas Advanced Computing Centers Stampede system was provided by an allocation from TACC and the NSF.