Geometry-aware neural solver for fast Bayesian calibration of brain tumor models
Ivan Ezhov, Tudor Mot, Suprosanna Shit, Jana Lipkova, Johannes C. Paetzold, Florian Kofler, Fernando Navarro, Chantal Pellegrini, Marcel Kollovieh, Marie Metz, Benedikt Wiestler, Bjoern Menze
11 Real-time Bayesian personalization via a learnablebrain tumor growth model
Ivan Ezhov*, Tudor Mot*, Suprosanna Shit, Jana Lipkova, Johannes C. Paetzold, Florian Kofler, FernandoNavarro, Marie Metz, Benedikt Wiestler and Bjoern Menze
Abstract —Modeling of brain tumor dynamics has the potentialto advance therapeutic planning. Current modeling approachesresort to numerical solvers that simulate the tumor progressionaccording to a given differential equation. Using highly-efficientnumerical solvers, a single forward simulation takes up to afew minutes of compute. At the same time, clinical applicationsof the tumor modeling often imply solving an inverse problem,requiring up to tens of thousands forward model evaluationswhen used for a Bayesian model personalization via sampling.This results in a total inference time prohibitively expensivefor clinical translation. Moreover, while recent data-drivenapproaches become capable of emulating physics simulation,they tend to fail in generalizing over the variability of theboundary conditions imposed by the patient-specific anatomy.In this paper, we propose a learnable surrogate with anatomyencoder for simulating tumor growth which maps the biophysicalmodel parameters directly to simulation outputs, i.e. the localtumor cell densities. We test the surrogate on Bayesian tumormodel personalization for a cohort of glioma patients. Bayesianinference using the proposed neural surrogate yields estimatesanalogous to those obtained by solving the forward model witha regular numerical solver. The near real-time computation costrenders the proposed method suitable for clinical settings. Thecode is available at https://github.com/IvanEz/tumor-surrogate.
Index Terms —Bayesian inference, deep learning, glioma, modelpersonalization, tumor growth
I. I
NTRODUCTION S IMULATION of brain tumor progression can providecomplementary information to medical imaging for ra-diotherapy planning. As shown in [1]–[12], tumor modelingcan be employed to define a personalized radio-treatment areausing biophysical models to estimate the most likely directions * The authors contributed equally.I. Ezhov and S. Shit are supported by the Translational Brain ImagingTraining Network under the EU Marie Sklodowska-Curie programme (GrantID: 765148). B. Menze, B. Wiestler and F. Kofler are supported throughthe SFB 824, subproject B12, by DFG through TUM International GraduateSchool of Science and Engineering, GSC 81, and by the Institute for AdvancedStudy, funded by the German Excellence Initiative. Fernando Navarro issupported by DFG-GRK 2274.I. Ezhov, T. Mot, S. Shit, J. Paetzold, F. Navarro, F. Kofler are withthe Department of Informatics, and with TranslaTUM - Central Insti-tute for Translational Cancer Research, TUM, Munich, Germany (e-mail:[email protected]).B. Wiestler and M. Metz are with the Neuroradiology Department ofKlinikum Rechts der Isar, TUM, Munich, GermanyJ. Lipkova is with Harvard Medical School, Brigham and Women’s Hospi-tal, Boston, United StatesB. Menze is with the Department of Informatics, TranslaTUM, TUM,Munich, Germany, and Department of Quantitative Biomedicine of UZH,Zurich, SwitzerlandThis work has been submitted to the IEEE for possible publication.Copyright may be transferred without notice, after which this version mayno longer be accessible. of tumor cell infiltration instead of solely targeting tumor areavisible on a scan. These methodologies mainly imply solvingan inverse problem: finding the parameters of the biophysicaltumor growth model resulting in a simulation output that bestmatches an empirical observation outlining the pathology.Existing approaches for inverse tumor modeling resort todeterministic [2], [4], [13] as well as probabilistic Bayesian[14]–[16] formalisms. All the approaches require multipleforward simulations to ensure convergence of the parametricestimation. The number of simulations ranges from severalthousand for approximate methods [5], [14], [16] to tens ofthousands in case of fully Bayesian analysis [3]. The forwardbrain tumor models are often based on the reaction-diffusionequation and are implemented using highly-efficient numericalsolvers. In [15], authors employ the Lattice Boltzmann methodwhich allows parallelized computing and takes ca. 80 secondson a 60 core machine for simulating the pathology growth. In[3], the forward model is implemented by means of a multi-resolution adapted grid solver with a simulation time of 1-3minutes using 2 cores. Despite the computational advances ofthe solvers, the minutes of a single forward model evaluationmultiplied by thousands of forward integrations necessary forthe inverse problem can result in weeks of total computingtime. This constrains the testing of more elaborate tumormodels (e.g., considering cell mixtures or multiple competingpatho-physiological processes [17]), and translation of thepersonalized radiotherapy planning into clinical practice [1]–[4].As recent years showed, speeding up heavy conventionalcomputation becomes feasible using end-to-end learning meth-ods. The data-driven methodology has also penetrated the fieldof numerical computing. Learnable surrogates were proposedfor various scientific computing tasks in the natural sciencesby exploiting fully-connected [18]–[20], convolutional [21]–[23], and hybrid [24]–[26] neural architectures. Among themare two methods that proved capable of learning even a directmapping from the space of parameters driving a simulator tothe space of the simulator solutions in a static geometry [22],[23]. Unfortunately, these promising methods are incapable ofdealing with inference in arbitrary geometries. This limits theirtransfer to model personalization that is crucially dependant onan adaption to the patient specific simulation domain.The contribution of the paper is the following: we introducea learnable method emulating a numerical tumor (glioma)growth forward solver. To the best of our knowledge, this is thefirst network-based approach in the computational pathologyfield that maps parameters of the biophysical model directly to a r X i v : . [ c s . C E ] S e p x2 Conv block D ⍴ T Conv block Conv block Conv k3s1m1x2 x28 x x x ReLu Convk3s1m128N=4
Conv blockx2
Conv block
Conv block Conv blockx2 x2 Conv block Conv k3s1m1288 x x x U pred Brain anatomy encodingBrain tumor decoding
WMCSFGM x x x Fig. 1:
Learnable brain tumor model surrogate. The network is composed of two main parts: a) brain anatomy encoder that maps theanatomy volumes (WM, GM, CSF) to a latent representation, b) brain tumor decoder that takes as an input a 1D vector of the parameters { D, ρ, T } , concatenated with the latent representation from (a), and maps the resulting tensor to the 3D tumor simulation volume. Theconvolutional block is composed of N=4 repetitions of convolutional operation (with kernel size k=3, stride s=1, number of filters m=128)and ReLu non-linearity. the simulation outputs while generalizing over the simulationgeometry. We achieve a 400 × speed-up comparing to anadvanced numerical solver by employing the tumor modelsurrogate with an anatomy encoder that enforces patient-specific boundary conditions. This enables a fast Bayesianmodel personalization that is consistent with the baselinenumerical solver. II. M ETHOD
1) Forward tumor model:
The simulations that we aimto emulate are generated by a 3D numerical solver relyingon a special partial-differential equation (PDE), the Fisher-Kolmogorov equation. The equation describes the evolutionof the pathology by considering diffusion and proliferation ofthe tumor cells under the Neumann boundary condition (B.C.): ∂ u ∂t = ∇ ( D ∇ u ) + ρ u (1 − u ) , (1) ∇ u · n = 0 B.C. (2)Here, u is the normalized 3D tumor cell density, D denotesthe diffusion tensor, ρ is the rate of cell proliferation, and WMCSF GMT1c
Fig. 2:
An example of an MRI T1c scan from the dataset andcorresponding segmentations obtained by registering the brain atlasto patient space. n is the normal vector to the boundary. We assume themigration of the tumor cells to occur only in white matter(WM) and grey matter (GM), considering isotropic diffusionwith D = D · I (where D ∈ { D w , D g } is a diffusioncoefficient in white or grey matter, and I is an identity matrix).Both WM and GM constitute the simulation domain whilecerebrospinal fluid (CSF) and skull determine the patient-specific boundary. The input to the solver is a set of parameters θ P = { D w , ρ, x, y, z, T } , where x, y, z define the position of the function u at time T = 0 , which is initialized as a pointsource.We simulate the tumors in multiple brain anatomies, whichare extracted from medical scans of patients diagnosed withthe tumor.
2) Learnable forward model surrogate with anatomy en-coder:
Our goal is to learn a surrogate that could map param-eters of the pathology model θ P to corresponding simulations u ( θ P ) . For this, we base our method on [23] which is designedto do the mapping for fluid simulations in a static spatial do-main. Different from [23], we need to consider patient-specificboundary conditions. To this end, we introduce an anatomyencoder that imposes anatomical boundary conditions, Fig. 1.The numerical solver’s output u ( θ P ) has a size of × × voxels. However, to provide a greater anatomicalvariability to the dataset on which we train the surrogate, wecrop all simulated outputs and corresponding brain anatomiesto × × volumes, centered at the initialization location x, y, z . The × × is greater than half the brain sizeand tumors bigger than this are incompatible with life.The architecture of the tumor model surrogate consists oftwo main parts:- Brain anatomy encoder which encodes 3D volumes ofthe brain tissues WM, GM, and CSF through a series ofconvolutional blocks. The blocks are composed of alternatingconvolution operations (with fixed parameters of kernel size3, stride 1, and the number of maps 128) and a non-linearityin the form of a linear rectifier. Each block is equipped witha skip connection linking the input and output of the blockvia an element-wise sum. Downsampling between the blocksis achieved by a convolutional operation with a stride of two.-
Brain tumor decoder that takes as an input a 1D vector ofthe parameters { D, ρ, T } alongside with the encoded anatomy.Note that we do not condition the decoder on the initializationlocation x, y, z , since as mentioned above we crop all trainingvolumes exactly at this location. Thus the network is taught toreproduce the tumor in the center for any volume. Before beingpassed to the decoder, the 1D vector of model parametersis mapped via a fully connected layer to a tensor of size × × × and is concatenated with the tensor of theencoded brain anatomy. The resulting tensor is graduallyupsampled through a series of convolutional blocks analogousto the encoder and nearest-neighbor upsampling (we refer thereader to [27] for a detailed discussion on why such type ofupsampling is preferred over the deconvolution operation). Atthe end of the series, a 4D tensor × × × isconvolved to the output prediction - a 3D tumor simulationvolume (height 64, width 64, depth 64). The decoder designis adopted from [23], [28].As a loss function, we compute the error between thepredicted ( u pred ) and simulated ( u sim ) cell concentrationfunctions under L norm separately in the CSF area and thearea of non-zero ground truth simulated tumor (GT): L total = (cid:107) u sim − u pred (cid:107) GT > + (cid:107) u sim − u pred (cid:107) CSF (3)As shown in the ablation analysis (Fig. 11), such lossdefinition significantly affects predictions performance in new patient geometry.
3) Bayesian model personalization:
To demonstrate theapplicability of the neural surrogate, we perform Bayesiantumor growth model personalization substituting the numericalsolver with the learnable one. As calibration data we use twotypes of imaging modalities: (a) T1 contrast-enhanced andFLAIR MRI modalities that allow estimating the morphologi-cal characteristic of the visible tumor; (b) FET-PET scans thatprovide information about tumor metabolic activity.Analogous to [14], [15], we relate the output of the tumorgrowth solver u ( θ P ) to imaging information via a probabilisticmodel, p ( D | u , θ ) = p (cid:0) y T c | u , θ P (cid:1) · p (cid:0) y F LAIR | u , θ I (cid:1) ·· p (cid:0) y P ET | u , θ I (cid:1) (4)where image observations D = { y T c , y F LAIR , y P ET } areassumed to be independent, and θ = { θ P , θ I } constituteparameters of the pathology model θ P and the probabilisticimaging model θ I . The latter is defined differently accordingto the type of modality: - MRI modalities provide information in the form of binarytumor segmentations ( y T c , y F LAIR ). Thus we assign for eachvoxel a discrete label y Mi ∈ { , } , M ∈ { T c, F LAIR } . Wemodel the probability of observing y T c , y F LAIR for a givenconcentration u ( θ P ) with Bernoulli distribution, p (cid:0) y T c,F LAIR | u , θ MI (cid:1) = (cid:89) i p (cid:0) y Mi | u i , θ MI (cid:1) == (cid:89) i α y Mi i,M · (1 − α i,M ) − y Mi , (5)where α i,M is the probability of observing tumor-inducedchanges defined as a double logistic sigmoid, α i,M ( u i , u c ) = 0 . . · sign (cid:0) u i − u Mc (cid:1) − e − ( ui − uMc ) σ α (6)With this formulation, we postulate that the tumor is not visibleon a scan below the threshold level u Mc . The parameter σ α isintroduced to take uncertainty in the threshold concentration u Mc into account. - FET-PET modality ( y P ET ) provides continuous informa-tion in each voxel and can be assumed to be proportional to
T1G d FLAIR Fet-pet
Fig. 3: An illustration of the imaging information used forBayesian model calibration: binary segmentations obtainedfrom T1Gd and FLAIR modalities, and FET-PET signal whichis proportional to the tumor density.
CSF
Fig. 4: Qualitative convergence plot for the network predic-tions over training epochs. For all columns the network inputparameters { D, ρ } and the brain anatomy are fixed, whereaseach column has a different value of T from 100 to 1000days with intervals of 100 days. The images were obtainedby inferring individual 3D volumes (in a single batch) for alltime values and taking a central 2D slice from each volume.The last rows framed in green correspond to the groundtruth simulation. The CSF delineation constraining the tumorgrowth is framed in red.Fig. 5: Convergence plot of the total loss. The total loss L total evaluated on a random batch of size 16 from the test set. Weused samples from 10 patients for training and samples fromthe remaining 5 patients for test.the tumor density [3]. In this case we model the likelihoodimaging function by a Gaussian distribution with unknownconstant of proportionality b : p (cid:0) y P ET | u , θ MI (cid:1) = (cid:89) i p (cid:0) y i | u i , θ MI (cid:1) = (cid:89) i N (cid:0) y i − bu i , σ (cid:1) (7)Here, σ models uncertainty in the PET signal, which isconsidered to be normalized y P ETi ∈ [0 , .In total there are eleven parameters (six pathology modelparameters θ P = { D w , ρ, x, y, z, T } and five imaging modelparameters θ I = { u T cc , u F LAIRc , σ α , b, σ } ) which we inferfrom the triplet of medical scans D = { y T c , y F LAIR , y P ET } using a Markov Chain Monte Carlo (MCMC) sampling algo-rithm [29].
4) Implementation:
The numerical tumor solver used forobtaining the simulation dataset is a highly parallelized gliomasolver returning a 3D tumor volume on a uniform spatial grid[3], [30].The surrogate network is trained using the Adam optimizer[31] with decay rates β = 0 . and β = 0 . for 30epochs, which was observed to be sufficient for convergence.The learning rate is cosine annealed from − to − over the training and the batch size is 16. We performedthe experiments on an NVIDIA Quadro P8000 using theTensorflow framework.For Bayesian MCMC inference we use an implementationof Transitional MCMC from [32] with 2048 samples periteration. III. E XPERIMENTS
1) Data:
The patient-specific tumor-free brain anatomies(WM, GM, CSF) in which we model the tumor were ob-tained from an MRI dataset of 15 glioblastoma patients byregistration-based approach proposed in [3], [33]. The datahave resolution 256x256x256 with isotropic voxels of 1mmside-length, but for simulation purposes, the calibration is per-formed in data downsampled to 128x128x128. The diffusioncoefficient in white matter D w is considered to be 10 timesgreater than in grey matter D g . The following ranges areused for random uniform sampling of the model parameters: D w ∈ [0 . , . mm /day , ρ ∈ [0 . , .
03] 1 /day , and T ∈ [50 , day with a step size of 50 days. The tumorlocation coordinates x, y , and z are sampled within the brainvolumes. Samples that have initial location { x, y, z } within theCSF do not trigger any tumor growth and were thus discarded.For samples with initial location closer than 32 pixels to anyof the volume borders, we do extra padding when croppingto 64x64x64 size. In total, we have a set of 20k parameters-simulation pairs for training.
2) Training the surrogate:
Conventional application of tu-mor modeling such as model personalization (via solving aninverse problem) implies estimating the model parameters bysampling over fixed, physiologically plausible ranges. Thus,we aim to employ the learnable surrogate which i) possessesan interpolation capacity for the parameters { D, ρ, T } , and ii)is capable to extrapolate for simulations in new geometries. Toprobe these properties, the test set of 2k samples was formedto have only the brain anatomies unseen by the network duringtraining, while the parametric { D, ρ, T } triplets were sampledfrom the same ranges as for the training. a) Quantitative convergence: Fig. 5 shows the conver-gence of the loss L total for both training and test sets. Fig. 6demonstrates the distribution of the mean absolute error withineach class of the brain tissues evaluated on the whole test set.Even though we observe samples with the error in the orderof − , the majority of the distribution lies within the orderof − . In Fig. 7 we depict the histograms of the DICEscore computed between the simulated and predicted tumorvolumes which are thresholded at different levels of the tumor (a) WM (b) GM (c) CSF Fig. 6: Histograms of the mean absolute error (cid:107) u pred − u sim (cid:107) computed on the whole validation set within each class of thebrain tissues (WM, GM, CSF). (a) DICE > > > Fig. 7: Histograms of the DICE score computed on the whole validation set for the tumor volumes u pred and u sim thresholdedat 0.2, 0.4, and 0.8 values of tumor cell concentration. The means of the distributions are: a) 0.821, b) 0.813, c) 0.815.cell concentration (such thresholding is exactly the operationthat we perform during the Bayesian inference to relate tumorconcentration with MRI signals, Eq. 6). The distributions arecentered close to the DICE 1.0. The smaller peak at DICE 0.0is due to the fact that thresholding of some simulations resultsin volumes containing a few (up to 5) non-zero voxels, whereasthresholding of the corresponding network predictions outputsvolumes of all-zero voxels (or vice versa). In our particularapplication of the tumor model personalization such volumesare significantly smaller that the binary MRI segmentationvolumes to which we calibrate the model and thus do notaffect the Bayesian inference outcome. b) Qualitative convergence: Fig. 4 illustrates visually thesurrogate’s predictions for simulations from the test set. Thepredicted distribution of the tumor concentration accuratelyemulates the ground truth simulations. The global constraintsimposed by the CSF anatomy (in which the tumor is con-strained to grow) are well captured by the network. Fig. 8depicts a comparison between a network prediction and asimulated ground truth in 3D.
3) Performance speed-up:
The computing time for a singlesimulation using the numerical solver is on average 120seconds using an Intel Xeon with 8 CPU cores and 64GBof RAM. The time required by using the neural solver forprocessing a batch of size 8 during inference is equal to2.4s, i.e. our surrogate is 400 times faster than the numericalsolver per single simulation. Such a speed-up reduces thecomputing time for applications such as model personalizationfrom weeks to minutes. Fig. 8: 3D comparison between the predicted (right) andground truth (left) tumor simulations embedded in the brainanatomy.We want to point out that a theoretical comparison betweenCPU- and GPU-based computations should be taken withcare for a few reasons. First, the neural surrogate does notneed to do any sequential computations and can batch processmultiple frames, whereas the numerical solver is successivelyintegrating the PDE equation in time and for a single sample.Second, the CPU-based solver acts on a grid of dimension128x128x128 whereas the network predicts the tumor ina cropped region of interest of 64x64x64 voxels. Finally,the convolutional layers are very efficiently realized in mostGPU-tailored deep learning frameworks. In this regard, a faircomparison could imply to contrast the GPU compute with aCPU cluster containing as many processing units as the GPU.However, the latter is rarely available for practical applications,and as our main goal is clinical translation, we compare CPU
P1 P2 P3 P4 P5 N e u r a l N u m e r i c a l P E T - F E T F L A I R T c Fig. 9: Results of the Bayesian inference for patients P1-P5. The three upper rows correspond to the imaging modalities usedfor tumor model calibration, namely T1c, FLAIR, and PET-FET. The two bottom rows show the simulations of glioma withmodel parameters inferred via the Bayesian inference using the numerical solver and proposed neural surrogate. and GPU platforms that are widely available. Also, we empha-size that when both solvers are tested on the application suchas Bayesian personalization, we treat the empirical imagingdata as well as the simulated domain identically keeping thesize and resolution same (even though the network predicts64x64x64 tumor volumes centered around the initial location { x, y, z } , they are then embedded in the original 128x128x128simulation domain for the personalization task).
4) Bayesian model personalization in patient data:
Weperformed the Bayesian brain tumor model calibration onpreoperative scans of 5 test patients using the proposed neuralsurrogate and numerical solver (the max-a-posterior (MAP)estimates of the tumor density are provided in Tab.1 inthe appendix). Tumor concentrations modeled with the MAPestimates are shown in Fig. 9. We observe an agreement ofthe glioma profiles obtained by the two methods and thevariability of the estimations is within the variability of theBayesian calibration. We also note that the learnable surrogateis trained on a dataset with continuous uniform distributionof the diffusion coefficient D and the proliferation rate ρ ,whereas the T parameter has a discrete interval. However,during the model calibration, the time parameter is sampledfrom a continuous interval. This implicitly suggests that thenetwork’s interpolation capacity is sufficient for parametricinference.Lastly, even though we tested the surrogate on the simplisticFisher-Kolmogorov model, a translation of the method to morecomplicated tumor growth descriptions should be straightfor-ward (through conditioning with additional model parameters)at no extra inference time cost.IV. D ISCUSSION
Despite the mentioned advantages, we would like to addressthe limitations of the proposed method. As depicted in Fig.7, even though the Dice distributions are peaked close to aperfect score, there is a room for improvement to reduce thedeviation. As seen in Fig. 4, the anatomy encoder capturesthe global constraints imposed by the CSF anatomy (in whichthe tumor is constrained to grow). However, conditioning ofthe network on the WM and GM maps improves moderatelythe performance as compared to training without such con-ditioning (see Fig. 7 and 10). This is somewhat expected asthe WM and GM geometries affect numerical simulations ina highly non-linear fashion as compared to the CSF, geometryof which is preserved pixel-wise in the output tumor volume.We experimented with several architectural and hyper-parametric configurations and selected the one providing thelowest error on the validation set (Fig. 10 and 11 show exam-ples of the ablation analysis). Given that the space of networkdesigns is close to infinite, future research could benefit fromusing automated methods for identifying an optimal design,e.g. exploiting neural architecture search analogous to [22].We do not have a clear answer to whether the mentionedlimitations and resulting approximation errors are acceptablefor clinical translation. As well as we do not know whetherthe Bayesian calibration itself under the simplistic Fisher-Kolmogorov formalism is suitable for the translation (both require a study on a large cohort of patients, post-surgeryanalysis, etc. [1]–[4]). However, we find the proposed methodto be a solid baseline in the search for optimal tumor modelsurrogate, which in turn can significantly speed up our searchfor a biophysical model descriptive enough for clinical trials.V. C
ONCLUSION
We present the first learnable surrogate with anatomy en-coder for tumor growth modeling that is capable of map-ping the model parameters to the corresponding simula-tions while respecting patient-specific anatomy. Our methodachieves close to real-time simulation allowing fast inferencein Bayesian settings. The technique can be readily adopted tomore complicated tumor growth models and similar 4D inversemodeling tasks. R
EFERENCES[1] J. Unkelbach and et al., “Radiotherapy planning for glioblastoma basedon a tumor growth model: improving target volume delineation,”
Physicsin Medicine and Biology , vol. 59, no. 3, pp. 747–770, 2014.[2] R. Rockne and et al., “A patient-specific computational model ofhypoxia-modulated radiation resistance in glioblastoma using 18f-fmiso-pet,”
JRS, Interface , vol. 12, 2015.[3] J. Lipkova and et al., “Personalized radiotherapy design for glioblastoma:Integrating mathematical tumor models, multimodal scans and bayesianinference,”
IEEE Transactions on Medical Imaging , pp. 1–1, 2019.[4] E. Konukoglu, O. Clatz, H. Delingette, and N. Ayache, “Personalizationof reaction-diffusion tumor growth models in mr images: Application tobrain gliomas characterization and radiotherapy planning,”
MultiscaleCancer Modeling , 12 2010.[5] M. Le and et al., “Personalized radiotherapy planning based on acomputational tumor growth model,”
IEEE Transactions on MedicalImaging , vol. 36, no. 3, pp. 815–825, mar 2017.[6] S. Subramanian, K. Scheufele, N. Himthani, and G. Biros, “Multiatlascalibration of biophysical brain tumor growth models with mass effect,” arXiv preprint arXiv:2006.09932 , 2020.[7] K. Scheufele, S. Subramanian, and G. Biros, “Automatic mri-drivenmodel calibration for advanced brain tumor progression analysis,” arXiv:Medical Physics , 2020.[8] K. Scheufele, S. Subramanian, A. Mang, G. Biros, and M. Mehl, “Image-driven biophysical tumor growth model calibration,” arXiv: QuantitativeMethods , 2019.[9] L. Zhang, L. Lu, R. M. Summers, E. Kebebew, and J. Yao, “Convolu-tional invasion and expansion networks for tumor growth prediction,”
IEEE transactions on medical imaging , vol. 37, no. 2, pp. 638–648,2017.[10] J. Petersen and et al., “Deep probabilistic modeling of glioma growth,” in
International Conference on Medical Image Computing and Computer-Assisted Intervention . Springer, 2019, pp. 806–814.[11] A. Elazab and et al., “Macroscopic cerebral tumor growth modeling frommedical images: A review,”
IEEE Access , vol. 6, pp. 30 663–30 679,2018.[12] A. Mang, S. Bakas, S. Subramanian, C. Davatzikos, and G. Biros,“Integrated biophysical modeling and image analysis: Application toneuro-oncology,”
Annual Review of Biomedical Engineering , vol. 22,pp. 309–341, 2020.[13] P. R. Jackson, J. Juliano, A. Hawkins-Daarud, R. C. Rockne, andK. R. Swanson, “Patient-specific mathematical neuro-oncology: Usinga simple proliferation and invasion tumor model to inform clinicalpractice,”
Bulletin of Mathematical Biology , vol. 77, no. 5, pp. 846–856, mar 2015.[14] B. Menze and et al., “A generative approach for image-based modelingof tumor growth,” in
IPMI , 2011, pp. 735–747.[15] M. Le and et al., “Bayesian personalization of brain tumor growthmodel,” in
MICCAI , 2015, pp. 424–432.[16] I. Ezhov and et al., “Neural parameters estimation for brain tumor growthmodeling,” in
MICCAI , 2019, pp. 787–795.[17] V. Cristini and J. Lowengrub,
Multiscale modeling of cancer: an inte-grated experimental and mathematical modeling approach . CambridgeUniversity Press, 2010. [18] M. Raissi, P. Perdikaris, and G. E. Karniadakis, “Physics-informedneural networks: A deep learning framework for solving forward andinverse problems involving nonlinear partial differential equations,”
Journal of Computational Physics , vol. 378, pp. 686–707, 2019.[19] V. Sitzmann and et al., “Implicit neural representations with periodicactivation functions,” arXiv preprint arXiv:2006.09661 , 2020.[20] B. Stevens and T. Colonius, “Finitenet: A fully convolutional lstmnetwork architecture for time-dependent partial differential equations,” arXiv preprint arXiv:2002.03014 , 2020.[21] N. Thuerey, K. Weißenow, L. Prantl, and X. Hu, “Deep learning methodsfor reynolds-averaged navier–stokes simulations of airfoil flows,”
AIAAJournal , vol. 58, no. 1, pp. 15–26, 2020.[22] M. Kasim and et al., “Up to two billion times acceleration of scien-tific simulations with deep neural architecture search,” arXiv preprintarXiv:2001.08055 , 2020.[23] B. Kim and et al., “Deep fluids: A generative network for parameterizedfluid simulations,”
Comput. Graph. Forum , vol. 38, pp. 59–70, 2018.[24] J. Hsieh, S. Zhao, S. Eismann, L. Mirabella, and S. Ermon, “Learn-ing neural pde solvers with convergence guarantees,”
ArXiv , vol.abs/1906.01200, 2019.[25] S. Shit and et al., “Implicit neural solver for time-dependent linear pdeswith convergence guarantee,” arXiv preprint arXiv:1910.03452 , 2019.[26] A. Sanchez-Gonzalez and et al., “Learning to simulate complex physicswith graph networks,” arXiv preprint arXiv:2002.09405 , 2020.[27] A. Odena, V. Dumoulin, and C. Olah, “Deconvolution and checkerboardartifacts,”
Distill , 2016. [Online]. Available: http://distill.pub/2016/deconv-checkerboard[28] D. Berthelot, T. Schumm, and L. Metz, “Began: Boundary equilib-rium generative adversarial networks,” arXiv preprint arXiv:1703.10717 ,2017.[29] J. Ching and Y.-C. Chen, “Transitional markov chain monte carlomethod for bayesian model updating, model class selection, and modelaveraging,”
Journal of engineering mechanics , vol. 133, no. 7, pp. 816–832, 2007.[30] J. Lipkova, https://github.com/JanaLipkova/GliomaSolver.[31] D. P. Kingma and J. Ba, “Adam: A method for stochastic optimization,”in
ICLR 2015 , 2015.[32] P. Hadjidoukas and et al., “ π Journal of Computational Physics , vol. 284, pp. 1–21, 2015.[33] F. Kofler and et al., “Brats toolkit: Translating brats brain tumorsegmentation algorithms into clinical and scientific practice,”
Frontiersin Neuroscience , vol. 14, 2020.
D ρ T x y z σ b u T cc u FLAIRc σ α P1 Numericalsolver 0.08994 0.01550 572.582 0.3098 0.6748 0.2819 0.2177 0.6384 0.7976 0.3718 0.0514Neuralsurrogate 0.08976 0.01657 553.019 0.3071 0.6657 0.2818 0.1702 0.7198 0.7344 0.3364 0.0500P2 Numericalsolver 0.08920 0.01993 503.146 0.6121 0.6706 0.3411 0.2358 0.7309 0.6502 0.4376 0.0688Neuralsurrogate 0.08945 0.02882 418.631 0.6217 0.6445 0.3316 0.2434 0.6324 0.6316 0.4337 0.0799P3 Numericalsolver 0.08998 0.01555 572.128 0.4129 0.3544 0.2900 0.2326 0.7023 0.6088 0.3608 0.0735Neuralsurrogate 0.08990 0.02897 418.946 0.4230 0.3678 0.2901 0.2498 0.6561 0.6294 0.3662 0.0772P4 Numericalsolver 0.08988 0.01372 609.455 0.3293 0.5657 0.2701 0.2496 0.6894 0.6010 0.3802 0.0758Neuralsurrogate 0.08920 0.02085 493.141 0.3578 0.5839 0.2761 0.2469 0.6287 0.6009 0.3358 0.0679P5 Numericalsolver 0.08938 0.01080 686.923 0.5530 0.6203 0.3147 0.2484 0.6000 0.6023 0.3775 0.0796Neuralsurrogate 0.08360 0.01066 687.841 0.5701 0.6259 0.3276 0.2411 0.7044 0.6098 0.3790 0.0759
TABLE I: MAP estimates of the marginal distribution from the Bayesian calibration on the patient data using the numericalsolver and proposed neural surrogate. The prior ranges are chosen as follows: D w ∈ [0 . , . mm /day , ρ ∈ [0 . , . /day , T ∈ [30 , days , σ ∈ [0 . , . , b ∈ [0 . , . , u T cc ∈ [0 . , . , u F LAIRc ∈ [0 . , . , and σ α ∈ [0 . , . . (a) DICE > > > Fig. 10: Ablation analysis I. Instead of inputting 3 volumes of different tissue types (WM, GM, CSF) only a single volumeof CSF tissue serves as an input to the network. The figure shows the histograms of the DICE score computed on the wholevalidation set for the tumor volumes u pred and u sim thresholded at 0.2, 0.4, and 0.8 values of tumor cell concentration. Themeans of the distributions are: a) 0.798, b) 0.799, c) 0.801. (a) DICE > > > Fig. 11: Ablation analysis II. Instead of splitting the L total into two terms as in Eq. 3, a standard L1 loss is used for computingthe error in the whole volume. The figure shows the histograms of the DICE score computed on the whole validation set forthe tumor volumes u pred and u simsim