Image-Driven Biophysical Tumor Growth Model Calibration

TL;DR

This paper introduces a novel MRI-based method for calibrating biophysical tumor growth models from a single snapshot by using an atlas as a pretumor reference, improving accuracy through joint registration and parameter estimation.

Contribution

The paper presents a new scheme that executes tumor growth simulation in the patient brain domain and combines image registration with parameter estimation, incorporating sparsity constraints within a PDE-constrained optimization framework.

Findings

More accurate tumor parameter reconstruction compared to previous methods

Effective integration of image registration and parameter estimation

Validated on real brain data with synthetic tumors

Abstract

We present a novel formulation for the calibration of a biophysical tumor growth model from a single-time snapshot, MRI scan of a glioblastoma patient. Tumor growth models are typically nonlinear parabolic partial differential equations (PDEs). Thus, we have to generate a second snapshot to be able to extract significant information from a single patient snapshot. We create this two-snapshot scenario as follows. We use an atlas (an average of several scans of healthy individuals) as a substitute for an earlier, pretumor, MRI scan of the patient. Then, using the patient scan and the atlas, we combine image-registration algorithms and parameter estimation algorithms to achieve a better estimate of the healthy patient scan and the tumor growth parameters that are consistent with the data. Our scheme is based on our recent work (Scheufele et al, "Biophysically constrained diffeomorphic image registration, Tumor growth, Atlas registration, Adjoint-based methods, Parallel algorithms", CMAME, 2018), but apply a different and novel scheme where the tumor growth simulation in contrast to the previous work is executed in the patient brain domain and not in the atlas domain yielding more meaningful patient-specific results. As a basis, we use a PDE-constrained optimization framework. We derive a modified Picard-iteration-type solution strategy in which we alternate between registration and tumor parameter estimation in a new way. In addition, we consider an $\ell_1$ sparsity constraint on the initial condition for the tumor and integrate it with the new joint inversion scheme. We solve the subproblems with a reduced-space, inexact Gauss-Newton-Krylov/quasi-Newton methods. We present results using real brain data with synthetic tumor data that show that the new scheme reconstructs the tumor parameters in a more accurate and reliable way compared to our earlier scheme.

Tables

Table 1. Table 1: Notation for the Moving Atlas joint tumor inversion and image registration formulation: 𝒎 𝒎 \boldsymbol{m} denotes the vector of anatomy labels (probability maps for white matter (WM), gray matter (GM) and cerebrospinal fluid (CSF)) defining the brain geometry, and c 𝑐 c denotes a probability map for tumor (TU). Subscripts A 𝐴 A and P 𝑃 P indicate variables in the atlas and patient space, respectively. Anatomy labels and tumor label with subscript D 𝐷 D denote the patient input imaging data (after prepossessing). The fields evolve in time and space; anatomy labels 𝒎 𝒎 \boldsymbol{m} for brain tissue evolve along a pseudo-time t R ∈ [ 0 , 1 ] subscript 𝑡 𝑅 0 1 t_{R}\in[0,1] associated to registration, and tumor labels evolve along the (normalized) time t T ∈ [ 0 , 1 ] subscript 𝑡 𝑇 0 1 t_{T}\in[0,1] of tumor growth. The integration of a simulated tumor map c ​ ( 1 ) 𝑐 1 c(1) into healthy brain tissue 𝒎 𝒎 \boldsymbol{m} is indicated by 𝒎 ′ superscript 𝒎 ′ \boldsymbol{m}^{\prime} and modeled via the formula 𝒎 ′ = 𝒎 ​ ( 1 − c ​ ( 1 ) ) superscript 𝒎 ′ 𝒎 1 𝑐 1 \boldsymbol{m}^{\prime}=\boldsymbol{m}(1-c(1)) . To judge the proximity of the predicted state of our mathematical model to the observed data 𝒎 D subscript 𝒎 𝐷 \boldsymbol{m}_{D} , and c D subscript 𝑐 𝐷 c_{D} , we define the ℓ 2 subscript ℓ 2 \ell_{2} -distance measures 𝒟 c subscript 𝒟 𝑐 \mathcal{D}_{c} and 𝒟 𝒎 subscript 𝒟 𝒎 \mathcal{D}_{\boldsymbol{m}} .

healthy atlas brain (input)	${𝒎}_A$
patient brain with tumor (input; target data)	${𝒎}_D$
patient tumor (input; target data)	$c_D$
healthy patient brain	${𝒎}_P\left(1\right)$
tumor initial condition	$c\left(0\right)$
simulated tumor	$c\left(1\right)$
approximated patient brain	${𝒎}_P^{\prime }\left(1\right)$
${\ell }_2$-misfit tumor	${𝒟}_c[c_1,c_2]:=\frac{1}{2}\parallel c_1-c_2{\parallel }_{L^2\left(\Omega \right)}^2$
${\ell }_2$-misfit anatomy labels	${𝒟}_{𝒎}[{𝒎}_1,{𝒎}_2]:=\frac{1}{2}\parallel {𝒎}_1-{𝒎}_2{\parallel }_{L^2{\left(\Omega \right)}^3}^2$

Table 2. Table 2: Summary of common parameters used in all test cases. We report values of the following parameters: n 𝑛 n denotes the image resolution with n i subscript 𝑛 𝑖 n_{i} , i = 1 , 2 , 3 𝑖 1 2 3 i=1,2,3 ; n 𝒑 subscript 𝑛 𝒑 n_{\boldsymbol{p}} is the number of Gaussian for the parametrization of the tumor initial condition, σ 𝜎 \sigma is the standard deviation of the associated Gaussian basis functions (and varies throughout test cases), δ 𝛿 \delta denotes the spacing in between centers of adjacent Gaussians; k 1 subscript 𝑘 1 k_{1} , k 2 subscript 𝑘 2 k_{2} and k 3 subscript 𝑘 3 k_{3} are the characteristic diffusion parameters for WM, GM, and CSF, k f subscript 𝑘 𝑓 k_{f} the overall scaling parameter for the isotropic part of the diffusion coefficient for net migration of cancerous cells into surrounding tissue; ρ 1 subscript 𝜌 1 \rho_{1} , ρ 2 subscript 𝜌 2 \rho_{2} , and ρ 3 subscript 𝜌 3 \rho_{3} are the characteristic reaction factors for WM, GM, and CSF, ρ f subscript 𝜌 𝑓 \rho_{f} is the overall reaction scaling factor; opttol R subscript opttol 𝑅 \text{opttol}_{R} , opttol T subscript opttol 𝑇 \text{opttol}_{T} are the convergence tolerances for registration and tumor inversion; maxit i = ( maxit i , N , maxit i , K ) subscript maxit 𝑖 subscript maxit 𝑖 𝑁 subscript maxit 𝑖 𝐾 \text{maxit}_{i}=(\text{maxit}_{i,N},\text{maxit}_{i,K}) denotes the maximum number of Newton iterations and Krylov iterations (for the KKT system) for the tumor inversion ( i = T 𝑖 𝑇 i=T ) and registration ( i = R 𝑖 𝑅 i=R ), respectively; β 𝒑 subscript 𝛽 𝒑 \beta_{\boldsymbol{p}} is the regularization parameter for the tumor inversion; β 𝒗 i ​ n ​ i ​ t superscript subscript 𝛽 𝒗 𝑖 𝑛 𝑖 𝑡 \beta_{\boldsymbol{v}}^{init} and β 𝒗 f ​ i ​ n ​ a ​ l superscript subscript 𝛽 𝒗 𝑓 𝑖 𝑛 𝑎 𝑙 \beta_{\boldsymbol{v}}^{final} are the initial and final values for the β 𝛽 \beta -continuation scheme, applied in image registration; ε ∇ subscript 𝜀 ∇ \varepsilon_{\nabla} is the bound on the variation of the deformation gradient det ​ ( ∇ 𝒚 ) det ∇ 𝒚 \mbox{det}(\nabla\boldsymbol{y}) used in the continuation scheme.

Description	Parameter(s)	Value(s)
image resolution	$n=n_1\cdot n_2\cdot n_3$	$n_i=128$
tumor initial condition parametrization	$n_{𝒑}$, $\sigma$, $\delta$	$\{125,343\}$, $\{\frac{2\pi }{30},\frac{2\pi }{64}\}$, $1.5$
tumor growth characteristic cell diffusivity	$\overline{k}=k_f^{\star }{(k_1,k_2,k_3)}^T$	$0.1\cdot {(1,0,0)}^T$
tumor growth characteristic net cell proliferation	$\overline{\rho }={\rho }_f^{\star }{({\rho }_1,{\rho }_2,{\rho }_3)}^T$	$15\cdot {(1,0,0)}^T$
optimizer tolerances	${\text{opttol}}_R$, ${\text{opttol}}_T$	$1\text{E}-3$, $1\text{E}-3$
maximum Newton / Krylov iterations, tumor inversion	${\text{maxit}}_T=({\text{maxit}}_{T,N},{\text{maxit}}_{T,K})$	$(50,-)$
maximum Newton / Krylov iterations, registration	${\text{maxit}}_R=({\text{maxit}}_{R,N},{\text{maxit}}_{R,K})$	$(50,80)$
tumor inversion regularization parameter	${\beta }_{𝒑}$	$1\text{E}-4$
registration regularization parameter (continuation)	${\beta }_{𝒗}^{init}$, ${\beta }_{𝒗}^{final}$	$1$, $1\text{E}-04$
bound on local volume change (variation $\text{det}\left(\nabla 𝐲\right)$)	${\epsilon }_{\nabla }$	$1\text{E}-2$

Table 3. Table 3: STRV test case variations. Parameter choices for the synthetic tumor with real velocity (STRV) test case; ground truth: ( ρ f ⋆ = 15 subscript superscript 𝜌 ⋆ 𝑓 15 \rho^{\star}_{f}=15 , ρ 1 ⋆ = 1 subscript superscript 𝜌 ⋆ 1 1 \rho^{\star}_{1}=1 , ρ 1 ⋆ = 0 subscript superscript 𝜌 ⋆ 1 0 \rho^{\star}_{1}=0 , k f ⋆ = 1.0 ​ E ​ − 1 subscript superscript 𝑘 ⋆ 𝑓 1.0E-1 k^{\star}_{f}=$1.010-1$ , k 1 ⋆ = 1 subscript superscript 𝑘 ⋆ 1 1 k^{\star}_{1}=1 , k 1 ⋆ = 0 subscript superscript 𝑘 ⋆ 1 0 k^{\star}_{1}=0 , 𝒑 = 𝒑 ⋆ 𝒑 superscript 𝒑 ⋆ \boldsymbol{p}=\boldsymbol{p}^{\star} ). For the initial condition parametrization, we use a regular grid of n 𝒑 = 125 subscript 𝑛 𝒑 125 n_{\boldsymbol{p}}=125 (large σ 𝜎 \sigma ), and n 𝒑 = 343 subscript 𝑛 𝒑 343 n_{\boldsymbol{p}}=343 (small σ 𝜎 \sigma ) Gaussian basis functions, respectively with standard deviation σ 𝜎 \sigma as outlined below, and a spacing of δ = 1.5 ​ σ 𝛿 1.5 𝜎 \delta=1.5\,\sigma . The grid is centered around the positions 𝒙 c 1 = 2 ​ π / 128 ​ ( 39 , 63 , 64 ) subscript 𝒙 subscript 𝑐 1 2 𝜋 128 39 63 64 \boldsymbol{x}_{c_{1}}=\nicefrac{{2\pi}}{{128}}(39,63,64) (STRV-C1) and 𝒙 c 2 = 2 ​ π / 128 ​ ( 61 , 89 , 64 ) subscript 𝒙 subscript 𝑐 2 2 𝜋 128 61 89 64 \boldsymbol{x}_{c_{2}}=\nicefrac{{2\pi}}{{128}}(61,89,64) (STRV-C2), respectively (cf. Fig.  4 and Fig.  6 for the location of tumor seed).

test case	setting #1	setting #2	setting #3
(STRV-C1)	$\sigma =2\pi /30$, $T=0.16$	$\sigma =2\pi /64$, $T=0.32$	$\sigma =2\pi /64$, $T=0.44$
(STRV-C2)	$\sigma =2\pi /64$, $T=0.16$	$\sigma =2\pi /64$, $T=0.32$

Table 4. Table 5: Comparison of the Moving Patient , and the new Moving Atlas scheme. We report results for the (STRV-C1) and (STRV-C2) test cases with centers 𝒙 c 1 subscript 𝒙 subscript 𝑐 1 \boldsymbol{x}_{c_{1}} and 𝒙 c 2 subscript 𝒙 subscript 𝑐 2 \boldsymbol{x}_{c_{2}} , and parameter setting # ​ 1 # 1 \#1 from Tab.  3 , respectively. Here we are not inverting for k f subscript 𝑘 𝑓 k_{f} , but only invert for the tumor initial condition and registration velocity. We analyze the sensitivity of the joint inversion with respect to perturbations in the characteristic diffusivity k f subscript 𝑘 𝑓 k_{f} by inverting with different choices of k f subscript 𝑘 𝑓 k_{f} ; bold numbers correspond to inversion with the ground truth value ( k f ⋆ = 1 ​ E ​ − 1 superscript subscript 𝑘 𝑓 ⋆ 1E-1 k_{f}^{\star}=$110-1$ ). We always use the ground truth proliferation rate ρ f = 15 subscript 𝜌 𝑓 15 \rho_{f}=15 . We report the average mismatch for the anatomy labels for the pathologic brain tissue labels μ B , L 2 subscript 𝜇 𝐵 superscript 𝐿 2 \mu_{B,L^{2}} , the healthy brain tissue labels μ B 0 , L 2 subscript 𝜇 subscript 𝐵 0 superscript 𝐿 2 \mu_{B_{0},L^{2}} and the tumor μ T , L 2 subscript 𝜇 𝑇 superscript 𝐿 2 \mu_{T,L^{2}} , the mean Dice coefficient for brain tissue DICE B subscript DICE 𝐵 \text{DICE}_{B} , for the healthy patient brain tissue DICE B 0 subscript DICE subscript 𝐵 0 \text{DICE}_{B_{0}} and for the tumor DICE T subscript DICE 𝑇 \text{DICE}_{T} (in cases without Dice score for the tumor reconstruction, the tumor probability map has values below 0.5 0.5 0.5 everywhere). Furthermore, we measure the ℓ 2 subscript ℓ 2 \ell_{2} -error e c ​ 0 , L 2 subscript 𝑒 𝑐 0 superscript 𝐿 2 e_{c0,L^{2}} of the reconstructed tumor initial condition to the ground truth.

ID			$k_f$	${\mu }_{B,L^2}$	${\text{DICE}}_B$	${\mu }_{B_0,L^2}$	${\text{DICE}}_{B_0}$	${\mu }_{T,L^2}$	${\text{DICE}}_T$	$e_{c0,L^2}$
$\#7$	STRV-C1 (#1)	MP	$5\text{E}-1$	$\mathrm{3.799\hspace{0.17em}000}\text{E}-01$	$\mathrm{8.618\hspace{0.17em}000}\text{E}-01$	$\mathrm{2.578\hspace{0.17em}000}\text{E}-01$	$\mathrm{8.786\hspace{0.17em}333}\text{E}-01$	$\mathrm{3.962\hspace{0.17em}000}\text{E}-01$	–	$\mathrm{4.993\hspace{0.17em}000}\text{E}-01$
$\#8$			$3\text{E}-1$	$\mathrm{3.763\hspace{0.17em}000}\text{E}-01$	$\mathrm{8.644\hspace{0.17em}333}\text{E}-01$	$\mathrm{2.566\hspace{0.17em}000}\text{E}-01$	$\mathrm{8.789\hspace{0.17em}000}\text{E}-01$	$\mathrm{2.878\hspace{0.17em}000}\text{E}-01$	$\mathrm{7.757\hspace{0.17em}000}\text{E}-01$	$\mathrm{8.684\hspace{0.17em}000}\text{E}-01$
$\#9$			$k_f^{\star }=1\text{E}-1$	$\mathrm{3.716\hspace{0.17em}000}\text{E}-01$	$\mathrm{8.721\hspace{0.17em}333}\text{E}-01$	$\mathrm{2.555\hspace{0.17em}000}\text{E}-01$	$\mathrm{8.787\hspace{0.17em}667}\text{E}-01$	$\mathrm{1.454\hspace{0.17em}000}\text{E}-01$	$\mathrm{9.209\hspace{0.17em}000}\text{E}-01$	$\mathrm{3.225\hspace{0.17em}000}\text{E}-01$
$\#10$			$1\text{E}-3$	$\mathrm{3.727\hspace{0.17em}000}\text{E}-01$	$\mathrm{8.749\hspace{0.17em}333}\text{E}-01$	$\mathrm{2.594\hspace{0.17em}000}\text{E}-01$	$\mathrm{8.729\hspace{0.17em}000}\text{E}-01$	$\mathrm{2.491\hspace{0.17em}000}\text{E}-01$	$\mathrm{8.746\hspace{0.17em}000}\text{E}-01$	$\mathrm{5.392\hspace{0.17em}000}\text{E}-01$
$\#11$		MA	$5\text{E}-1$	$\mathrm{2.396\hspace{0.17em}000}\text{E}-01$	$\mathrm{9.455\hspace{0.17em}000}\text{E}-01$	$\mathrm{1.617\hspace{0.17em}000}\text{E}-01$	$\mathrm{9.468\hspace{0.17em}667}\text{E}-01$	$\mathrm{5.645\hspace{0.17em}000}\text{E}-01$	–	$\mathrm{1.372\hspace{0.17em}000}\text{E}00$
$\#12$			$3\text{E}-1$	$\mathrm{2.400\hspace{0.17em}000}\text{E}-01$	$\mathrm{9.371\hspace{0.17em}667}\text{E}-01$	$\mathrm{1.632\hspace{0.17em}000}\text{E}-01$	$\mathrm{9.477\hspace{0.17em}667}\text{E}-01$	$\mathrm{2.863\hspace{0.17em}000}\text{E}-01$	$\mathrm{3.658\hspace{0.17em}000}\text{E}-01$	$\mathrm{9.412\hspace{0.17em}000}\text{E}-01$
$\#13$			$k_f^{\star }=1\text{E}-1$	$\mathrm{2.345\hspace{0.17em}000}\text{E}-01$	$\mathrm{9.427\hspace{0.17em}667}\text{E}-01$	$\mathrm{1.623\hspace{0.17em}000}\text{E}-01$	$\mathrm{9.475\hspace{0.17em}000}\text{E}-01$	$\mathrm{6.773\hspace{0.17em}000}\text{E}-02$	$\mathrm{9.563\hspace{0.17em}000}\text{E}-01$	$\mathrm{8.021\hspace{0.17em}000}\text{E}-02$
$\#14$			$1\text{E}-3$	$\mathrm{2.368\hspace{0.17em}000}\text{E}-01$	$\mathrm{9.481\hspace{0.17em}000}\text{E}-01$	$\mathrm{1.607\hspace{0.17em}000}\text{E}-01$	$\mathrm{9.481\hspace{0.17em}000}\text{E}-01$	$\mathrm{5.402\hspace{0.17em}000}\text{E}-01$	–	$\mathrm{7.043\hspace{0.17em}000}\text{E}-01$
$\#15$	STRV-C2 (#1)	MP	$3\text{E}-1$	$\mathrm{3.764\hspace{0.17em}000}\text{E}-01$	$\mathrm{8.792\hspace{0.17em}667}\text{E}-01$	$\mathrm{2.553\hspace{0.17em}000}\text{E}-01$	$\mathrm{8.799\hspace{0.17em}333}\text{E}-01$	$\mathrm{3.657\hspace{0.17em}000}\text{E}-01$	–	$\mathrm{4.957\hspace{0.17em}000}\text{E}-01$
$\#16$			$2\text{E}-1$	$\mathrm{3.761\hspace{0.17em}000}\text{E}-01$	$\mathrm{8.794\hspace{0.17em}667}\text{E}-01$	$\mathrm{2.551\hspace{0.17em}000}\text{E}-01$	$\mathrm{8.801\hspace{0.17em}000}\text{E}-01$	$\mathrm{2.638\hspace{0.17em}000}\text{E}-01$	–	$\mathrm{5.253\hspace{0.17em}000}\text{E}-01$
$\#17$			$k_f^{\star }=1\text{E}-1$	$\mathrm{3.762\hspace{0.17em}000}\text{E}-01$	$\mathrm{8.796\hspace{0.17em}333}\text{E}-01$	$\mathrm{2.552\hspace{0.17em}000}\text{E}-01$	$\mathrm{8.801\hspace{0.17em}333}\text{E}-01$	$\mathrm{1.931\hspace{0.17em}000}\text{E}-01$	$\mathrm{8.341\hspace{0.17em}000}\text{E}-01$	$\mathrm{4.240\hspace{0.17em}000}\text{E}-01$
$\#18$			$1\text{E}-2$	$\mathrm{3.823\hspace{0.17em}000}\text{E}-01$	$\mathrm{8.747\hspace{0.17em}667}\text{E}-01$	$\mathrm{2.594\hspace{0.17em}000}\text{E}-01$	$\mathrm{8.750\hspace{0.17em}333}\text{E}-01$	$\mathrm{2.566\hspace{0.17em}000}\text{E}-01$	$\mathrm{6.857\hspace{0.17em}000}\text{E}-01$	$\mathrm{6.938\hspace{0.17em}000}\text{E}-01$
$\#19$		MA	$3\text{E}-1$	$\mathrm{2.408\hspace{0.17em}000}\text{E}-01$	$\mathrm{9.459\hspace{0.17em}333}\text{E}-01$	$\mathrm{1.629\hspace{0.17em}000}\text{E}-01$	$\mathrm{9.465\hspace{0.17em}667}\text{E}-01$	$\mathrm{4.197\hspace{0.17em}000}\text{E}-01$	–	$\mathrm{5.992\hspace{0.17em}000}\text{E}-01$
$\#20$			$2\text{E}-1$	$\mathrm{2.397\hspace{0.17em}000}\text{E}-01$	$\mathrm{9.462\hspace{0.17em}000}\text{E}-01$	$\mathrm{1.624\hspace{0.17em}000}\text{E}-01$	$\mathrm{9.467\hspace{0.17em}667}\text{E}-01$	$\mathrm{2.539\hspace{0.17em}000}\text{E}-01$	–	$\mathrm{7.955\hspace{0.17em}000}\text{E}-01$
$\#21$			$k_f^{\star }=1\text{E}-1$	$\mathrm{2.365\hspace{0.17em}000}\text{E}-01$	$\mathrm{9.477\hspace{0.17em}667}\text{E}-01$	$\mathrm{1.604\hspace{0.17em}000}\text{E}-01$	$\mathrm{9.480\hspace{0.17em}667}\text{E}-01$	$\mathrm{6.072\hspace{0.17em}000}\text{E}-02$	$\mathrm{9.697\hspace{0.17em}000}\text{E}-01$	$\mathrm{6.329\hspace{0.17em}000}\text{E}-02$
$\#22$			$1\text{E}-2$	$\mathrm{2.395\hspace{0.17em}000}\text{E}-01$	$\mathrm{9.462\hspace{0.17em}333}\text{E}-01$	$\mathrm{1.622\hspace{0.17em}000}\text{E}-01$	$\mathrm{9.468\hspace{0.17em}000}\text{E}-01$	$\mathrm{3.390\hspace{0.17em}000}\text{E}-01$	$\mathrm{6.209\hspace{0.17em}000}\text{E}-01$	$\mathrm{6.776\hspace{0.17em}000}\text{E}-01$

Equations

$$\text{mbox}{min}_{\mathbf{p},\mathbf{v}}\Vert \mathcal{R}\left(\mathcal{T}({\mathbf{m}}_A,\mathbf{p}),\mathbf{v}\right)-{\mathbf{m}}_D\Vert ,$$

$$\text{mbox}{min}_{\mathbf{p},\mathbf{v}}\Vert \mathcal{T}\left(\mathcal{R}({\mathbf{m}}_A,\mathbf{v}),\mathbf{p}\right)-{\mathbf{m}}_D\Vert .$$

$$\mathbf{m}(\mathbf{x},t)={\left(m_i(\mathbf{x},t)\right)}_{i=WM,GM,CSF}\in {\mathbf{R}}^3$$

$$\underset{\mathbf{v},\mathbf{p},k_f,w}{\min }{\mathcal{D}}_c[c(1),c_D]+{\mathcal{D}}_{\mathbf{m}}[{\mathbf{m}}_P(1)(1-c(1)),{\mathbf{m}}_D]+\mathcal{S}[\mathbf{p},\mathbf{v},w]$$

$$\mathcal{T}:$$

$$\Phi ({\mathbf{m}}_P(1))\mathbf{p}$$

$$\mathcal{R}:$$

$${\mathbf{m}}_A$$

$$\text{mbox}div\mathbf{v}$$

$$\mathbf{k}(\mathbf{m}):=k_f\sum _{i=1}^3{k}_i{m}_i\mathbf{I}=\overline{k}\mathbf{m}\mathbf{I}\text{mbox}and\rho (\mathbf{m}):={\rho }_f\sum _{i=1}^3{\rho }_i{m}_i=\overline{\rho }\mathbf{m}$$

$${\mathcal{S}}_{\mathbf{p}}^1[\mathbf{p}]=\frac{1}{2}\Vert \mathbf{p}{\Vert }_{L^1(\Omega )}^2\text{mbox}and{\mathcal{S}}_{\mathbf{p}}^2[\mathbf{p}]=\frac{1}{2}\Vert \mathbf{p}{\Vert }_{\text{W}(\Omega )}^2:=\sum _{i=1}^{p_n}{w}_i\cdot p_i^2,$$

$${\mathcal{S}}_{\mathbf{v}}[\mathbf{v}]=\frac{1}{2}\Vert \mathbf{v}{\Vert }_{H^1(\Omega )}^2=\frac{1}{2}{\int }_{\Omega }\sum _{i=1}^3|\nabla v^i(\mathbf{x}){|}^2\mathrm{d}\Omega .$$

$${\partial }_{t}c-\text{mbox}div\mathbf{k}({\mathbf{m}}_P(1))\nabla c-f(c,\rho ({\mathbf{m}}_P(1)))$$

$$c(0)-\Phi ({\mathbf{m}}_P(1))\mathbf{p}$$

$$-{\partial }_t\alpha -\text{mbox}div\mathbf{k}({\mathbf{m}}_P(1))\nabla \alpha +{\partial }_c{f}^{\star }(c,\rho ({\mathbf{m}}_P(1)))\alpha$$

$$c_D-c(1)+({\mathbf{m}}_P(1){)}^T({\mathbf{m}}_P(1)(1-c(1))-{\mathbf{m}}_D)-\alpha (1)$$

$${\partial }_t{\mathbf{m}}_P+\nabla {\mathbf{m}}_P\mathbf{v}$$

$${\mathbf{m}}_P(0)-{\mathbf{m}}_A$$

$$-{\partial }_t\mathbf{\lambda }-\text{mbox}div(\mathbf{\lambda }\otimes \mathbf{v})$$

$$(c(1)-1)({\mathbf{m}}_P(1)(1-c(1))-{\mathbf{m}}_D)+{\partial }_{{\mathbf{m}}_P(1)}\Phi (\mathbf{p})\alpha (0)-$$

$$\displaystyle\hbox to0.0pt{\hss$\displaystyle\int_{0}^{1}(\nabla c)^{T}\nabla\alpha\,\partial_{\boldsymbol{m}_{P}(1)}\boldsymbol{k}^{\star}(\boldsymbol{m}_{P}(1))+\partial_{\boldsymbol{m}_{P}(1)}f^{\star}(c,\rho(\boldsymbol{m}_{P}(1)))\alpha\,dt-\boldsymbol{\lambda}(1)$}=\boldsymbol{0}$$

$${\beta }_{\mathbf{p}}{\nabla }_{\mathbf{p}}{\mathcal{S}}_{\mathbf{p}}[\mathbf{p}]-{\Phi }^T\alpha (0)$$

$${\int }_0^1{\int }_{\Omega }{\mathbf{m}}_P(1)\left((\nabla c{)}^T\nabla \alpha \right)dxdt$$

$${\beta }_{\mathbf{v}}{\nabla }_{\mathbf{v}}{\mathcal{S}}_{\mathbf{v}}[\mathbf{v}]+\mathcal{K}[{\int }_0^1(\nabla {\mathbf{m}}_P{)}^T\mathbf{\lambda }\mathrm{d}t]$$

$${\partial }_{{\mathbf{m}}_P(1)}\mathbf{k}$$

$$(\boldsymbol{p}^{j+1},k_{f}^{j+1})=\operatorname*{arg\,min}_{\boldsymbol{p},k_{f}}~{}\underbrace{\mathcal{D}_{c}[c(1),c_{D}]+\mathcal{D}_{\boldsymbol{m}}[\boldsymbol{m}_{P}(1)^{j}(1-c(1)),\boldsymbol{m}_{D}]+\beta_{\boldsymbol{p}}\mathcal{S}_{\boldsymbol{p}}[\boldsymbol{p}]}_{\mbox{$=:\mathcal{J}_{T}$}}$$

$$({\mathbf{v}}^{j+1},w^{j+1})=\underset{\mathbf{v},w}{\mathrm{argmin}}{\mathcal{D}}_{\mathbf{m}}[{\mathbf{m}}_P(1)(1-c(1{)}^{j+1}),{\mathbf{m}}_D]+{\int }_{\Omega }{\mathbf{q}}^T{\mathbf{m}}_P(1)\mathrm{d}\mathbf{x}+{\beta }_{\mathbf{v}}{\mathcal{S}}_{\mathbf{v}}[\mathbf{v}]+{\beta }_w{\mathcal{S}}_w[w]$$

$$\mathbf{q}(\mathbf{x})={\int }_0^1\overline{k}((\nabla c(1{)}^{j+1}{)}^T\nabla {\alpha }^{j+1})+\overline{\rho }c(1{)}^{j+1}(1-c(1{)}^{j+1}){\alpha }^{j+1}\mathrm{d}t.$$

$$w_{i}=\left\{\begin{array}[]{ll}w^{\text{large}}\,\mbox{ if }\,|p_{i}^{\ell_{1}}|<\tau_{S}\|\boldsymbol{p}\|_{L^{\infty}(\Omega)}\\ w^{\text{small}}\,\mbox{ otherwise.}\end{array}\right.$$

$${\mathbf{p}}^{k+1}$$

$${\mathbf{p}}_{sd}^k$$

$${\mathbf{p}}^{k+1}$$

$$|{\mathcal{J}}_T({\mathbf{p}}^{k+1})-{\mathcal{J}}_T({\mathbf{p}}^k)|$$

$$\Vert {\mathbf{p}}^{k+1}-{\mathbf{p}}^k{\Vert }_{L^{\infty }(\Omega )}$$

$${\mathcal{H}}_s(\mathbf{p})=\frac{n_{\mathbf{p}}-\text{nicefrac}\Vert \mathbf{p}{\Vert }_{L^1(\Omega )}\Vert \mathbf{p}{\Vert }_{L^{\infty }(\Omega )}}{n_{\mathbf{p}}-1}$$

$$\forall \mathbf{x}\in \Omega :c^{(\cdot ,1)}(\mathbf{x})+\sum _{i=4}^4{m}_i^{(\cdot ,1)}(\mathbf{x})\stackrel{!}{=}1.$$

$${\mu }_{B,L^2}:=\frac{\Vert {\mathbf{m}}_P^{\prime }(1)-{\mathbf{m}}_D{\Vert }_{L^2(\Omega {)}^3}}{\Vert {\mathbf{m}}_A-{\mathbf{m}}_D{\Vert }_{L^2(\Omega {)}^3}},$$

$$H(\mathbf{u}):=\{u_i\ge 0.5\},$$

$${\text{DICE}}_{l,B}:=2\frac{|H(m_{P,l}(1))\cap H(m_{D,l})|}{|H(m_{A,l})|+|H(m_{D,l})|},$$

$$\Vert \mathbf{g}{\Vert }_{\text{rel}}:=\Vert {\mathbf{g}}^j{\Vert }_{L^2(\Omega )}/\Vert {\mathbf{g}}^0{\Vert }_{L^2(\Omega )}.$$

Full text

Preprint submitted to the SIAM Journal on Scientific Computing (SISC).

1 Introduction

Glioblastoma multiforme (GBM) tumor is a terminal primary brain cancer—the most aggressive one. Biophysical models are increasingly used to help the analysis of GBM MRI scans for epidemiological studies but also for assisting clinical decision making [54, 59], survival estimation, diagnosis, and preoperative and treatment planing [41, 47, 14, 55]. The key step in integrating biophysical models with clinical information is to calibrate them with patient MRI scans. After calibration, we can either use the estimated parameters as biomarkers, or we can evolve the calibrated PDE to estimate short-term tumor infiltration. Tumor growth models are typically nonlinear reaction-diffusion PDEs and their calibration is challenging. Not only do we need to estimate reaction and diffusion parameters, but also the tumor concentration, which is only partially (and implicitly) observed in the MRI scan. If we had an MRI scan of the patient without tumor, we could solve an inverse problem for the initial condition, the reaction and diffusivity parameters of the tumor. Unfortunately, a “healthy patient” MRI scan is rarely available. To resolve this conundrum, Davatzikos’s group [41] pioneered the idea of using an MRI scan of another, healthy, individual as a proxy for the tumor-free patient scan. In practice, we use a standardized average brain of several individuals, also known as a statistical atlas. But it turns out that naively using such an atlas would result in erroneous results. The second key idea addresses this by simultaneously using image registration to deform the atlas towards to the patient scan as far as possible. In this paper, we present a novel formulation of such a joint registration and inversion problem and a numerical solver scheme.

In the image registration problem, the goal is to estimate spatial point-correspondences between a template image $m_{T}$ (the atlas) and a reference image $m_{R}$ (the patient scan). To solve the registration problem, we use an optimal control formulation, in which we seek a (stationary) velocity field $\boldsymbol{v}$ (parameterizing a deformation map $\boldsymbol{y}$), such that the transported template image intensities match the intensities in the reference image, i.e., $m_{T}\circ\boldsymbol{y}\approx m_{R}$. In the tumor inversion problem, we want to estimate tumor-growth parameters (such as tumor origin, infiltration and proliferation rates) of our PDE model so that if we start growing the tumor at $t=0$ we match the partial tumor observations at $t=1$. (The astute reader is probably wondering how do we know the time horizon. We don’t, but $t=1$ is related to a non-dimensional form of the tumor growth PDE.) We simultaneously solve for both $\boldsymbol{v}$ and the tumor parameters.

Let us try to explain the setting using a somewhat high-level notation. In our work we are not using the original MRI scan intensities, but assume that we have segmentation labels describing different brain tissue (white matter, gray matter etc.) distributions for both atlas and patient images. Let $\boldsymbol{m}_{A}$ be the healthy atlas labels and $\boldsymbol{m}_{D}$ the patient labels (both vector functions in the unit cube in $\mathbb{R}^{3}$). Let $\mathcal{T}$ and $\mathcal{R}$, respectively, be abstract forward operators for a tumor simulation and registration mapping component. $\mathcal{T}(\boldsymbol{m},\boldsymbol{p})$ takes a label image $\boldsymbol{m}$ without tumor and tumor model parameters $\boldsymbol{p}$ and creates a labeled image that now has tumor labels in addition to the original healthy tissue labels (e.g., gray matter, white matter, etc.). $\mathcal{R}(\boldsymbol{m},\boldsymbol{v})$ takes a label image $\boldsymbol{m}$ and a velocity $\boldsymbol{v}$ that parameterizes the deformation and creates a deformed image.

Then, we can summarize the general idea of our previous Moving-Patient and the new Moving-Atlas scheme by defining two optimization problems:

[TABLE]

Notice that the main difference is switching the order of the registration and the tumor operators; see Fig. 1 for an illustration. This seemingly simple change has significant impact in the solution of this problem. In [49, 48], we used the Moving-Patient scheme and only inverted for the tumor initial condition. Note that in MP, to evaluate the objective function, we first grow a tumor in the atlas and then deform the patient image to match the resulting tumor-bearing atlas. As long as the tumor-bearing atlas is topologically similar to the tumor-bearing patient brain, image registration can always yield good tumor reconstruction, even with completely wrong tumor parameters. Thus this approach enables efficient non-diffeomorphic registration between an atlas and a tumor-bearing patient, but is limited in terms of meaningful tumor inversion.

In this paper, we propose the Moving Atlas scheme. To evaluate the objective function given, we first deform the healthy atlas (using $\boldsymbol{v}$) so that it matches the patient. Then we grow a tumor in the transformed atlas. It turns out that the new formulation is more appropriate for the biophysical modeling since the tumor growth takes place in images that resemble the actual patient.

Contributions.   The main contributions of this paper are as follows:

1. (i)

We present a new optimization problem formulation (Moving Atlas) for tumor growth model calibration based on patient individual single-snapshot data. 2. (ii)

The new formulation prevents fitting of patient input data to a possibly poor tumor reconstruction. The computed biophysical parameters “live” in the patient space, i.e., the tumor parameters are estimated assuming tumor growth in (an approximation of) the healthy patient brain (as opposed to the atlas brain for the former scheme). This renders the new scheme better for biophysical inversion.111The presented results are to be seen as a proof of concept for the developed methodology. To enable predictive capabilities, we require a more complex tumor model. 3. (iii)

We derive a Picard-iteration-type solution scheme that alternates between the image registration and the inverse tumor-growth problem.222A modified objective function for the registration sub-component allows to fulfill the strongly coupled first order optimality conditions of the joint optimization problem formulation 4. (iv)

We enhance the tumor inversion component. The tumor solver used in [49] inverted for the tumor initial condition only, penalized with an $\ell_{2}$ constraint. In this paper, we use an $\ell_{1}$ constraint (which restricts the tumor to more plausible initial condition), and also invert for the diffusion parameter in the tumor growth PDE (modeling tumor infiltration). 5. (v)

We conduct numerical experiments to evaluate the nem scheme using synthetically grown tumor data from real clinical brain imaging data. We demonstrate that the Moving Atlas scheme yields better results in terms of the accuracy and robustness of tumor parameter reconstruction.333For all considered test cases the tumor growth ground truth parameters are known, to allow for comparison of the respective reconstruction accuracy. For a synthetic example, we showcase malfunction444The patient input data is fitted to a poor reconstruction of the tumor in the atlas (due to large anatomical differences); consequently, the target data is corrupted in the further course of the inversion. of the Moving Patient scheme, whereas the new Moving Atlas scheme scheme yields good results.

Limitations.   The reaction-diffusion tumor-growth model we consider is purely phenomenological. Important biophysical, bio-mechanical and biochemical effects and phenomena such as mitosis, apoptosis, chemotaxis, deformation of brain parenchyma (mass-effect), and the modeling of edema, necrosis, and angiogenesis are neglected. It is, however, the most common model used for clinical analysis [53, 55, 5, 15, 22, 17]. The effect of tumor growth induced deformation of brain parenchyma is important and not accounted for in our current model. This is ongoing work. Note that the main work here lies in the enhancement of the tumor component. Integrating it in the coupled formaulation and the Picard iteration is straight-forward due to our modular approach. The inverse tumor-growth problem with inversion for all growth parameters is ill-posed. We have to introduce additional prior information. Currently, if we want to invert for the proliferation rate $\rho$ and the diffusion rate, we run the solver for multiple values of $\rho$ and compare the quality of reconstruction results. We do not account for uncertainties in input imaging data, employed algorithms and solution approaches, and model parameters (such as tumor initial condition and characteristic net cell migration into surrounding tissue). A Bayesian framework is subject to future work; we currently invert for the maximum a posteriori estimate of the parameters. Our approach does not directly process MRI data but requires prior segmentation of the imaging data. Note that reliably integrating models with MRI is an open problem, since one needs to estimate cell-density from data or link model output to intensities. Lastly, we emphasize that within this work we present a proof-of-concept for a new, promising methodology; future work will include its application to real clinical data.

Related Work.   One important application of biophysical models is to enable non-diffeomorphic image registration. Classical image registration [39, 50] assumes that the input images (template and reference) are topologically equivalent. But registering a healthy atlas (no tumor) to a tumor-bearing patient image violates this assumption. We refer to this problem as the normal-to-abnormal registration. We refer to our previous work in [49] for an extensive review in non-diffeomorphic image registration.

Here, we are interested in fitting models to patient-individual imaging data. Several groups have tackled this problem using derivative-free optimization approaches [4, 20, 22, 29, 35, 38, 57, 58], or address the parameter estimation problem within a derivative-free Bayesian framework [16, 24, 25, 27, 26, 37, 45]. There has been only limited work in the direction of PDE-constrained optimization for model-based image analysis [6, 17, 10, 20, 21, 28, 29, 36, 46]. Most of these approaches rely on longitudinal patient data. In particular, they require knowledge of the healthy patient brain before tumor occurrence. To allow for model-inversion based on single-snapshot data only, deformable inter-subject registration becomes necessary to artificially generate a second snap-shot in time. Such a combination of registration and biophysical inversion has been previously targeted in [14, 17, 41, 61, 60, 62]. The work in [41, 61, 60, 62], however, uses a purely mechanical model for tumor progression and falls short in providing information about progression and infiltration of cancerous cells into surrounding healthy tissue. In [2, 23, 17, 14], the authors propose a framework for joint segmentation, registration and tumor modeling, which is very similar to our approach. In our previous work [49], we significantly reduce the time-to-solution by employing second-order derivative information in conjunction with highly scalable and efficient numerics, improved algorithms and powerful preconditioners. We eliminate the need for manual tumor seeding, and improve inter-subject registration performance.

Outline.   The outline of the paper is as follows: In §2 we present the mathematical formulation for the Moving Atlas joint registration and biophysical inversion problem. We present the PDE-constrained optimization problem and derive the first order optimality conditions in §2.1 and propose a Picard-iteration-type solution strategy, outlined in §2.2. Details on continuation schemes, adaptivity and convergence criteria are given in §2.3. More details on the two main sub-components, diffeomorphic medical image registration and biophysical tumor-growth inversion are summarized in §3. In §4, we perform numerical experiments to analyze our scheme, and compare it to the Moving Patient scheme [49].

2 Formulation and Picard Iteration for the Moving Atlas Coupled Problem

To apply tumor inversion and image registration in a joint approach, we consider an optimal control formulation, which results in a PDE-constrained, nonlinear optimization problem. In this section, we present the new Moving Atlas formulation along with an iterative fixed-point coupling scheme. This scheme allows us to solve the joint optimization problem based on two separate components: a tumor-growth inversion solver and a modified diffeomorphic image registration component. Before presenting the formulation and the iterative solver, we shortly introduce the notation we are going to use throughout the paper.

Notation.   $m_{X}(\boldsymbol{x},t)\in[0,1]$, $(\boldsymbol{x},t)\in\Omega\times[0,1]$, with $X\in\{WM,GM,CSF\}$ represents probability maps of different brain tissue types, namely, white matter (WM), gray matter (GM), and cerebrospinal fluid (CSF; which includes ventricles). The probability maps for each tissue label are computed from (segmented) MR imaging data in a pre-processing step. In the following, we refer to $m_{X}(\boldsymbol{x},t)$ as anatomy labels. $\Omega=[0,2\pi]^{3}$ denotes the normalized spatial domain. We gather these probability maps in a space-time vector field

[TABLE]

The normalized tumor cell concentration $c\colon\Omega_{B}\times[0,1]\rightarrow[0,1],(\boldsymbol{x},t)\mapsto c(\boldsymbol{x},t)$ (tumor map) is given as a fourth field and interpreted as probability map for cancerous tissue. $\Omega_{B}\subset\Omega=[0,2\pi)^{3}$ is embedded555We use a penalty approach to discretize the tumor equations in $\Omega$ and approximate the Neumann boundary conditions on $\partial\Omega_{B}$. in the simulation domain $\Omega$ and denotes the domain occupied by brain tissue; it is bounded by the skull. All anatomy and tumor labels evolve in space and time. To simplify notation, we introduce the space-time domain $U\mathrel{\mathop{:}}=\Omega\times(0,1]$, and $\bar{U}\mathrel{\mathop{:}}=\Omega\times[0,1)$. We omit the spatial dependency, and indicate temporal evolution by the pseudo time $t_{R}\in[0,1]$ for the advection problem of the registration, and the normalized time $t_{T}\in[0,1]$ for tumor growth, cf. Tab. 1. To judge the proximity of the predicted state of our mathematical model to the observed data $\boldsymbol{m}_{D}$, $c_{D}$, we define the $\ell_{2}$-distance measures $\mathcal{D}_{c}$ and $\mathcal{D}_{\boldsymbol{m}}$ (see Tab. 1), which drive the optimization process.

2.1 Moving Atlas Coupled Formulation

The Moving Atlas coupled problem formulation is given by the minimization problem

[TABLE]

That is, we seek parameters and space time fields $(\boldsymbol{v},\boldsymbol{p},k_{f},w)$, such that the predicted state of our model becomes similar to the observed data. The distance between data and predicted state is quantified by the two $\ell_{2}$-distance measures $\mathcal{D}_{c}[c(1),c_{D}]$ and $\mathcal{D}_{\boldsymbol{m}}[\boldsymbol{m}_{P}(1)(1-c(1)),\boldsymbol{m}_{D}]$. $\mathcal{D}_{c}[c(1),c_{D}]$ measures the discrepancy between the simulated tumor $c(1)$ computed from (3a)–(3b) in the approximated patient space and the input tumor data $c_{D}$. $\mathcal{D}_{\boldsymbol{m}}[\boldsymbol{m}_{P}(1)(1-c(1)),\boldsymbol{m}_{D}]$ quantifies the discrepancy between the patient’s anatomy labels and the warped-to-patient atlas anatomy labels in regions not affected by the tumor666The simulated pathology is embedded into the advected atlas anatomy labels when caculating the registration misfit. This can be seen as a “masking” of the tumor region, i.e., the registration does not alter the tumor concentration nor the discrepancy measure $\mathcal{D}_{c}$. computed from (3c)–(3e)). $\mathcal{S}[\boldsymbol{p},\boldsymbol{v},w]\mathrel{\mathop{:}}=\beta_{\boldsymbol{p}}\mathcal{S}_{\boldsymbol{p}}[\boldsymbol{p}]+\beta_{\boldsymbol{v}}\mathcal{S}_{\boldsymbol{v}}[\boldsymbol{v}]+\beta_{w}\mathcal{S}_{w}[w]$ is a regularization term.

To simulate tumor-growth over time, we use a diffusion-reaction model with diffusion tensor $\boldsymbol{k}$ and a logistic reaction term $f(c,\rho)=\rho c(1-c)$ with reaction rate $\rho$. For the tumor initial conditions, we use a parametrization $\Phi\boldsymbol{p}$, i.e., $c(0)$ lives in an $n_{\boldsymbol{p}}$-dimensional space spanned by Gaussian basis functions ($\boldsymbol{p}\in\mathbf{R}^{n_{\boldsymbol{p}}}$, $\Phi\boldsymbol{p}:=\sum_{i=1}^{n_{\boldsymbol{p}}}\Phi_{i}p_{i}$). Within CSF, the Gaussians are set to zero to prevent spurious diffusion of cancerous cells into CSF. We model the diffusion and reaction coefficient as

[TABLE]

with $\overline{k}=(k_{1}k_{f},k_{2}k_{f},k_{3}k_{f})^{T}$ and $\overline{\rho}=(\rho_{1}\rho_{f},\rho_{2}\rho_{f},\rho_{3}\rho_{f})^{T}$, i.e., both diffusion and reaction parameters vary in space depending on the value of the anatomy labels $\boldsymbol{m}$ describing the brain geometry. For simplicity, we only consider isotropic diffusion777Note, that our approach can be generalized in a straight-forward way to general diffusion tensors. To regularize the tumor growth initial condition parametrization, we use combinations of $\ell_{1}$- and $\ell_{2}$-regularization, i.e.,

[TABLE]

where $\beta_{\boldsymbol{p}}>0$ is the regularization parameter and $\|\cdot\|_{\textbf{W}(\Omega)}$ is a weighted $\ell_{2}$-norm. For more details on how we switch between the two penalty terms in the solver algorithm, see §2.3.

For image registration, we use an advection model with a stationary velocity field $\boldsymbol{v}$. The image registration velocity is regularized via the $H^{1}$ Sobolev semi-norm:

[TABLE]

$\mathcal{S}_{w}[w]=\|w\|_{L^{2}(\Omega)}$ is a regularization term for $w=\mbox{div}\,{\boldsymbol{v}}$, which allows us to control volume change. We use periodic boundary conditions on $\partial\Omega$.

Optimality Conditions.   As a general optimization approach, we use the method of Lagrange multipliers. Taking variations yields the strong form of the first order optimality conditions:

[TABLE]

with adjoint variables $\alpha$, $\boldsymbol{\lambda}$, $\nu$ associated to the state variables $c(1)$, $\boldsymbol{m}_{P}$, and $\boldsymbol{m}_{P}(1)$. We eliminate $w$ from the optimization problem and use the resulting operator $\mathcal{K}$ in (7k). For $\mbox{div}\,\boldsymbol{v}=w=0$, $\mathcal{K}$ would be the Leary projection \mathcal{K}(\boldsymbol{u})\mathrel{\mathop{:}}=\boldsymbol{u}+\nabla\rotatebox[origin={c}]{180.0}{\nabla}^{-1}\mbox{div}\,\boldsymbol{u}; for a non-zero $w$, the projection operator becomes slightly more involved; we refer to [30, 31] for additional details. The gradients of the diffusion and reaction terms in (7h) with respect to $\boldsymbol{m}_{P}(1)$ can be derived from (4), and are given by:

[TABLE]

Discussion.   We summarize the advantages of the moving atlas scheme over the previously used moving patient scheme [49]:

(i) image registration operates on the brain anatomy only such that the tumor solver is forced to actually produce a biophysically more meaningful888Assuming a biophysical tumor progression model which accurately reflects the physiological processes; the model used here does not fulfill this assumption. tumor that is similar to the observed patient tumor.

(ii) Tumor growth parameters, i.e., $\boldsymbol{p}$ and $k_{f}$, are estimated directly in an approximation of the patient brain such that their quality does not strongly deteriorate with differences between the atlas and the patient brain.

This comes at a price, for which we account for by employing suitable solver strategies presented in §3:

(a) the decomposition of the moving atlas formulation into tumor and registration components requires modifications of the registration formulation as shown in §2.2,

(b) exploiting the full potential of the moving atlas scheme requires sparse localization of the tumor initial condition, encouraged by an $\ell_{1}$-regularization in the tumor component as outlined in §3.1 (cf. also [52, 48]).

2.2 Picard Iteration Scheme

We solve the optimization problem by iterating over the registration velocity $\boldsymbol{v}$, the tumor growth initial guess parametrization $\boldsymbol{p}$ and (if switched on) the diffusion coefficient $k_{f}$. We use a superscript $j$ to mark variables associated with the $j$th Picard iteration. Each Picard iteration is decomposed two into sub-steps:

1. (1)

Solve inverse tumor problem, i.e., given $\boldsymbol{v}^{j}$ solve for $(\boldsymbol{p}^{j+1},k_{f}^{j+1})$ such that

[TABLE]

subject to the tumor forward problem (3a) and (3b) for the healthy patient anatomy $\boldsymbol{m}_{P}(1)=\boldsymbol{m}_{P}(1)^{j}$ kept fixed (result of the previous Picard iteration or initial guess). This means, we fulfill the first order optimality conditions (7a)–(7d), (7i) and (7j) for $\boldsymbol{m}_{P}(1)=\boldsymbol{m}_{P}(1)^{j}$. Aside from new iterates for $\boldsymbol{p}$ and $k_{f}$, this step yields a new simulated tumor $c(1)^{j+1}$ and tumor adjoint $\alpha^{j+1}$. 2. (2)

Solve modified registration problem, i.e., given $(\boldsymbol{p}^{j+1},k_{f}^{j+1})$, solve for $\boldsymbol{v}^{j+1}$ such that

[TABLE]

subject to the image advection problem (3c)–(3e). $\boldsymbol{q}=(q_{1},q_{2},q_{3})^{T}$ is defined as

[TABLE] 3. (3)

iterate until convergence.

The particular modification of (10) by introducing the term $\int_{\Omega}\boldsymbol{q}^{T}\boldsymbol{m}_{P}(1)\,\mathop{}\!\mathrm{d}\boldsymbol{x}$ ensures that the registration problem reproduces (7e) – (7h) and (7) as first order optimality conditions. Note that we freeze the variables $c(1)$, and $\alpha$ from the previous iteration, i.e., neglect the indirect impact of the registration velocity $\boldsymbol{v}$ on $c(1)$ and $\alpha$ via changes in $\boldsymbol{m}_{P}(1)$. Besides new iterates for $\boldsymbol{v}$, the registration step also provides a new approximation $\boldsymbol{m}_{P}(1)^{j+1}$ of the healthy patient brain. The converged solution of the respective Picard or fixed-point iteration fulfills all first order optimality conditions of (3) with the exception of (7h). For the latter, the term $\partial_{\boldsymbol{m}_{P}(1)}\Phi(\boldsymbol{p})\alpha(0)$ is neglected and not enforced explicitly. We found experimentally, that our scheme works well and produces satisfying results despite the simplification.

2.3 Adaptivity and Convergence Criteria

In our Picard scheme, we use several strategies to control convergence and to adapt the tumor and image registration regularization over the Picard iterations.

For the tumor solver, the regularization term $\beta_{\boldsymbol{p}}\mathcal{S}_{\boldsymbol{p}}[\boldsymbol{p}]$ is used to attain better conditioning, but also to encourage sparse localization of the initial condition. Thus, we use $\mathcal{S}_{\boldsymbol{p}}[\boldsymbol{p}]=\nicefrac{{1}}{{2}}\|\boldsymbol{p}\|_{L^{1}(\Omega)}^{2}$ to compute a sparse solution of the initial condition and switch to a weighted $\ell_{2}$-regularization $\mathcal{S}_{\boldsymbol{p}}[\boldsymbol{p}]=\|\boldsymbol{p}\|^{2}_{\boldsymbol{W}(\Omega)}=\nicefrac{{1}}{{2}}\sum_{i=1}^{p_{n}}w_{i}\cdot p_{i}^{2}$ in a second phase of the tumor inversion solver, where the weights $w_{i}$ are calculated according to

[TABLE]

Here, $\tau_{S}$ is a user defined tolerance. The $\ell_{2}$-phase uses $\boldsymbol{p}^{\ell_{1}}$ achieved in the $\ell_{1}$-phase as initial approximation. In the $\ell_{1}$-phase we invert for $\boldsymbol{p}$ only and keep the diffusivity $k_{f}$ (estimated in the $\ell_{2}$-phase) fixed;999We found experimentally, that diffusion inversion in the $\ell_{1}$-phase has vanishing impact on the proposed sparsity of the solution, and the $\ell_{2}$ estimation for $k_{f}$ is better. initially, we start with $k_{f}=0$. For more details on numerics of the tumor inversion solver, we refer to [49].

For the registration solver, convergence and robustness are improved by regularization. We perform a parameter continuation. That is, we start with a large value $\beta_{\boldsymbol{v}}^{\text{init}}$ and subsequently reduce $\beta_{\boldsymbol{v}}$ by a factor of ten after each Picard iteration until we reach the prescribed value $\beta_{\boldsymbol{v}}^{\text{final}}$. Image registration calculates solutions and possibly increases $\beta_{\boldsymbol{v}}$ internally based on a check of lower admissible bounds for the determinant of the deformation gradient. The deviation of the determinant of the deformation gradient from one indicates local volume changes (expansion for a value larger then one and contraction for a value smaller then one; see [31, 49] for more details).

As a stopping criterion for the Picard iterations, we consider a mixed criterion that is fulfilled either if we reach $\beta_{\boldsymbol{v}}^{\text{final}}$ or if the image registration fails due to violation of the user defined lower bound on the determinant of the deformation gradient. We always execute two Picard iterations with the final $\beta_{\boldsymbol{v}}^{\text{final}}$.

For both tumor inversion and image registration, we use fixed tolerances for the norm of the reduced gradients as given in (7i), (7j), and (7k), respectively. Note that we use norms relative to the gradient norm at the beginning of the first Picard iteration for all gradients. We combine these gradient tolerances with a prescribed maximal number of Newton iterations for the reduced space KKT systems (see also [49]).

3 Numerical Methods for Tumor Inversion and Image Registration

In this section, we shortly describe the numerical and algorithmic features of the tumor and the image registration solvers. We give only a brief summary along with references to previous publications since these methods have not been changed from our previous work.

3.1 Regularization and Nonlinear Solver Components for Tumor Inversion

The tumor inversion problem is solved based on an optimize-then-discretize approach. To solve the respective system of optimality conditions, we use different solvers for $\ell_{1}$ or weighted $\ell_{2}$-regularization, respectively. For $\ell_{1}$-regularization, the General Iterative Shrinkage Thresholding (GIST) [13] is applied. In every iteration, it determines a solution $\boldsymbol{p}^{k+1}$ close to the result $\boldsymbol{p}_{sd}^{k}$ of a gradient descent step in the direction $\nabla\mathcal{D}_{c}(c(1),c_{D})$ with step length $\alpha^{k}$:

[TABLE]

The stopping criteria are given by

[TABLE]

We perform parameter continuation in $\beta_{\boldsymbol{p}}$ to estimate its value. We start with a large $\beta_{\boldsymbol{p}}$, bounded by $\|\nabla\mathcal{D}_{c}(c(1),c_{D})\|_{L^{\infty}(\Omega)}$ and perform a binary search in subsequent GIST iterations based on the Hoyer sparsity measure [19], $\mathcal{H}_{s}(\boldsymbol{p})\in[0,1]$:

[TABLE]

where $n_{\boldsymbol{p}}$ is the dimensionality of $\boldsymbol{p}$. Solutions with larger values of $\mathcal{H}_{s}$ feature higher sparsity, i.e., have less non-zero entries. We use a tolerance of $0.95$ to identify sparse solutions. If the sparsity is large enough (larger than the tolerance), we search for a smaller regularization parameter and, alternatively, increase it if the sparsity is insufficient. GIST features slow first order convergence, which is, however, acceptable for our purpose as only a rough estimate for $c(0)$ is to be calculated with $\ell_{1}$-regularization.

In the $\ell_{2}$-phase, we employ a reduced space method. This means that we assume the state and adjoint equations to be fulfilled exactly; they can be eliminated from the KKT system. For the tumor inversion solver, the reduced gradients are calculated based on the current iterates for $\boldsymbol{p}$ and $k_{f}$ as follows:

(i) solve the state equations (7a) – (7b),

(ii) use the result to solve the adjoint equations (7c) – (7d),

(iii) insert the result in (7i) – (7j) to calculate the reduced gradient with respect to $\boldsymbol{p}$ and $k_{f}$.

As a solver for the reduced KKT system, we use a quasi-Newton approach (LBFGS, [44, p. 135ff]) yielding a matrix-free approximation of the matrix-vector product of the gradient with the respective inverse Hessian and, thus, making an inner linear solver obsolete. We use Moré-Thuente line-search [42] for globalization.101010Note that we do not use line-search in the outer Picard iterations. The $\ell_{2}$-phase is terminated when the relative change of the norm of the reduced gradient is below a user defined threshold $\text{opttol}_{T}>0$. The reference gradient is the gradient obtained for the zero initial guess, $\boldsymbol{p}^{0}$.

3.2 Nonlinear Solver Components for Image Registration

We use an optimize-then-discretize approach for the image registration solver [34]. We employ a reduced space method: The current iterate for the velocity $\boldsymbol{v}$ is used to

(i) solve the state equations (7e) – (7f) (forward advection),

(ii) use the result to solve the adjoint equations (7g) – (7h),

(iii) insert the result in (7k) to calculate the reduced gradient.

We use a Newton-Krylov approach with an inexact Gauss-Newton linearization and a matrix-free PCG method as inner linear solver (see [30, 31, 11, 33, 32] for more details) to compute the search direction.

For the convergence of the optimizer, we use a combination of the relative change of

(i) the norm of the gradient in (7k), (ii) the objective in (10) and (iii) the control variable $\boldsymbol{v}$, all controlled by a single parameter $\text{opttol}_{R}>0$, as a stopping criterion. More details on the stopping conditions can be found in [30, 40]; see [12, 305 ff.] for a discussion. For the inner PCG solver, we perform inexact solves [7, 8] with a tolerance that is proportional to the norm of the reduced gradient of our problem. This prevents over-solving of the Hessian system if far from the optimum (large gradient norm). See also [44, p. 165ff.] for details. We use Armijo line-search for globalization.

For both tumor inversion and registration solver, we specify a maximum number ($\text{maxit}_{N,T}$ / $\text{maxit}_{N,R}$) of Newton, a maximum number ($\text{maxit}_{K,T}$ / $\text{maxit}_{K,R}$) of Krylov iterations, and a lower bound of $110-6$ for the absolute norm of the gradient as a safeguard against a prohibitively high number of iterations.

3.3 Discretization in Space and Time for the State and Adjoint Equations

For all state and adjoint equations in (7), we use the following numerical ingredients that have been published in more detail in  [30, 31, 11, 10, 1, 33, 32]:

(i) spectral elements for spatial discretization combined with a regular mesh;

(ii) 3D Fourier transforms to compute spatial derivatives;

(iii) an unconditionally stable semi-Lagrangian time-stepping scheme to avoid stability issues and small time-steps for image advection;

(iv) an approximation of boundary conditions for the tumor equations at the brain surface based on a penalty approach combined with periodic boundary conditions at the boundary of the overall domain $\bar{\Omega}=[0,2\pi]^{3}$;

(v) an unconditionally stable, second order Strang-splitting approach [51, 18] with analytical solution of the reaction terms and an implicit Crank-Nicholson method for the diffusion terms in the tumor equations.

4 Numerical Results

The experiments serve as a proof of concept for the newly developed joint registration and biophysical inversion strategy (Moving Atlas) and demonstrate that the previously introduced moving patient scheme [49] can fail dramatically in terms of both tumor reconstruction and inversion for initial conditions and tumor growth parameters. We quantitatively and qualitatively compare the new Moving Atlas strategy to the Moving Patient with an emphasis on model inversion and parameter estimation and address the scheme’s eligibility to recover true model parameters. This is an important step towards reconstruction of meaningful biophysical characteristics. We consider the following classes of test cases:

SYN – Synthetic Brain.   This first academic test case shows that the moving patient scheme can dramatically fail in terms of biophysical inversion. More precisely, we demonstrate that the moving patient scheme yields largely wrong reconstruction of the tumor, the tumor initial condition, and the diffusion coefficient. We use two simple elliptical anatomies consisting of white matter and gray matter only for atlas and healthy patient brain, respectively, and grow a tumor with known diffusion and reaction parameters in the patient anatomy. The resulting pathologic patient brain serves as target input data for our joint inversion approach. In order to achieve significant results, we create a case with large inter-subject differences between atlas and patient anatomy.

STRV – Synthetic Tumor, Real Velocity.   This second class of test cases is based on real brain geometries and serves as a means to examine the potential of the moving atlas formulation for biophysical inversion in realistic geometries, i.e, to examine the quality of reconstruction for the known ground truth initial conditions, $k_{f}$ and healthy patient geometry. We use two different real healthy brain geometries and grow a synthetic tumor in one of them, which we use as the patient brain. The second healthy brain serves as a normal brain template (atlas brain). In contrast to real patient tumors, in this case we know the ground truth for tumor initial conditions, diffusion coefficient $k_{f}$ and healthy patient geometry. Note that we do not know the ground truth for the registration velocity $\boldsymbol{v}$.

In our test cases, we compare three different solution strategies:

I. Tumor Standalone (T).   For the tumor standalone solution, the inversion is carried out assuming a statistical atlas brain geometry as healthy patient brain. No inter-subject registration is applied. The error is a function of the inter-subject variability between atlas and patient brain.

II. The Joint Registration and Biophysical Inversion Moving Patient Scheme (MP).   The Moving Patient scheme presented in [49] exhibits some inherent shortcomings with respect to the suitability for biophysical inversion and model calibration. The registration may fit the input data to a poor tumor reconstruction. We use the method described in [49] enhanced with a sparsity constraint for the tumor inversion solver §3.1.

III. The Joint Registration and Biophysical Inversion Moving Atlas Scheme (MA).   The Moving Atlas scheme presented in this work is designed to remedy some of the shortcomings of the Moving Patient scheme.

In the following, we shortly describe data and parameters for the experiments in §4.1 as well as the hardware used and the general setup in §4.2. We introduce performance measures used in the evaluation of our numerical results in §4.3. In §4.4, and §4.5 we show results for the two classes of test cases.

4.1 Data and Parameters

We use a common data basis for brain geometries as well as some fixed parameter settings for all test cases as described below.

Brain Imaging Data.   For the STRV cases, we use normal brain MR imaging data obtained at the Perelman School of Medicine at the University of Pennsylvania. The three-dimensional imaging data have an image resolution of $256^{3}$ voxels. We use binary segmentations of MRI scans for white matter (WM), gray matter (GM), and cerebrospinal fluid with ventricles (CSF). In a pre-processing step, these labels are smoothed and rescaled to probability maps, i.e, values between zero and one. To ensure partition of unity across all anatomy labels for each $\boldsymbol{x}$ in $\Omega$, we introduce a fourth probability map for background based on

[TABLE]

Common Model and Numerical Parameters.   The tumor model and solver parameters are summarized in Tab. 2. The solver parameters are based on experiments. We briefly discuss some of the more involved choices: For the parametrization of the tumor initial condition, we select a regular grid of $n_{\boldsymbol{p}}$ Gaussian basis functions with a variable standard variation $\sigma\in\{\nicefrac{{2\pi}}{{30}},\nicefrac{{2\pi}}{{64}}\}$ and a grid spacing of $\delta=1.5\sigma$ around the center of mass of the tumor. The number $n_{\boldsymbol{p}}$ is chosen a priori on a case-by-case basis, such that the pathological domain is covered sufficiently.

For all experiments, we consider tumor growth only in white matter ($k_{1}=1;\rho_{1}=1$), and, consequently, set the characteristic diffusivity and net cell proliferation in gray matter and CSF to zero, i.e., $k_{2}=k_{3}=0$, and $\rho_{2}=\rho_{3}=0$. To grow the synthetic tumor, we choose $\rho_{f}^{\star}=15$ as cell proliferation rate (reaction coefficient), and $k_{f}^{\star}=$110-1 as rate of cell migration into surrounding tissue (diffusion coefficient), respectively. For the tumor evolution simulation, we choose a time step of $\Delta t=0.01$ and various time horizons $T\in\{0.16,0.32,0.44\}$. The regularization parameter $\beta_{\boldsymbol{p}}$ for the $\ell_{2}$-regularized inverse tumor problem has been determined experimentally from an L-curve study similar to [[10](#bib.bib10)], using $n_{\boldsymbol{p}}=125$ Gaussian basis functions and an image resolution of $n=128^{3}$. We fix the relation $\nicefrac{{\sigma}}{{\delta}}=1.5$ for the spacing of the Gaussian basis functions111111This leads to an invariant condition number of $\Phi^{T}\Phi$ for the corresponding interpolation problem for varying standard variations $\sigma$, i.e., width of the Gaussians.. The lower bound $\beta_{\boldsymbol{v}}^{final}$ on the regularization parameter $\beta_{\boldsymbol{v}}$ has been determined based on extensive numerical experiments for different synthetic and real brain data sets [[31](#bib.bib31), [32](#bib.bib32)]. By experimental analysis we found that reducing the regularization weight by a factor of $10$ in every Picard-iteration is sufficient121212Using a smaller reduction or more Picard-iterations per fixed regularization weight did not improve the overall result.. For the optimizer, we use tolerances of $\mbox{opttol}_{R}=\mbox{opttol}_{T}=$110-3 for registration and tumor inversion, respectively131313Note that, in combination with the applied stopping conditions, this results in a required relative gradient of $110-3$ for the tumor inversion, but a gradient reduction of only about two orders of magnitude for the registration.. Further reduction of the gradient did not improve the final reconstruction quality, since the overall error is bound by an $\mathcal{O}(1E-3)$ error introduced by the solution of the hyperbolic advection problem (cf. [49]). To prevent prohibitively large runtimes, we define upper bounds on the number of Newton (and Krylov) iterations.

4.2 Hardware and Setup

All numerical experiments were executed on the supercomputer HazelHen  at the High Performance Computing Center HLRS in Stuttgart (www.hlrs.de), a Cray XC40 system with a peak performance of $7.42$ Petaflops comprising $7,712$ nodes with Xeon E5-2680 v3 processors and $24$ cores on two sockets and 128 GB memory per node. Our SIBIA framework is written in C++ and uses MPI for parallelism. It is compiled using the Intel 17 compiler. We use PETSc’s implementations for linear algebra operations and PETSc’s TAO package for the nonlinear optimization [3, 43], AccFFT for Fourier transforms [10, 9], and PnetCDF for I/O [56]. We use $3$ nodes with $64$ MPI tasks for (down-sampled) data sizes of size $n_{i}=128$, $i=1,2,3$ for all runs reported in this study.

4.3 Performance Metrics

The task of our numerical experiments is to assess

(i) the convergence towards solutions with low mismatch between the prediction of our model in (3a)–(3e) and the observed patient input data, both in the brain geometry $\boldsymbol{m}$ and the tumor $c$ (for all test cases),

(ii) the reconstruction quality for the (in practice inaccessible) healthy patient brain141414For our semi-synthetic test case setting, where the patient target data is generated from synthetic tumor progression simulation in the healthy brain, the healthy patient brain is known and the approximation quality can be measured exactly. as a direct output of the moving atlas formulation,

(iii) the quality of inversion for the biophysical parameters, i.e., tumor growth initial conditions and diffusion coefficient, in terms of correctness of the parameters compared to the known ground truth.

We report all performance measures in the patient space, i.e., with respect to the patient anatomy and not the atlas anatomy, since the patient space is the relevant space from an applications point of view. We use the following metrics: We measure the relative mismatch/residual between patient anatomy and reconstructed anatomy, between the reconstructed healthy patient anatomy and the ground truth, and between patient tumor and simulated tumor in the patient domain (compare Fig. 1):

[TABLE]

where $\boldsymbol{m}_{P}^{\star}(1)$ denotes the ground truth healthy patient brain. Based on the cardinality $|\,\cdot\,|$ of a set and a selection function $H$ with threshold $0.5$, i.e.,

[TABLE]

we calculate Dice coefficients for the individual label maps associated with the anatomy labels for $l\in\{WM,GM,CSF\}$ and their average across labels for the patient and the atlas anatomy

[TABLE]

Analoguously, we report values for the Dice coefficient computed for the healthy patient geometry, denoted by $\text{DICE}_{B_{0}}$, and for the anatomy label of the tumor, denoted by $\text{DICE}_{T}$. To monitor convergence of the Picard iterations, we report the relative value of an approximation to the gradient for the coupled problem (see §2.1) for the final iteration $k$:

[TABLE]

where $\boldsymbol{g}^{j}$ is an approximation to the gradient of the coupled optimization problem (3) (given by (7i) – (7k)) after the $i$th Picard iteration and $\boldsymbol{g}^{0}$ the corresponding approximate gradient for the initial guess. The gradient is only approximate, since we neglect the term $\partial_{\boldsymbol{m}_{P}(1)}\Phi(\boldsymbol{p})\alpha(0)$ in the final condition (7h) of the adjoint equation for the registration problem. This term is guaranteed to be small if the norm of the remaining gradient is small as the latter requires $\alpha(0)$ and the change in the estimated patient geometry to be small. Finally, we calculate the relative $\ell_{2}$-error for the initial condition:

[TABLE]

where $c(0)=\Phi\boldsymbol{p}$ and $\Phi\boldsymbol{p}^{\star}$ denotes the parametrization of the ground truth tumor initial condition.

4.4 Results for SYN – Synthetic Brain

This first, purely synthetic experiment uses two fairly different, very simplistic anatomies for the atlas and patient brain and an artificially grown tumor in one of them (called the patient brain).

Purpose.   This case is designed to cause the Moving Patient scheme to fail due to large differences in the brain anatomy between patient and atlas in the tumor region. We demonstrate that the Moving Atlas scheme retains good reconstruction quality for the tumor anatomy label, the tumor initial condition, and the characteristic diffusivity of tumor cells, while the Moving Patient scheme falls short as it inverts tumor growth in the ’wrong’ (atlas) geometry.

Setup.   For both, atlas and patient, we generate a purely synthetic brain anatomy, composed of an elliptical white matter inclusion surrounded by gray matter. The white matter inclusions for atlas and patient differ greatly in size. The target data $c_{D}$ are generated synthetically from a tumor-progression simulation (with $\boldsymbol{p}=\boldsymbol{p}^{\star}$) using our reaction-diffusion model with $\rho_{f}^{\star}=15$ and $k_{f}^{\star}=$110-2$$ from an initial condition with two nearby Gaussians enabled in the patient geometry. The healthy atlas geometry is depicted in the first row, the tumor bearing patient geometry in the last row of Fig. 2. We use the same set of $n_{\boldsymbol{p}}=343$ Gaussians with standard deviation $\sigma=\nicefrac{{2\pi}}{{30}}$ for the target data generation and inversion. For simplicity, tumor growth is enabled only in white matter. We invert for the registration velocity $\boldsymbol{v}$, the parameters $\boldsymbol{p}$ for the tumor initial condition, and for the characteristic diffusivity $k_{f}$ in white matter.

Observation.   In this case, the inclusions of white matter in atlas and patient vary significantly in size (cf. Fig. 2). Thus, registering the atlas to the patient requires strong local expansion of volume. Looking at Fig. 2 (right tableau), in the first iteration, to prevent the registration from directly matching the target tumor unaided by the tumor solver, we choose a higher regularization parameter $\beta_{\boldsymbol{v}}$ for the registration problem. Therefore, warping the target tumor to the atlas space basically corresponds to a copy operation. (before $k=2$). This warped-to-atlas tumor is used for tumor inversion in the atlas space. The subsequent registration step (with lower regularization) computes a velocity that registers the patient (small tumor in large white matter blob) to the atlas (small tumor in small white matter blob). This results in a compression of the patient white matter, and consequently yields a very small target tumor for the inversion in the atlas ($k=3$).151515Note, the displayed results are in the patient space, i.e., the inverse of the deformation is applied to display atlas labels in patient space, and the very small reconstructed tumor in the atlas is expanded when warped to the patient space. In subsequent iterations (with further reduced registration regularization), the warped-to-atlas target tumor is compressed further, ultimately causing the Moving Patient scheme to fail and returning a relative error of almost $100\%$ for the reconstruction of the tumor label, the reconstruction of the tumor initial condition, and the reconstruction of the characteristic diffusivity $k_{f}$ in white matter cf. Fig. 2). The new Moving Atlas scheme, however, maintains a good reconstruction quality and results in only $3.3\%$ relative errors for the tumor label, $18.4\%$ relative errors for the reconstruction of the tumor initial condition and $23.2\%$ relative error in the reconstruction of the characteristic diffusivity rate $k_{f}$ in white matter.

4.5 Results for STRV – Synthetic Tumor & Real Velocity

With the results for the STRV test cases, we show the general applicability of the new registration and tumor inversion coupling scheme as well as its improved properties in terms of correct reconstruction of initial tumor and healthy brain geometry. We evaluate the Moving Atlas scheme for a real brain geometry with synthetically generated target tumor such that the true tumor model parameters are known and we have access to the true healthy patient brain.

Purpose.   With the STRV test case experiments, we pursue three main purposes:

(i) We give experimental evidence for the convergence of the new Moving Atlas scheme for several settings;

(ii) we compare our coupled image registration and tumor inversion schemes Moving Patient and Moving Atlas to standalone tumor inversion, replacing the missing second snapshot (healthy patient) by a statistical atlas without involving image registration;

(iii) we show the superiority of the new Moving Atlas scheme over the previously presented Moving Patient scheme in terms of reconstruction quality of the (initially unknown) healthy patient brain anatomy and meaningful reconstruction of biophysical tumor-growth characteristics, in particular shape and sparsity of the ground truth tumor initial condition and prediction of the grown tumor.

Setup.   For both atlas and patient, we use real brain MR imaging data from two (healthy) individuals. The target tumor $c_{D}$ is generated synthetically from a tumor-progression simulation (with $\boldsymbol{p}=\boldsymbol{p}^{\star}$) using our reaction-diffusion model. Accordingly, the ground truth velocity field is unknown (or may not even exist due to a possibly wrong model), but the true tumor model parameters are known. We consider two different tumor locations (STRV-C1) and (STRV-C2) and varying time horizons for the tumor-progression simulation. An overview is given in Tab. 3. For (STRV-C1) we also employ different settings of the parametrization of the tumor initial condition: For parameter setting $\#1$, we use a set of $n_{\boldsymbol{p}}=125$ Gaussian basis functions with standard deviation $\sigma=\nicefrac{{2\pi}}{{30}}$, whereas, for settings $\#2$ and $\#3$, we use smaller (but more; $n_{\boldsymbol{p}}=343$) Gaussians with $\sigma=\nicefrac{{2\pi}}{{64}}$ to generate sparser ground truth initial conditions. The same set of Gaussian basis functions is used for target data generation and inversion. The target data are generated from a forward simulation with $\rho_{f}^{\star}=15$ and $k_{f}^{\star}=$110-1$$ and an initial condition with two nearby Gaussians enabled. For simplicity, tumor growth is enabled only in white matter. We invert for the tumor initial condition and the characteristic diffusivity $k_{f}$ in white matter.

Misfit and Gradient Reduction of the Moving Atlas Scheme.   To experimentally assess convergence of the Moving Atlas scheme over the Picard iterations, we report the $\ell_{2}$-mismatch and Dice overlay coefficients for the reconstruction of the (pathologic) brain anatomy ($\mu_{B,L^{2}}$ and $\text{DICE}_{B}$), the healthy patient anatomy ($\mu_{B_{0},L^{2}}$ and $\text{DICE}_{B_{0}}$), and the reconstructed (grown) tumor ($\mu_{T,L^{2}}$ and $\text{DICE}_{T}$) in Fig. 3 for the (STRV-C1) test case. Furthermore, we monitor the $\ell_{2}$-error $e_{c0,L^{2}}$ for the reconstructed tumor initial condition with respect to the ground truth and the relative norm $\|\boldsymbol{g}\|_{\text{rel}}$ of the (approximated) reduced gradient of the coupled formulation (3).

Observations.   Fig. 3 shows a monotonic reduction of data-misfit values and the corresponding gain in tissue overlay Dice scores for the brain anatomy, the healthy patient anatomy, the grown tumor, and the initial tumor condition. Although we do not have a proof for the convergence of the Moving Atlas Picard-iteration-type solution strategy (outlined in §2.2), we monitor the relative norm $\|\boldsymbol{g}\|_{\text{rel}}$ of the (approximated) reduced gradient for the Moving Atlas formulation (3), indicating convergence to an optimum. The norm of the approximated gradient continually decreases throughout our Picard-iteration-type solution strategy. Our solution scheme, and hence the gradient, neglect the term $\partial_{\boldsymbol{m}_{P}(1)}\Phi(\boldsymbol{p})\alpha(0)$ in (7h). This term, however, becomes very small as the misfit and, thus, $\alpha(0)$ decrease.

Standalone Tumor Inversion versus Moving Patient versus Moving Atlas.   To compare the different inversion schemes, we show results obtained using the newly introduced Moving Atlas (MA) strategy to the results for the Moving Patient (MP) solution scheme and those for a tumor standalone (T) solver without inter-subject registration, respectively, in Tab. 4 for both test cases (STRV-C1) and (STRV-C2).

Observations.   Solving the single-snapshot tumor inversion problem without inter-subject registration using a tumor standalone solver for some arbitrary normal brain anatomy (atlas) (cf. [10]) exhibits poor reconstruction quality. Our joint registration and biophysical inversion Moving Atlas approach (run $\#3$ and $\#6$ in Tab. 4) results in significantly higher reconstruction quality and outperforms the other strategies (Moving Patient, run $\#2$ and $\#5$, and tumor standalone, run $\#1$ and $\#4$) by a large margin, both in terms of data-misfit and Dice coefficient, as well as in reconstructing the true initial condition.

Comparing against the tumor standalone approach, we improve the reconstruction quality of the input data (grown tumor) from a relative data-misfit of $2.9010-1$ (run $\#1$) to $6.4910-2$ (run $\#3$) using the new Moving Atlas joint inversion scheme, and likewise improve the reconstruction quality for the tumor initial condition (compared to the ground truth) from a relative $\ell_{2}-$error of $4.2510-1$ for tumor standalone (run $\#1$) to $7.1310-2$ (run $\#3$) for the Moving Atlas scheme. Furthermore, the estimation of the characteristic diffusivity in white matter $k_{f}$ is significantly improved for the joint inversion scheme (going from a relative error of $4.8010-1$ (tumor standalone; run $\#4$) to a relative error of $7.7010-2$ (Moving Atlas; run $\#6$)).

Moving Patient versus Moving Atlas.   For a more detailed comparison of the Moving Atlas and the Moving Patient scheme, we show qualitative results of the reconstructed tumor initial condition and the grown tumor compared to the target data and ground truth initial condition in Fig. 4 (for STRV-C1) and Fig. 6 for (STRV-C2). In Fig. 7, we outline the initial and final residuals for the different brain tissue labels (white matter (WM), gray matter (GM), cerebrospinal fluid (CSF), grown tumor (TU($t=1$)), and tumor initial condition TU($t=0$); for axial slice 64 of a 3D-volume) before and after the joint registration and biophysical inversion for the test case variants (STRV-C1) and (STRV-C2). Furthermore, we display structural differences between the utilized normal brain anatomies for both inversion schemes. We show qualitative results for a longer time horizon of tumor evolution for the (STRV-C1) test case (with a very sparse initial condition) in Fig. 5.

We also study the sensitivity of the Moving Atlas and Moving Patient solution strategies with respect to perturbations in the tumor model parameters. In Tab. 5, we vary the value of the characteristic diffusivity in white matter from the ground truth, and monitor the obtained reconstruction performance.

Observations.   The new Moving Atlas scheme outperforms the Moving Patient counterpart [49] in various ways. Tab. 4 shows that, for the reconstruction of the brain anatomy (model prediction for anatomy labels of brain tissue labels WM, GM, and CSF), we reach a relative data-misfit ($\ell_{2}$-error) of $2.3310-1$ (run $\#3$) compared to $3.7210-1$ (run $\#2$) for the Moving Patient scheme (this translates to a Dice score of $9.4310-1$ (averaged over all labels) for the Moving Atlas compared to a Dice score of $8.7210-1$ for the old scheme). Similarly, for the approximation of the actual healthy patient anatomy, we improve from a Dice score of $8.7910-1$ (relative $\ell_{2}$-error of $25.6\%$; run $\#3$) obtained from the Moving Patient solution to a Dice score of $9.4710-1$ (relative $\ell_{2}$-error of only $16.2\%$; run $\#2$) using the new scheme. Fig. 7 illustrates the structural differences between the atlas and patient brain, the initial residuals and the residuals after joint inversion using the Moving Atlas and the Moving Patient scheme, respectively, for an exemplary slice (axial slice 64) of the 3D volume. We observe smaller errors for the Moving Atlas solution. The improved anatomy reconstruction is characteristic for the Moving Atlas strategy; we observe similar trends for various solver and parameter configurations (compare also Tab. 5 with inversion under perturbed or wrong tumor model parameters).

The Moving Atlas scheme furthermore results in improved reconstruction quality of the predicted grown tumor and higher similarity to the target data. For (STRV-C1), the Moving Atlas scheme results in a relative error of $6.4\%$ (Dice score of $9.5010-1$; run $\#3$) compared to a relative error of $14.4\%$ (Dice score of $9.2510-1$; run $\#2$) for the Moving Patient solution. The numbers for (STRV-C2) are very similar with $6.0\%$ relative error (Dice score of $9.7010-1$; run $\#6$) versus $19.0\%$ relative error (Dice score of $7.9610-1$; run $\#5$).

The primary objective for the derivation of the Moving Atlas scheme was to allow for a more informative inversion for biophysical model parameters. We primarily consider the estimation of the tumor initial condition, but also invert for the characteristic diffusivity $k_{f}$ in white matter. The relative errors161616Note, that for the Moving Atlas scheme, the reconstructed initial condition naturally “lives” in the patient space; for the Moving Patient solution strategy, we invert for the tumor initial condition in the atlas space. For a fair comparison, the reconstructed initial condition for the Moving Patient scheme is advected to the patient space before computing the relative error to the ground truth. $e_{c0,L^{2}}$ for the inversion of the ground truth tumor initial condition in Tab. 4 obtained for each scheme clearly attest the Moving Atlas scheme to be more reliable and sound in recovering biophysical parameters171717Our analysis assumes that tumor-growth is perfectly described by a reaction-diffusion model, which is quite certainly not true. When we say “biophysically meaningful” this is to be seen under the aforementioned assumption. Note furhermore, that the tumor model can easily be exchanged in the modular setting of our Picard iteration approach.. Using the improved Moving Atlas scheme, we achieve excellent reconstruction of the true initial condition with a relative error of only $7.1\%$ (for (STRV-C1); run $\#3$) and $6.3\%$ (for (STRV-C2); run $\#6$), as opposed to a relative error of $32.5\%$ (for (STRV-C1); run $\#2$) and $42.4\%$ (for (STRV-C2); run $\#5$), respectively, when using the Moving Patient scheme instead. An illustration of the error is given in Fig. 7 (for an axial cut at slice 64). We explain the improved inversion properties of the Moving Atlas scheme by its general idea to seek for a good approximation of the healthy patient brain anatomy first, and, thus, to carry out the inversion in the “correct” space. For the Moving Patient scheme, the inversion in the “wrong” anatomy (atlas space) has the potential to induce large errors for the estimation of biophysical parameters as shown in §4.4 since the inter-subject deformation map can aid matching the pathologic brains. In particular, there is no implication that this deformation map produces reasonable results when applied to the tumor initial condition (or other quantities that describe the time point of tumor genesis) in order to translate the latter into the individual patient application space. Visual inspection of our simulation results for longer time horizons in Fig. 5 support this point: The Moving Patient solution recovers a wrong position of the tumor initial condition (compare the green (ground truth) and purple (reconstruction) iso-volume representations of the initial tumor condition. A wrong position, wrong shape, or wrong sparsity of the initial condition can, in a second step, also cause a wrong estimation of the characteristic proliferation rate or tumor cell infiltration rate. For our experiments, the difference for the estimation of the characteristic diffusivity $k_{f}$ is, however, not significant; both schemes yield comparable results. One possible reason for this could be that the diffusion part of the gradient initially is large but then flattens out compared to the gradient component for the initial condition parametrization. The theoretical and numerical analysis of this effect is ongoing work.

5 Conclusion

We present a new inverse problem formulation to solve the patient-specific tumor inversion problem given a single data-snapshot in time only. We remedy the lack of healthy patient imaging data via registration to a healthy atlas brain, resulting in a joint inter-subject registration and biophysical inversion scheme. We simultaneously solve for an estimation of the inter-subject deformation map, and (a subset of) the biophysical model parameters for tumor progression. More specifically, we solve for an estimation of a sparse tumor initial condition, and for the characteristic diffusivity rate of tumor cells in white matter. The presented scheme conceptually improves on our earlier coupled solution approach [49]. The latter uses the registration to warp patient target data towards the atlas space, and inverts for tumor model parameters in the atlas brain (proxy for the unknown healthy patient brain). This may introduce large errors due to unwanted data fitting and model inversion in a wrong anatomy. We demonstrate this for synthetic cases. The new scheme computes an approximation of the healthy patient and performs the tumor inversion in the partient space. We furthermore encourage a sparse localization of the initial tumor (presented in our work [52]) to allow for a more reliable reconstruction of the model parameters. This biophysically motivated constraint can help to understand tumor genesis and its role in the subsequent tumor evolution. We derive a Picard-iteration-type solution strategy by dividing the strongly coupled set of first order optimality conditions into two solvers.

For our new method, we observe

(i) improved similarity of the reconstructed probability maps of brain tissue labels (WM, GM, CSF) with respect to the patient data, (ii) significantly smaller errors for the reconstruction of the grown tumor, (iii) a better approximation of the actual healthy patient brain anatomy, and, (iv) more reliable reconstruction of biophysical parameters, such as the tumor initial condition.

The reconstruction result obtained from the Moving Atlas scheme seems to be more “natural” as we invert in the “correct” brain anatomy of the (healthy) patient, and the registration does not act on the tumor probability map (as is the case for the Moving Patient strategy). In particular, using the Moving Atlas scheme, we were able to improve the Dice score for the brain anatomy from $8.7210-1$ for the Moving Patient strategy to $9.4310-1$; the Dice score for the approximation of the healthy patient brain anatomy from $8.7910-1$ for Moving Patient to $9.4710-1$; and the Dice coefficient for the tumor reconstruction from $8.3410-1$ for Moving Patient to $9.7010-1$. We observe an excellent reconstruction of the true initial condition with a relative error of only $7.1\%$ (STRV-C1) and $6.3\%$ (STRV-C2) (compared to $32.3\%$ and $42.4\%$ for Moving Patient) using the more sophisticated Moving Atlas scheme. We further conclude, that a sparse localization of the initial condition is essential if targeting biophysical parameter estimation from patient MRI. The previously used $\ell_{2}$-regularization results in rich initial conditions. Lastly, the Moving Atlas solution scheme seems to be more useful in identifying wrong parameters (such as cell proliferation rate and cell migration rate of brain tissue) than the Moving Patient counterpart.

Acknowledgements

This material is based upon work supported by a DAAD (German Academic Exchange Service) scholarship granted to the first author and the Tinsley Oden Fellowship, granted to the last author. The work is further supported by AFOSR grant FA9550-17-1-0190; by Simons Foundation grant (Award Number 586055, AM); by NSF grant CCF-1337393; by the U.S. Department of Energy, Office of Science, Office of Advanced Scientific Computing Research, Applied Mathematics program under Award Numbers DE-SC0010518 and DE-SC0009286; by NIH grant 10042242; by DARPA grant W911NF-115-2-0121; and by the University of Stuttgart, Institute for Parallel and Distributed Systems. 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, DARPA, and NSF. Computing time on the High-Performance Computing Centers (HLRS) Hazel Hen system (Stuttgart, Germany) was provided by an allocation of the federal project application ACID-44104.