Inverse Problem Regularization for 3D Multi‐Species Tumor Growth Models

TL;DR

The paper introduces a multi-species tumor growth model for brain cancer and a method to calibrate it using MRI scans, improving tumor estimation accuracy.

Contribution

A novel two-stage regularization approach for calibrating a multi-species tumor growth model from single-time MRI data.

Findings

The proposed regularization improves tumor Dice score by 5%–10% compared to single-species reconstructions.

Model parameter reconstruction errors are reduced by 4%–80% with regularization.

The model can estimate infiltrative tumor cells using observable tumor species.

Abstract

We present a multi‐species partial differential equation (PDE) model for tumor growth and an algorithm for calibrating the model from magnetic resonance imaging (MRI) scans. The model is designed for glioblastoma multiforme (GBM) a fast‐growing type of brain cancer. The modeled species correspond to proliferative, infiltrative, and necrotic tumor cells. The model calibration is formulated as an inverse problem and solved by a PDE‐constrained optimization method. The data that drives the calibration is derived by a single multi‐parametric MRI image. This is a typical clinical scenario for GBMs. The unknown parameters that need to be calibrated from data include 10 scalar parameters and the infinite dimensional initial condition (IC) for proliferative tumor cells. This inverse problem is highly ill‐posed as we try to calibrate a nonlinear dynamical system from data taken at a single time. To address this ill‐posedness, we split the inversion into two stages. First, we regularize the IC reconstruction by solving a single‐species compressed sensing problem. Then, using the IC reconstruction, we invert for model parameters using a weighted regularization term. We construct the regularization term by using auxiliary 1D inverse problems. We apply our proposed scheme to clinical data. We compare our algorithm with single‐species reconstruction and unregularized reconstructions. Our scheme enables the stable estimation of non‐observable species and quantification of infiltrative tumor cells. Our regularization improves the tumor Dice score by 5%–10% compared to single‐species model reconstruction. Also, our regularization reduces model parameter reconstruction errors by 4%–80% in cases with known initial condition and brain anatomy compared to cases without regularization. Importantly, our model can estimate infiltrative tumor cells using observable tumor species.

We present a multi‐species partial differential equation (PDE) model for tumor growth and an algorithm for calibrating the model from magnetic resonance imaging (MRI) scans. The model is designed for glioblastoma multiforme (GBM) a fast‐growing type of brain cancer.

Figures

A schematic representation of the multi‐species model dynamics in 1D. The tumor species growth process is represented, with vascularization indicated by oxygen.

Ill‐posedness analysis of multi‐species model in 1D. The spectrum of Hessian eigenvalues for 3500 parameter combinations is presented in Figure 2a. The singular values calculated in Algorithm 1 are shown in Figure 2b. We depict the absolute value of last four directions of the computed regularization term (UR) in Figure 2c corresponding to the four smallest singular values. The exponential decrease of eigenvalues in Figure 2a indicates high ill‐posedness of the problem. Additionally, the wide range of small eigenvalues suggests varying degrees of ill‐posedness among samples, which is accounted for in Algorithm 1. Figure 2b demonstrates how different directions impact the ill‐posedness, with smaller singular values indicating greater ill‐posedness. Thus, using the inverse of the singular directions can lead to improved inversion results. The results in Figure 2c clearly indicates ill‐posedness for high values γ0 and a combination of high values of ρ and high values of δc contributes to ill‐posedness.

1D inversion results of synthetic data with additive Gaussian noise (5%,10%,20%,30%). Here we show the reconstructed species concentration at T=1. First to fourth rows show relative noise levels. Plots in first to third columns show proliferative, necrotic, and infiltrative species, respectively. The green, gray, magenta, and blue lines represent the underlying ground truth species (before observation denoted as “Ground Truth”), ground truth species with noise (“Noisy Ground Truth”), reconstructed species without regularization (“Non‐Reg. Inv.”), and reconstructed species with regularization (“Reg. Inv.”), respectively. Regularized inversion is observed to better match the underlying species even when the data is perturbed by noise, particularly in the case of infiltrative cells where deviations are greatly reduced with regularization.

Evaluating 3D inversion using synthetic data with known IC and brain anatomy. We start by generating a tumor test‐case and showing tumor species including proliferative or enhancing (EN), necrotic (NEC) and edema (ED) and segmentation in the first row (a). We then add 5%, 10%, 20%, and 30% Gaussian noise to the species and apply an observation operator to get noisy data. We perform inversion with and without regularization on this data. The inversion results are displayed in rows 2–4, with the enhancing, necrotic, and infiltrative concentrations in columns 1–3, and the segmentation in the last column. The regularized inversion yields acceptable segmentation for all noise levels, whereas the non‐regularized case shows more variability for different noise levels.

3D inversion result with synthetic test‐cases generated with an unknown IC and brain anatomy. Different normal brains with no mass deformation are used to generate the synthetic cases. In each row, we present a single test‐case, with tumor species (proliferative or enhancing (EN), necrotic (NEC), edema (ED) and infiltrative cells) from the synthetic cases and inversion displayed in the first, second, and third columns, respectively. The given data and segmentation computed by the inversion algorithm are depicted in the fourth column. Moreover, we depicted the projected ground truth IC and estimated IC in the segmentation for data and inversion, respectively, Although the IC is an estimate, the reconstruction from the inversion is good as evidenced by the qualitative observation. The challenge of reconstructing multi‐species arises due to the unknown initial brain anatomy and isotropic migration of tumor in all directions due to assigned white matter label.

We conduct a multi‐species inversion of the BraTS20_Training clinical dataset, selecting two patients, using the MRI scan segmentation as input. Our results are presented in a tabular format, with each row showing the inversion and data for each patient. The first three columns depict the inverted species (proliferative or enhancing (EN), necrotic (NEC), and infiltrative cells), the fourth column shows the observed species from the inversion, and the last column displays the data. Our results demonstrate that, despite only having a single snapshot from each patient, we can qualitatively reconstruct the species distribution with reasonable accuracy.

Tables

TABLE 1. Common notations used to describe the forward and inverse problem for the multi‐species tumor growth model.

Notation	Description	Range
x	Spatial coordinates	(02π3 for 3D and 02π in 1D)
t	Time	0,1
p	Proliferative tumor cells	0,1
i	Infiltrative tumor cells	0,1
n	Necrotic tumor cells	0,1
c	Tumor cell (c=i+p+n)	0,1
o	Oxygen concentration	0,1
g	Gray matter cells	0,1
w	White matter cells	0,1
f	CSF/VT	0,1
p0	Initial condition (IC) for proliferative cells (px,0) (See Equation (1g))	0,1
ℛ	Tumor growth operator (logistic growth operator)	—
D	Tumor migration operator (diffusion)	—
ρm	Inhomogeneous reaction rate	—
k	Inhomogeneous diffusion rate	—
αo	Transition function from p to i cells	—
βo	Transition function from i to p cells	—
γo	Transition function from i and p to n cells	—
ρ	Reaction coefficient (See Equation (2))	5.0,2.50×101
κ	Diffusion coefficient (See Equation (4))	1×10−3,1×10−1
α0	Transitioning coefficient from p to i cells (See Equation (6))	1×10−1,1.00×101
β0	Transitioning coefficient from i to p cells (See Equation (7))	1×10−1,1.50×101
γ0	Transitioning coefficient from i and p to n cells (See Equation (8))	1.0,2.00×101
δc	Oxygen consumption rate (See Equation (1d))	1.0,2.00×101
δs	Oxygen supply rate (See Equation (1d))	1.0, 8.0
oh	Hypoxic oxygen threshold (See Equations (2) and (8))	1×10−3,8×10−1
oi	Invasive oxygen threshold (See Equations (2), (6), and (7))	2×10−1,1.0
ie	Threshold coefficient for edema (See Equation (12))	1×10−3,3×10−1

TABLE 2. Numerical parameters for tumor inversion and forward solver. Initial guess is mid‐range for growth coefficients. We refer to [3] for IC inversion parameters in more details.

Parameter	Value
Spatial discretization for 1D inversion	256
Spatial discretization for 3D inversion	1603
Time‐step of the forward solvers in 1D and 3D	2×10−2
Shape factor ω in Oc, On,Op,α,β and γ	32
Shape factor ω in ℋi−ie	256
Epsilon used in finite differencing of Hessian, ϵ	1×10−6
Regularization parameter in 1D inversion, λ	1×10−4
Regularization parameter in 3D inversion, λ	1×10−5
Relative function tolerance	1×10−7
Absolute function tolerance	1×10−3
Sample size per iteration for the CMA‐ES solver	16

TABLE 3. Evaluation of 1D inversion scheme for the model parameters with synthetic data assuming known initial condition for the proliferative tumor cells.

Case	Noise (%)	πn	πp	πl	μn,L2	μp,L2	μi,L2	eq
Non‐Reg	5	1.65×10−2	9.36×10−1	4.07×10−2	7.74×10−2	2.43×10−1	5.99×10−1	6.64×10−1
Reg	5	9.86×10−1	9.36×10−1	5.80×10−1	8.01×10−2	1.45×10−1	3.69×10−1	5.23×10−1
Non‐Reg	10	1.77×10−2	9.34×10−1	9.89×10−2	4.91×10−2	1.11×10−1	8.73×10−1	1.36×100
Reg	10	9.83×10−1	9.29×10−1	4.26×10−1	5.79×10−2	8.12×10−2	7.59×10−1	6.86×10−1
Non‐Reg	20	1.50×10−2	8.92×10−1	5.59×10−2	6.04×10−2	2.89×10−1	8.07×10−1	6.76×10−1
Reg	20	9.70×10−1	8.99×10−1	5.13×10−1	9.86×10−2	3.25×10−1	6.87×10−1	7.98×10−1
Non‐Reg	30	6.85×10−2	8.22×10−1	4.43×10−2	2.41×10−1	2.64×10−1	4.96×10−1	4.65×10−1
Reg	30	9.17×10−1	8.32×10−1	3.71×10−1	3.17×10−1	2.56×10−1	4.14×10−1	4.74×10−1

TABLE 4. Numerical evaluation of 3D inversion using synthetic data with known anatomy and IC.

Case	Noise (%)	πn	πp	πl	μn,L2	μp,L2	μi,L2	eq
Non‐Reg	5	1.98×10−1	5.72×10−1	5.87×10−1	9.95×10−1	9.36×10−1	4.21×10−1	1.30×100
Reg	5	9.70×10−1	9.81×10−1	9.84×10−1	4.41×10−2	5.66×10−2	2.91×10−2	2.57×10−1
Non‐Reg	10	9.23×10−1	9.53×10−1	9.52×10−1	2.14×10−1	2.51×10−1	1.07×10−1	7.58×10−1
Reg	10	9.35×10−1	9.58×10−1	9.61×10−1	2.07×10−1	2.83×10−1	8.90×10−2	6.69×10−1
Non‐Reg	20	9.57×10−1	9.72×10−1	9.73×10−1	1.01×10−1	1.44×10−1	1.39×10−1	9.42×10−1
Reg	20	9.49×10−1	9.67×10−1	9.67×10−1	1.18×10−1	1.36×10−1	1.21×10−1	4.81×10−1
Non‐Reg	30	9.24×10−1	9.53×10−1	9.44×10−1	1.94×10−1	3.40×10−1	1.42×10−1	3.76×10−1
Reg	30	9.54×10−1	9.71×10−1	9.70×10−1	1.06×10−1	1.83×10−1	1.13×10−1	3.30×10−1

TABLE 5. Quantitative evaluation of 3D inversion with unknown IC using synthetic data.

Case	πtcms	πtcss	πp	πn	πl	μn,L2	μp,L2	μi,L2	μ0,L2	μ0,ℓ1	μ0,ℓ1x	μ0,ℓ1y	μ0,ℓ1z	eq
Case 1	8.43×10−1	7.73×10−1	6.46×10−1	7.59×10−1	7.85×10−1	9.78×10−1	5.08×10−1	9.00×10−1	1.18×100	1.60×100	1.14×100	1.12×100	1.05×100	4.28×100
Case 2	7.95×10−1	7.35×10−1	5.43×10−1	7.74×10−1	5.74×10−1	1.28×100	8.78×10−1	9.48×10−1	1.10×100	1.33×100	1.06×100	6.82×10−1	1.06×100	1.36×100
Case 3	8.19×10−1	6.78×10−1	6.13×10−1	7.93×10−1	7.20×10−1	1.26×100	7.64×10−1	8.94×10−1	1.09×100	1.27×100	1.06×100	1.04×100	1.07×100	8.28×10−1

TABLE 6. The results of the multi‐species model inversion using clinical data from the BraTS20 dataset: The performance of the model is measured by the Dice score, which serves as a quantitative measure of the accuracy of the reconstructed tumor species.

Patient's ID	πtcms	πtcss	πn	πp	πl
BraTS20_Training_039	8.23×10−1	8.37×10−1	6.94×10−1	6.18×10−1	4.24×10−1
BraTS20_Training_042	9.01×10−1	8.92×10−1	7.25×10−1	5.41×10−1	0.00×100

Full text

Introduction

1

Biophysical models of tumor growth can help quantify tumor infiltration that is not clearly visible. Quantifying infiltration is used in patient stratification, preoperative planning, and treatment planning. As every subject's anatomy and disease is unique, a main challenge is constructing subject‐specific biophysical models. This requires fitting free biophysical model parameters to data. Such parameters include the healthy brain anatomy, the location in which the tumor started, and several scalar growth coefficients modeling spreading, proliferation, and clearance. In this paper, we present an inversion methodology that integrates classical inverse problem theory with a multi‐species brain tumor growth model. Our approach addresses some challenges of single time snapshot inversion and our numerical results show strengths and weaknesses of the scheme.

Contributions

1.1

In our previous study [1], we developed a multi‐species tumor growth model integrated with tumor mass deformation but we did not consider its calibration from data. Calibrating this model is challenging due to the complexity of non‐linear coupled partial differential equations (PDEs) and the sparsity of the data. We are interested in reconstruction using a multiparameteric MRI scan at a single time point. We split the reconstruction into two stages. In the first stage we reconstruct the tumor initial condition (IC) assuming a simpler single species model. In the second stage, we reconstruct 10 scalar PDE coefficients using the IC from the first stage. For the first task we use our previous work on inverting for the IC using a single‐species PDE [2, 3] and use it to reconstruct model parameters. Our focus in this paper is two‐fold: (1) analyze and derive algorithms for the second stage inverse problem; and (2) evaluate the combined inversion using synthetic data. We summarize our contributions below.

1. Using a 1D model, we analyze the Hessian of the objective function (numerically), and discuss the severe ill‐posedness of the inversion problem for the 10‐parameter multi‐species PDE tumor growth model. We also use the 1D formulation to construct a weighted regularization term that we then use for 3D reconstruction. Since we do not have the ground truth parameters for clinical data, we evaluate this scheme with synthetic data. We empirically demonstrate that this regularization term improves growth coefficient estimation in the inversion setting.
2. We then estimate the tumor's IC using our previous works [2, 3]. We combine IC inversion with growth coefficients inversion. We evaluate this full scenario on synthetic data and highlight the strengths and weaknesses of the method.
3. We parallelize our execution to efficiently solve the forward problem using parallel architecture and graphics processing units (GPUs), allowing for 3D inversion on a 1603 imaging resolution to take approximately 5 h for any synthetic case on 4 GPUs.

As an example, we also report results of the reconstruction using a clinical dataset. A full clinical evaluation of the overall methodology and its comparison with a single‐species model is a subject of ongoing investigation and will be reported elsewhere.

Limitations

1.2

The proposed multi‐species model includes edema because it is clearly visible in multi‐parametric magnetic resonance imaging (mpMRI) scans. Our edema model is a simple algebraic relation, similar to our previously proposed model in [1]. In this analysis, we also ignore the so‐called mass effect, the normal brain tissue mechanical deformation due to the presence of the tumor. We have designed schemes to account for mechanical deformation in single‐species models ([4, 5]), but this is beyond the scope of this study. Our growth model does not incorporate diffusion tensor imaging, which can help capture more complex invasion patterns. Our single‐snapshot multiparametric MRI does not include more advanced imaging modalities like perfusion. In this paper, we use a simplified 1D model to compute the regularization, which differs from the 3D model in some terms. We estimate the tumor IC using a single‐species reaction diffusion model proposed in our previous studies [2, 3] and do not invert for the IC using the multi‐species model. Our scheme is primarily focused on brain tumors, but it is also applicable to tumors in other organs such as the pancreas, breast, prostate, and kidney.

Related Works

1.3

The modeling of tumor growth studies can be divided into two categories: 1. Forward models 2. Inverse problems and solution algorithms. Forward models mainly focus on the biophysical model that describes the biological process of brain tumors. Inverse problems and solution algorithms, on the other hand, involve aspects such as noise models, regularization, observation operators, and calibration of forward models in order to find the parameters to reconstruct data. Despite numerous attempts to simulate tumor growth models, less work has been done on the inverse problems themselves.

The most commonly used forward model for tumor growth is the single‐species reaction–diffusion model, as seen in various studies [6, 7, 8, 9]. Other models incorporating mass effect, chemotaxis, and angiogenesis are also explored in [10, 11, 12]. Despite offering insight into biological processes, these models lack calibration due to mathematical and computational complexities. Our group has worked to address these challenges, enabling clinical analysis of the models [5, 13].

Inverse problems have been extensively studied to calibrate single‐species reaction–diffusion models and quantify infiltrative tumor cells [14, 15, 16]. In reference [17], the inverse problem is formulated by including anisotropy diffusion derived from diffusion tensor imaging. In reference [18] the authors utilize Bayesian inversion for personalized inversion of model parameters for invasive brain tumors. Our previous works [2, 3] addressed the challenges of calibrating the growth coefficients and IC inversion simultaneously, proposing a ℓ0 constraint to localize IC inversion. In reference [4], the model is combined with mass deformation and clinical assessment is performed for a large number of MRI scans [5], improving survival prediction for patients. However, the single‐species model does not use all available information from MRI scans and considers necrotic and enhancing cells as single‐species while ignoring the edema cells, ignoring the biological process between species and resulting in a model output that differs from observed segmented tumor from MRI scans. Recent studies have incorporated data‐driven methods with follow‐up MRI validation [19, 20]. While these methods leverage larger datasets, they do not consider the coupled multi‐physics nature of tumor growth. Our approach focuses on modeling the complete multi‐physics system, though with a more limited patient cohort.

Recently [21], a dual‐species model using longitudinal MRI scans to calibrate the growth coefficients has been developed. However, longitudinal data is rarely available in clinical practice. In our previous work [1], we formulated a multi‐species model from [11] and combined it with mass deformation of the brain as a forward solver. The solver developed in this study was applied to analyze multiple clinical cases, as detailed in [22]. The present paper, however, focuses on addressing the inverse problem, with particular emphasis on examining its ill‐posed nature. Furthermore, we propose a regularization scheme to mitigate the challenges associated with ill‐posedness, an aspect not explored in our previous work [22]. This investigation aims to enhance the robustness and reliability of the solver in clinical applications.

Outline

1.4

In Section 2, we outline the forward problem and mathematical formulation of the inverse problem. We then present the 1D model, derive the adjoint and gradient formulations, and estimate the Hessian of the problem. We then propose an algorithm to compute a weighted regularization operator, which we then use in 3D reconstructions. The solution algorithm for the 3D inversion problem is also discussed, including the estimation of IC and growth coefficients in 3D. In Section 3, we analyze the ill‐posedness of the 1D problem and present inversion results using our computed regularization. The inversion in 3D is also analyzed for two cases: (1) estimating tumor growth parameters given IC and brain anatomy, and (2) estimating both the brain anatomy and IC while inverting for the tumor growth parameters. Results of the full inversion algorithm are presented on two clinical cases, with a comparison of the reconstructed observed species.

Methodology

2

In the single‐species model, we use the aggregate variable cx,t to represent all the tumor cells. Here x is a point in ℝd (d=3 for 3D model and d=1 for 1D model) and t∈0,T is the time. We denote T as the time‐horizon (the time since the tumor growth onset) and we set T=1 to non‐dimensionalize PDE forward model for the single‐snapshot inversion. In the multi‐species model, we have three species: necrotic nx,t, proliferative px,t, infiltrative ix,t. We can relate the single‐species model to the multi‐species model by setting cx,t=nx,t+px,t+ix,t. We consider probability maps related to MRI scans instead of precise tissue types in our model blue where each voxel is assigned a probability value between 0 and 1 rather than a binary tissue classification. The brain is described by separate maps for gray matter (gx,t), white matter (wx,t), and cerebrospinal fluid (CSF)‐ventricles (VT), denoted by fx,t. We use diffeomorphic image registration to segment the healthy brain into gray matter, white matter, CSF and VT [13].

Formulation (Forward Problem)

2.1

Our model is a go‐or‐grow multi‐species tumor model that accounts for tumor vascularization. Tumor cells are assumed to exist in one of two states: proliferative or invasive. In an environment with sufficient oxygen, tumor cells undergo rapid mitosis by consuming oxygen. However, under low oxygen concentrations (hypoxia), the cells switch from proliferative to invasive. The invasive cells migrate to richer, higher‐oxygen areas and then switch back to the proliferative state. When the oxygen levels drop below a certain level, the tumor cells become necrotic (Figure 1). Necrotic cells are dead cells that remain in this state when oxygen levels stay below the hypoxia threshold. We show the schematic progression of the multi‐species process in Figure model (presented in 1D). Initially, only proliferative cells are present, with the oxygen concentration at its maximum initial value of one. Subsequently, the proliferative cells consume the oxygen, causing them to transition into infiltrative cells. These infiltrative cells then migrate in search of additional oxygen. Our model is summarized as.

$$[eqn]$$ $$[eqn]$$ $$[eqn]$$ $$[eqn]$$ $$[eqn]$$ $$[eqn]$$ $$[eqn]$$ $$[eqn]$$ $$[eqn]$$ $$[eqn]$$ $$[eqn]$$ $$[eqn]$$ $$[eqn]$$

A schematic representation of the multi‐species model dynamics in 1D. The tumor species growth process is represented, with vascularization indicated by oxygen.

The paper uses notations listed in Table 1. The model ignores mass effects from proliferative and necrotic cells and uses conservation equations with transitional (γ,β and α), source (δc, δs, and ℛ), and diffusion (D) terms. See [1, 11] for model choice and description. The reaction operator (ℛ) for species s is modeled as an inhomogeneous nonlinear logistic growth operator as,

$$[eqn]$$

where om and oi are the mitosis and invasive oxygen thresholds, respectively. We restrict the growth conditions with different oxygen levels. We use om=oi+oh2 from [1, 11]. ρm is an inhomogeneous spatial growth term defined as,

$$[eqn]$$

where ρ and ρg are scalar coefficient quantifying proliferation in white matter and gray matter, respectively. We assume 20% growth in gray matter (ρg=0.2ρ) compared to white matter [18, 23, 24]. The growth terms converge to zero when the total tumor concentration (c) reaches 1. The inhomogeneous isotropic diffusion operator D is defined as follows.

$$[eqn]$$

where k defines the inhomogeneous diffusion rate in gray and white matter and we define it as,

$$[eqn]$$

where κ and κg are scalar coefficient quantifying diffusion in white matter and gray matter, respectively. We set κg=0.2κ [18, 23, 24]. Note that the migration and proliferation occur only into white/gray matter and not in CSF/VT. The parameters α,β, and γ functions between are set as follows,

$$[eqn]$$ $$[eqn]$$ $$[eqn]$$

where the scalar coefficients controlling species transitioning are α0,β0, and γ0, and ℋ is the smoothed Heaviside function. There are no movement equation for ventricles since our model does not include mass deformation due to tumor growth. Note that our model differs from the one in [1] in the following ways:

• The transition term αp is changed to αp1−i to prevent transitioning once i≈1 (See Equations (1a) and (1b)).
• The transition term βi is changed to βi1−p to prevent transitioning once p≈1 (See Equations (1a) and (1b)).
• The necrosis term γi and γp are changed to γi1−n and γp1−n to prevent death once n≈1 (See Equations ((1a), (1b), (1c))).
• A term in reaction operator ℛ on a species s is changed from 1−s to 1−c, where c is the total tumor concentration c=i+p+n. This change ensures that there will be no reaction or increase of tumor concentration once the total tumor concentration c≈1 (See Equations (1a), (1b), (1e), and (1f)).
• In Equation (1c), the term γg+w has been removed to ensure that only infiltrative and proliferative cells can transition into the necrotic state.
• We modify the β function to be the complement of the α function (See Equation (7)) in oxygen terms.
• We have incorporated the edema model into the observation operator instead of using a separate model (See Equation (12)).

Observation Operators

2.2

We now discuss how we relate problem predictions to MRI images. Currently, MRI cannot fully reveal tumor infiltration. Instead, we use a relatively standard preprocessing workflow [25] in which the MRI is segmented to different tissue types. For each patient, the preprocessing workflow includes affine registration of the four mpMRI scans to a template atlas, followed by tumor region segmentation using a neural network trained on the BraTS dataset with expert radiologist labels. We define the tumor region and segment healthy substructures (GM, WM, VT, CSF) using ANTs [26] with ensemble atlas‐based segmentation, with these segmentation labels serving as proxies for tissue concentration data. To reconcile the species concentrations generated by our multi‐species model with the observed species segmentation obtained from MRI scans, we need to introduce additional modeling assumptions related to the so‐called observation operator. These operators map the model‐predicted fields to a binary (per tissue type) segmentation, which can then be compared to the MRI segmentation. In particular, we define observation operators for the total tumor region (c) and each individual species, allowing us to determine the most probable species at each point x. The observation operators are defined as follows:

$$[eqn]$$ $$[eqn]$$ $$[eqn]$$ $$[eqn]$$

where Oc,Op,On, and Ol are the observed tumor (c), proliferative (p), necrotic (n) and edema cells (l), respectively. Note that the species used in the above setup are the species at t=1. In this setup, we first determine the most probable location for the tumor (c) using Equation (9), and then we determine the most probable species within Oc using Equations (10) and (11). To model edema, we adopt a simplified approach similar to [1]. In this model, locations with infiltrative concentration above a threshold (ie) are considered as edema if they are not detected as necrotic or proliferative. This model treats edema as a label that is independent of the other PDEs, and can thus be computed in the observation operator Ol using Equation (12). To model the Heaviside function ℋ, we use a sigmoid approximation function given by,

$$[eqn]$$

where ω is the shape factor, determining the smoothness of the approximation. This parameter is set based on the spatial discretization. Thus, the proposed observation operators do not have any free parameters. Next, we discuss the inverse problem to invert for the model parameters.

Inverse Problem

2.3

Our model requires three parameters to characterize tumor progression:

• Healthy precancerous brain segmentation (Ω).
• IC of proliferative tumor cells (p0).
• Growth coefficients vector denoted by q consisting of the following terms: diffusion coefficient κ in Equation (4), reaction coefficient ρ in Equation (2), transition coefficients α0, β0, and γ0 in Equations ((6), (7), (8)), oxygen consumption rate δc in Equation (1d), oxygen supply rate δs in Equation (1d), hypoxic oxygen threshold oh in Equation (2) and (8), invasive oxygen threshold in Equations (2), (6), and (7) and threshold coefficient for edema ie in Equation (12).

We discuss in Sections 2.5.1 and 2.5.2 how we estimate Ω and p0. Here, we focus on growth coefficients q inversion. Given the Ω and p0, we frame the inverse problem as a constrained optimization problem,

$$[eqn]$$ $$[eqn]$$

We aim to minimize the objective function J by optimizing the scalar coefficients vector q in the forward model ℱp0q. This inverse problem is a generalization of the single‐species reaction–diffusion PDE model that has already been shown to be exponentially ill‐posed [2, 27, 28, 29]. We circumvent the ill‐posedness by designing a problem regularization that involves a combination of sparsity constraints, a max constraint in the initial condition and ℓ2 penalty on the unknown parameters, κ and ρ. As we have more undetermined parameters, we cannot use the regularization scheme in [2, 3]. To address this ill‐posedness, we first perform analysis on the Hessian of the optimization problem and then generalize it to the full 3D inversion model. Next, we perform an analysis for the 1D inverse problem and compute a regularization term to assist in estimation of growth coefficients q.

Coefficients Inversion Analysis for a 1D Model

2.4

To demonstrate the ill‐posedness of our model, we analyze a multi‐species model in 1D. To simplify the problem, we assume that the only normal tissue type is white matter (removing the gray matter, ventricles (VT), and CSF from the model). In this setting, the forward model with periodic boundary conditions is described by:

$$[eqn]$$ $$[eqn]$$ $$[eqn]$$ $$[eqn]$$ $$[eqn]$$ $$[eqn]$$ $$[eqn]$$ $$[eqn]$$ $$[eqn]$$ $$[eqn]$$

The optimization problem has exactly the same form as Equations (14) and (15),

$$[eqn]$$ $$[eqn]$$

where pd,nd, and ld are the proliferative, necrotic, and edema data. The corresponding Lagrangian of this problem is given by:

$$[eqn]$$

where ap,ai,an,ao, and aw are the adjoint variables with respect to p,i,n,o, and w.

The first order optimality conditions can be derived by requiring stationary of the Lagrangian with respect to the state variables (p,i,n,o, and w). Therefore, we obtain,

$$[eqn]$$ $$[eqn]$$ $$[eqn]$$ $$[eqn]$$ $$[eqn]$$

Therefore, the adjoint equations are

$$[eqn]$$ $$[eqn]$$ $$[eqn]$$ $$[eqn]$$ $$[eqn]$$ $$[eqn]$$ $$[eqn]$$ $$[eqn]$$ $$[eqn]$$ $$[eqn]$$

And finally, the inversion (gradient) equations for the parameters can be evaluated as,

$$[eqn]$$ $$[eqn]$$ $$[eqn]$$ $$[eqn]$$ $$[eqn]$$ $$[eqn]$$ $$[eqn]$$ $$[eqn]$$ $$[eqn]$$ $$[eqn]$$

where ℋ′ is the derivative of the Heaviside function.

To compute the gradient, it is necessary to first solve the forward problem, as outlined in Equations ((16a), (16b), (16c), (16d), (16e), (16f), (16g), (16h), (16i), (16j)), then solve the adjoint equations, as per Equations ((21a), (21b), (21c), (21d), (21e), (21f), (21g), (21h), (21i), (21j)), and then compute the gradient using Equations ((22a), (22b), (22c), (22d), (22e), (22f), (22g), (22h), (22i), (22j)). To empirically study the ill‐posedness of the 1D inversion problem, we compute the Hessian of the optimization setup through central differencing of the gradients computed in Equations ((22a), (22b), (22c), (22d), (22e), (22f), (22g), (22h), (22i), (22j)).

In Section 3.3.1, the results of the ill‐posedness analysis are presented and the spectrum of eigenvalues of the Hessian for samples (q) from coefficient space is illustrated in Figure 2a. To overcome this challenge, a regularization method is proposed in Algorithm 1 to penalize these ill‐posed directions in the objective function. The method begins by uniformly sampling from the growth coefficient space and estimating the Hessian at the minimum objective function value for each sample. The samples with smaller eigenvalues are deemed to contain less valuable information and, thus, their corresponding eigenvectors are weighted accordingly. A singular value decomposition is then performed to evaluate the weight of ill‐posed directions among different coefficient combinations. The singular values are depicted in Figure 2b. The right singular directions are identified and weighted based on the inverse of the singular values (UR) to mitigate directions with small singular values in the objective function. These directions with smaller singular values are more ill‐posed and, therefore, require more penalization in the objective function.

Ill‐posedness analysis of multi‐species model in 1D. The spectrum of Hessian eigenvalues for 3500 parameter combinations is presented in Figure 2a. The singular values calculated in Algorithm 1 are shown in Figure 2b. We depict the absolute value of last four directions of the computed regularization term (UR) in Figure 2c corresponding to the four smallest singular values. The exponential decrease of eigenvalues in Figure 2a indicates high ill‐posedness of the problem. Additionally, the wide range of small eigenvalues suggests varying degrees of ill‐posedness among samples, which is accounted for in Algorithm 1. Figure 2b demonstrates how different directions impact the ill‐posedness, with smaller singular values indicating greater ill‐posedness. Thus, using the inverse of the singular directions can lead to improved inversion results. The results in Figure 2c clearly indicates ill‐posedness for high values γ0 and a combination of high values of ρ and high values of δc contributes to ill‐posedness.

Specifically, we uniformly sample N=3500 parameter values and compute the Hessian matrix Hq at each objective function minimum. We then perform singular value decomposition on these Hessian matrices (H=UΣVT) to identify the principal directions of ill‐posedness. From this decomposition, we construct the regularization operator URTq, which we incorporate into the regularized objective function:

$$[eqn]$$

where λ controls the regularization strength. The fast decay of the Hessian eigenvalues shown in Figure 2a indicates high ill‐posedness of the problem. Additionally, the wide range of small eigenvalues suggests varying degrees of ill‐posedness among samples, which is accounted for in our regularization approach. Figure 2b demonstrates how different directions impact the ill‐posedness, with smaller singular values indicating greater ill‐posedness.

ALGORITHM 1Algorithm to compute the regularization.1: procedure (n, bl, bu)         ⊳ Get the number of samples n, lower bound vector bl and upper bounds vector bu 2:     qii=1n∼Ublbu                          ⊳ Take uniform samples of parameters3:     D←0 4:     for i=1…n do 5:         H←Hessqi                                    ⊳ Compute the Hessian6:        Λ,Q←EvalDecompH                    ⊳ Compute Evals. and Evecs. of Hessian7:        D←D∪Λ12Q                         ⊳ Store the weighted ill‐posed directions8:     end for 9:     Σ,U←SVDD                                 ⊳ SVD on ill‐posed directions10:    UR←UΣ−1                                ⊳ Weighting the ill‐posed directions11: end procedure

3D Inversion

2.5

As mentioned, the inverse problem requires three parameters to characterize the tumor progression: (a) Healthy brain anatomy Ω to describe wx,0,gx,0,fx,0. (b) Initial condition of proliferative tumor cells p0 and (c) growth coefficients vector q. In the following, we explain how we estimate each of these parameters.

Brain Anatomy

2.5.1

The process of deriving the healthy brain anatomy (Ω) involves multiple steps. First, an MRI scan is affine registered to a template atlas. Then, a neural network [18] trained on the BraTS dataset [25] is used to segment the tumor into proliferative, necrotic, and edema regions. Next, 10 normal brain scans with known healthy segmentation are registered to the patient's scan using [26]. An ensemble‐based approach is then used to combine the atlas‐based segmentations for the different healthy tissue labels [5, 13]. We combine the healthy tissue labels and tumor segmentation as a single brain segmentation of the patient. Finally, the tumorous regions are replaced with white matter to estimate the healthy brain anatomy Ω.

IC Inversion

2.5.2

To estimate IC of proliferative tumor cells, we use the single‐species IC inversion proposed in [3]. This algorithm uses an adjoint based method for reaction–diffusion PDE model constrained with sparsity and max constraints. Our current inverse IC scheme deviates from the [3] algorithm by limiting the support region to the necrotic areas only, as our simulations in [1] revealed that the IC always present in the necrotic regions.

Coefficients Inversion

2.5.3

To invert for the growth coefficients q in 3D, we use the regularization term computed from the 1D model as sampling the Hessian is very expensive in 3D. As the structure of the nonlinear operators is the same, our hypothesis is the 1D regularization operator can be used to stabilize the 3D reconstruction. One contribution of this paper is to verify the hypothesis using numerical experiments with synthetically constructed data. Given the IC (p0) and brain anatomy (Ω) estimated from Sections 2.5.1 and 2.5.2, we can solve the following optimization problem with constraints:

$$[eqn]$$ $$[eqn]$$

Results

3

In this study, we evaluate our algorithm on synthetic cases and clinical dataset. We carry out the 3D inversions using the Frontera and Lonestar system at the Texas Advanced Computing Center (TACC) at The University of Texas at Austin. Our solver is written in both C++ and Python. The forward solver uses C++ with PETSc, AccFFT, PnetCDF. The inverse problem is solved using the covariance matrix adaptation evolution strategy (CMA‐ES) in Python [30, 31]. Each 3D inversion takes 32 h on a single NVIDIA Quadro RTX 5000 GPU. For 1D test‐cases, we use Python for the forward problem with the same CMA‐ES optimizer.

Numerical Parameters

3.1

We list all numerical parameters used in our solver in Table 2. We now explain a justification for our choice:

• For the multi‐species forward solver, we use the same setting and numerics of [3].
• For IC inversion, we use the solver from our previous studies [2, 3]. The single‐species reaction diffusion model with ℓ0 constraint is used to solve a single‐snapshot inverse problem iteratively. The IC is estimated within necrotic regions and the sparsity parameter is fixed at 5 to allow up to 5 voxels to be considered for IC. The experiments on choosing this value can be found in [3].
• We choose 1603 as spatial discretization for 3D inversion to balance computational cost with detail preservation, reducing computation from standard 2563 while maintaining tumor characteristics and volume measurements.
• We determine the optimal value for the regularization parameter λ by evaluating the quality of data reconstruction with 20% additive noise. We select this noise percentage because previous studies have demonstrated that the tumor segmentation tool [25] achieves an accuracy of approximately 80%.
• For the inversion of growth coefficients q, we employ CMA‐ES with a relative function tolerance of 1×10−8 and an absolute function tolerance of 1×10−3. Our analysis indicates that reducing these tolerance values further does not result in improved outcomes.
• The sample size per iteration for the CMA‐ES is 16 as a balance between evaluation costs and having a diverse population for covariance sampling.
• We set the initial guess for the CMA‐ES solver to the midpoint of the parameter bounds and the initial variance to half the parameter range. Despite experimenting with higher variances and different initial guesses, we did not observe any improvement in solver performance.
• We use a shape factor ω of 32 for all the Heaviside functions except for ℋi−ie to achieve a smooth transition. For the edema region, a smooth Heaviside function in ℋi−ie can lead to high spatial diffusivity of infiltrative cells i. To prevent this, we use a larger shape factor (ω=256) for the edema region compared to the others. We did not observe sensitivity to this parameter in our results.

Performance Measures

3.2

Our main goal of the numerical experiments is to show how the proposed regularization can stabilize the multi‐species inversion and then use the same regularization to perform a full clinical inversion. First, we show how the regularization helps the inversion in synthetic cases and then combine it with IC inversion on clinical scans. In this regard, we define the following metrics to quantify our reconstructions,

• relative error in growth coefficient ι:

$$[eqn]$$

where ι* is ground truth coefficient and ιrec is the reconstructed value. We denote the growth coefficients vector with q and eq denotes the average relative error of growth coefficients.

• relative error in final species reconstruction st=1 where s can be n, p or i:

$$[eqn]$$

where srec1 is the final reconstructed species, s1 is the final species concentration grown with ground truth growth coefficients q.

• relative error in initial tumor concentration:

$$[eqn]$$

where prec0 is the reconstructed IC and p*0 is the ground truth IC.

• relative ℓ1 error in initial tumor parameterization (p0):

$$[eqn]$$

where p0rec is the reconstructed tumor IC and p0* is the ground truth IC.

• relative ℓ1 reconstruction for projected IC along a axis:

$$[eqn]$$

where a represents the axis which can be x,y or z.

• reconstruction Dice score for an observed species s:

$$[eqn]$$

where ∣srec∣ and ∣s*∣ are the cardinality number of the observed species from inversion and data, respectively and ∣srec∩s*∣ is the cardinality number of the intersection between inverted species srec and species from data s*.

• reconstruction Dice score for tumor core (Op∪On) denoted as tc

$$[eqn]$$

where ∣tcrec∣ and ∣tc*∣ are the cardinality number of tumor core from inversion and data, respectively and ∣tcrec∩tc*∣ is the cardinality number of the intersection between inverted tumor core tcrec and tumor core from data tc*. Note that we denote the Dice score from single‐species model with πtcss and Dice score from multi‐species model with πtcms.

1D Test‐Case

3.3

In this section, first, we discuss the ill‐posedness of the problem in 1D and then, we present the inversion results. We present the analysis of the 1D inversion assuming that the IC is known and we only need to invert for the growth coefficients.

III‐Posedness Analysis

3.3.1

To address the ill‐posedness of the problem, we evaluate the Hessian for 3500 samples of q uniformly derived from coefficients range and perform an eigenvalue decomposition. The eigenvalue spectrum, displayed in Figure 2a, exhibits an exponential decrease, indicating that the problem is highly ill‐posed. The wide range of eigenvalues suggests that different combinations of coefficients contribute to the problem's ill‐posedness with varying degrees of difficulty to reconstruct the coefficients.

Following Algorithm 1, we further analyze the ill‐posed directions by performing singular value decomposition (SVD) on weighted eigenvectors. The result of the SVD is displayed in Figure 2b through the plot of the singular values. The smaller singular value suggests higher ill‐posedness of the corresponding coefficient combinations. In Figure 2c, we have illustrated the last four directions of the computed UR, which correspond to the 4 smallest singular values from Figure 2b. It is crucial to consider the weight for penalizing the γ0 parameter in the first direction. This penalization corresponds to the necrotic species and accounts for the fact that if a location is assigned to the necrotic region in both the data and observation, increasing the γ0 coefficient will not reduce the error as the objective function is insensitive to it. Furthermore, we have observed through our simulations that high values of δc and ρ in the second and third directions contribute to the ill‐posedness.

Inversion Analysis in 1D With Noise

3.3.2

In this section, we present results from inverting 1D synthetic cases, where exact IC is known and focus is on finding growth coefficients q. To test inversion stability, we generate species concentrations (s*) using ground truth growth coefficients q*, and we generate the noisy species s^=s1+noise where noise has a zero‐mean Gaussian distribution with variance σ (N0σ). We then apply the observation operators to generate the noisy observed data. We vary σ to generate different average noise levels among all species p,n and i. This noisy data is then input for the inversion algorithm and tested with and without regularization to evaluate stability and performance of reconstruction.

In Figure 3, we show reconstructed species concentrations (proliferative p, infiltrative i, necrotic n) for synthetic cases with different levels of independent Gaussian noise added to each signal. The results with and without regularization are displayed. For a quantitative comparison, see Table 3 where regularization generally improves reconstruction performance compared to non‐regularized results. In Table 3, we also observe that regularization avoids parameter inversion problems, as seen in the case with 10% noise where some parameters diverge by comparing the average relative errors for parameters reconstruction. Table A1 reports inverted coefficients for each parameter with and without regularization. Overall, the inversion scheme with regularization is more stable in terms of species reconstruction and parameter inversion.

1D inversion results of synthetic data with additive Gaussian noise (5%,10%,20%,30%). Here we show the reconstructed species concentration at T=1. First to fourth rows show relative noise levels. Plots in first to third columns show proliferative, necrotic, and infiltrative species, respectively. The green, gray, magenta, and blue lines represent the underlying ground truth species (before observation denoted as “Ground Truth”), ground truth species with noise (“Noisy Ground Truth”), reconstructed species without regularization (“Non‐Reg. Inv.”), and reconstructed species with regularization (“Reg. Inv.”), respectively. Regularized inversion is observed to better match the underlying species even when the data is perturbed by noise, particularly in the case of infiltrative cells where deviations are greatly reduced with regularization.

Further, we examine the influence of varying coefficient combinations on the performance of the inversion process. To achieve this objective, we conduct experiments on 150 synthetic data sets contaminated with 20% Gaussian noise and the results are presented in Table A2. The reconstruction measures are calculated both for regularized and non‐regularized inversion cases. Our results indicate that the use of regularization leads to a noticeable reduction in the average relative error of the coefficients (eq), with an improvement of approximately 30%. This improvement is particularly noticeable for the coefficients γ0,β0, and α0.

3D Test‐Cases

3.4

In this section, we analyze the 3D inversion. Results for cases with known brain anatomy (Ω) and IC (p0) is presented first, where our goal is to estimate growth coefficients q only. The next case involves inverting for both IC and growth coefficients. Finally, results for real clinical data from BraTS20 training dataset [25] are presented.

Artificial Tumor With Known IC and Anatomy

3.4.1

In this section, we evaluate inversion for growth parameters only. Given a ground truth q*, we solve the forward problem to generate data. We test different noise levels of 5%,10%,20%, and 30% and present the reconstruction with and without regularization in Figure 4. We report the results in Table 4 and individual model parameters and their relative errors in Table A3.

Evaluating 3D inversion using synthetic data with known IC and brain anatomy. We start by generating a tumor test‐case and showing tumor species including proliferative or enhancing (EN), necrotic (NEC) and edema (ED) and segmentation in the first row (a). We then add 5%, 10%, 20%, and 30% Gaussian noise to the species and apply an observation operator to get noisy data. We perform inversion with and without regularization on this data. The inversion results are displayed in rows 2–4, with the enhancing, necrotic, and infiltrative concentrations in columns 1–3, and the segmentation in the last column. The regularized inversion yields acceptable segmentation for all noise levels, whereas the non‐regularized case shows more variability for different noise levels.

As seen in Table 4 and Figure 4, regularization noticeably improves the reconstruction error eq in 5% noise, while improving the reconstruction at higher noise levels. Although Dice scores and segmentation do not show substantial improvement, the underlying species reconstruction, not directly observed from data, is improved.

Artificial Tumor With Unknown IC and Anatomy

3.4.2

In this section, we present the results for two‐stage inversion process for both IC (p0) and growth coefficients (q). Our solution involves first the tumorous regions are filled with white matter and the healthy brain is estimated and we invert for the IC using a single‐species model, as detailed in [2, 3]. Once the IC is estimated, we estimate the growth coefficients.

Three cases are tested with generated data and different parameter combinations. The inversion results are depicted in Figure 5 and the corresponding numerical measures are reported in Table 5. The individual growth coefficients and their relative error are also included in Table A4.

3D inversion result with synthetic test‐cases generated with an unknown IC and brain anatomy. Different normal brains with no mass deformation are used to generate the synthetic cases. In each row, we present a single test‐case, with tumor species (proliferative or enhancing (EN), necrotic (NEC), edema (ED) and infiltrative cells) from the synthetic cases and inversion displayed in the first, second, and third columns, respectively. The given data and segmentation computed by the inversion algorithm are depicted in the fourth column. Moreover, we depicted the projected ground truth IC and estimated IC in the segmentation for data and inversion, respectively, Although the IC is an estimate, the reconstruction from the inversion is good as evidenced by the qualitative observation. The challenge of reconstructing multi‐species arises due to the unknown initial brain anatomy and isotropic migration of tumor in all directions due to assigned white matter label.

As we can observe, we demonstrate that even though the IC estimation is derived from another model, it is still capable of reconstructing the observed tumor species from the data. Our test cases, as shown in Figure 5, indicate that the IC estimation provides a relatively good estimate for the initial condition. The estimation of the IC for a reaction–diffusion tumor growth model is exponentially ill‐posed, making it a challenging task, as noted in prior studies such as [2, 3].

As shown in Table 5, our solver not only achieves a reconstruction comparable to that of the single‐species model but also provides additional information on the underlying species. Our method enables the characterization of observed species and the estimation of non‐observable ones from the given data. As demonstrated in Section 3.4.1, our algorithm exhibits good convergence once a suitable estimate of IC is obtained. Notably, the performance of our method is influenced by the brain tumor anatomy, with the diffusion patterns substantially impacted by the tissue types in the test cases. Our model simplifies the problem by estimating the tumor growth model, given the challenge of accurately determining brain anatomy due to mass deformation. Furthermore, as illustrated in Figure 5, the diffusion of proliferative cells is primarily in all directions due to the allocation of white matter in the brain anatomy, while the data is influenced by the gray matter.

Inversion on Two Clinical Data

3.5

In this section, we demonstrate the application of our inversion scheme to two patients from the BraTS20 dataset [25]. The procedure for reconstruction is outlined in Section 3.4.2 and involves first reconstructing IC of the tumor, followed by estimation of the growth parameters. The reconstructed species are visualized in Figure 6, with the projected IC estimate shown in the inverted segmentation. We report that the results can be found in Table 6 and Table A5. Note that the only information available from the patients is the segmentation of the observed species from MRI scans, and direct comparison of the underlying species concentrations is not possible.

We conduct a multi‐species inversion of the BraTS20_Training clinical dataset, selecting two patients, using the MRI scan segmentation as input. Our results are presented in a tabular format, with each row showing the inversion and data for each patient. The first three columns depict the inverted species (proliferative or enhancing (EN), necrotic (NEC), and infiltrative cells), the fourth column shows the observed species from the inversion, and the last column displays the data. Our results demonstrate that, despite only having a single snapshot from each patient, we can qualitatively reconstruct the species distribution with reasonable accuracy.

Our solver achieves a reconstruction accuracy comparable to that of the single‐species model. Similar to the findings in Section 3.4.2. Furthermore, we are able to reconstruct species data and quantify infiltrative tumor cells. However, the reconstruction may be affected by IC estimation using the single‐species model, particularly for BraTS20_Training_039, as evidenced by the presence of a considerable amount of necrotic tissue in the initial condition location. Furthermore, the anisotropic progression of edema in BraTS20_Training_039 presents a challenging task for the solver to identify infiltrative cells. A similar trend can be observed for BraTS20_Training_042, where the edema does not uniformly diffuse. Our edema model is a simplified representation used in [1, 11]. Additionally, the results may be influenced by the species' volume in the data, as our unweighted objective function prioritizes reducing errors for the dominant species, that is, necrotic tissue in BraTS20_Training_042.

Conclusion

4

We presented an inverse solver for tumor growth model from one time snapshot and a regularization scheme to allow stable reconstruction of model parameters. The solver estimates the initial condition of the multi‐species model and growth coefficients. The input data is the MRI‐scan based tumor segmentation. Our solver uses just a single multiparametric scan, and we design a single‐snapshot inverse solver for the multi‐species model.

Our results demonstrate the efficacy of this regularization in improving the estimation of the growth coefficients in one‐dimensional and three‐dimensional settings, even in the presence of high relative noise. Moreover, we combine the regularization terms with our prior work on tumor initial condition estimation and apply the inversion algorithm to both synthetic and clinical cases, revealing the reconstruction's stability for the underlying species concentration. Our solver not only estimates the tumor core in a range comparable to that of the single‐species model, but also provides additional information on the reconstruction of other species. In this way, we can estimate the non‐observable species for the tumor growth model.

The solver is implemented using a scalable sampling optimization scheme, demonstrating a reasonable computational time. For initial condition estimation, Bayesian inference on the multispecies model with prior regularization can be used. However, such approaches may require many evaluations of the forward problem and its adjoint and are computationally prohibitive. Although the estimation of the initial condition could be improved (e.g., by using Bayesian inference with prior regularization, despite its high computational costs) to enhance the growth coefficients and reconstructions, our results provide valuable insights into the inverse problem for multi‐species tumor growth models.

In the future, we plan to conduct a comprehensive evaluation of the inversion scheme on a large number of clinical images to observe the correlation between the growth coefficients and clinical data outcomes, such as survival rate. As the model includes more detailed biophysical processes than the single‐species model, there are numerous opportunities for future extensions, including the inclusion of anisotropic diffusion and mass deformation for a more accurate modeling of GBM growth.

Ethics Statement

The authors have nothing to report.

Conflicts of Interest

The authors declare no conflicts of interest.