{'fraction_non_alphanumeric': 0.032071226844095545, 'fraction_numerical': 0.015413335886266814, 'mean_word_length': 4.769676884838442, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Real-time Bayesian personalization via a learnable brain tumor growth model

Modeling of brain tumor dynamics has the potential to advance therapeutic planning. Current modeling approaches resort to numerical solvers that simulate the tumor progression according to a given differential equation. Using highly-efficient numerical solvers, a single forward simulation takes up to a few minutes of compute. At the same time, clinical applications of the tumor modeling often imply solving an inverse problem, requiring up to tens of thousands forward model evaluations when used for a Bayesian model personalization via sampling. This results in a total inference time prohibitively expensive for clinical translation. Moreover, while recent data-driven approaches become capable of emulating physics simulation, they tend to fail in generalizing over the variability of the boundary conditions imposed by the patient-specific anatomy. In this paper, we propose a learnable surrogate with anatomy encoder for simulating tumor growth which maps the biophysical model parameters directly to simulation outputs, i.e. the local tumor cell densities. We test the surrogate on Bayesian tumor model personalization for a cohort of glioma patients. Bayesian inference using the proposed neural surrogate yields estimates analogous to those obtained by solving the forward model with a regular numerical solver. The near real-time computation cost renders the proposed method suitable for clinical settings. The code is available at https://github.com/IvanEz/tumor-surrogate.

of tumor cell infiltration instead of solely targeting tumor area visible on a scan. These methodologies mainly imply solving an inverse problem: finding the parameters of the biophysical tumor growth model resulting in a simulation output that best matches an empirical observation outlining the pathology.
Existing approaches for inverse tumor modeling resort to deterministic [2], [4], [13] as well as probabilistic Bayesian [14]- [16] formalisms. All the approaches require multiple forward simulations to ensure convergence of the parametric estimation. The number of simulations ranges from several thousand for approximate methods [5], [14], [16] to tens of thousands in case of fully Bayesian analysis [3]. The forward brain tumor models are often based on the reaction-diffusion equation and are implemented using highly-efficient numerical solvers. In [15], authors employ the Lattice Boltzmann method which allows parallelized computing and takes ca. 80 seconds on a 60 core machine for simulating the pathology growth. In [3], the forward model is implemented by means of a multiresolution adapted grid solver with a simulation time of 1-3 minutes using 2 cores. Despite the computational advances of the solvers, the minutes of a single forward model evaluation multiplied by thousands of forward integrations necessary for the inverse problem can result in weeks of total computing time. This constrains the testing of more elaborate tumor models (e.g., considering cell mixtures or multiple competing patho-physiological processes [17]), and translation of the personalized radiotherapy planning into clinical practice [1]- [4].
As recent years showed, speeding up heavy conventional computation becomes feasible using end-to-end learning methods. The data-driven methodology has also penetrated the field of numerical computing. Learnable surrogates were proposed for various scientific computing tasks in the natural sciences by exploiting fully-connected [18]- [20], convolutional [21]- [23], and hybrid [24]- [26] neural architectures. Among them are two methods that proved capable of learning even a direct mapping from the space of parameters driving a simulator to the space of the simulator solutions in a static geometry [22], [23]. Unfortunately, these promising methods are incapable of dealing with inference in arbitrary geometries. This limits their transfer to model personalization that is crucially dependant on an adaption to the patient specific simulation domain.
The contribution of the paper is the following: we introduce a learnable method emulating a numerical tumor (glioma) growth forward solver. To the best of our knowledge, this is the first network-based approach in the computational pathology field that maps parameters of the biophysical model directly to  Fig. 1: Learnable brain tumor model surrogate. The network is composed of two main parts: a) brain anatomy encoder that maps the anatomy 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. The convolutional 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 simulation geometry. We achieve a 400× speed-up comparing to an advanced numerical solver by employing the tumor model surrogate with an anatomy encoder that enforces patientspecific boundary conditions. This enables a fast Bayesian model personalization that is consistent with the baseline numerical solver.
II. METHOD 1) Forward tumor model: The simulations that we aim to emulate are generated by a 3D numerical solver relying on a special partial-differential equation ( ∇u · n = 0 B.C. Here, u is the normalized 3D tumor cell density, D denotes the diffusion tensor, ρ is the rate of cell proliferation, and WM CSF GM T1c n is the normal vector to the boundary. We assume the migration of the tumor cells to occur only in white matter (WM) and grey matter (GM), considering isotropic diffusion with D = D · I (where D ∈ {D w , D g } is a diffusion coefficient in white or grey matter, and I is an identity matrix). Both WM and GM constitute the simulation domain while cerebrospinal fluid (CSF) and skull determine the patientspecific 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 point source.
We simulate the tumors in multiple brain anatomies, which are extracted from medical scans of patients diagnosed with the tumor.
2) Learnable forward model surrogate with anatomy encoder: Our goal is to learn a surrogate that could map parameters of the pathology model θ P to corresponding simulations u(θ P ). For this, we base our method on [23] which is designed to do the mapping for fluid simulations in a static spatial domain. Different from [23], we need to consider patient-specific boundary conditions. To this end, we introduce an anatomy encoder that imposes anatomical boundary conditions, Fig. 1.
The numerical solver's output u(θ P ) has a size of 128 × 128 × 128 voxels. However, to provide a greater anatomical variability to the dataset on which we train the surrogate, we crop all simulated outputs and corresponding brain anatomies to 64 × 64 × 64 volumes, centered at the initialization location x, y, z. The 64 × 64 × 64 is greater than half the brain size and tumors bigger than this are incompatible with life.
The architecture of the tumor model surrogate consists of two main parts: -Brain anatomy encoder which encodes 3D volumes of the brain tissues WM, GM, and CSF through a series of convolutional blocks. The blocks are composed of alternating convolution operations (with fixed parameters of kernel size 3, stride 1, and the number of maps 128) and a non-linearity in the form of a linear rectifier. Each block is equipped with a skip connection linking the input and output of the block via an element-wise sum. Downsampling between the blocks is achieved by a convolutional operation with a stride of two.
-Brain tumor decoder that takes as an input a 1D vector of the parameters {D, ρ, T } alongside with the encoded anatomy. Note that we do not condition the decoder on the initialization location x, y, z, since as mentioned above we crop all training volumes exactly at this location. Thus the network is taught to reproduce the tumor in the center for any volume. Before being passed to the decoder, the 1D vector of model parameters is mapped via a fully connected layer to a tensor of size 8 × 8 × 8 × 3 and is concatenated with the tensor of the encoded brain anatomy. The resulting tensor is gradually upsampled through a series of convolutional blocks analogous to the encoder and nearest-neighbor upsampling (we refer the reader to [27] for a detailed discussion on why such type of upsampling is preferred over the deconvolution operation). At the end of the series, a 4D tensor 64 × 64 × 64 × 128 is convolved to the output prediction -a 3D tumor simulation volume (height 64, width 64, depth 64). The decoder design is adopted from [23], [28].
As a loss function, we compute the error between the predicted (u pred ) and simulated (u sim ) cell concentration functions under L 1 norm separately in the CSF area and the area of non-zero ground truth simulated tumor (GT): As shown in the ablation analysis (Fig. 11), such loss definition significantly affects predictions performance in new patient geometry.
3) Bayesian model personalization: To demonstrate the applicability of the neural surrogate, we perform Bayesian tumor growth model personalization substituting the numerical solver with the learnable one. As calibration data we use two types of imaging modalities: (a) T1 contrast-enhanced and FLAIR MRI modalities that allow estimating the morphological characteristic of the visible tumor; (b) FET-PET scans that provide information about tumor metabolic activity.
Analogous to [14], [15], we relate the output of the tumor growth solver u(θ P ) to imaging information via a probabilistic model, where image observations D = {y T 1c , y F LAIR , y P ET } are assumed to be independent, and θ = {θ P , θ I } constitute parameters of the pathology model θ P and the probabilistic imaging model θ I . The latter is defined differently according to the type of modality: -MRI modalities provide information in the form of binary tumor segmentations (y T 1c , y F LAIR ). Thus we assign for each voxel a discrete label y M i ∈ {0, 1}, M ∈ {T 1c, F LAIR}. We model the probability of observing y T 1c , y F LAIR for a given concentration u(θ P ) with Bernoulli distribution, where α i,M is the probability of observing tumor-induced changes defined as a double logistic sigmoid, With this formulation, we postulate that the tumor is not visible on a scan below the threshold level u M c . The parameter σ α is introduced to take uncertainty in the threshold concentration u M c into account.
-FET-PET modality (y P ET ) provides continuous information in each voxel and can be assumed to be proportional to   the tumor density [3]. In this case we model the likelihood imaging function by a Gaussian distribution with unknown constant of proportionality b: Here, σ models uncertainty in the PET signal, which is considered to be normalized y P ET i ∈ [0, 1]. In total there are eleven parameters (six pathology model parameters θ P = {D w , ρ, x, y, z, T } and five imaging model , σ}) which we infer from the triplet of medical scans D = {y T 1c , y F LAIR , y P ET } using a Markov Chain Monte Carlo (MCMC) sampling algorithm [29]. 4) Implementation: The numerical tumor solver used for obtaining the simulation dataset is a highly parallelized glioma solver 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 β 1 = 0.5 and β 2 = 0.999 for 30 epochs, which was observed to be sufficient for convergence. The learning rate is cosine annealed from 10 −4 to 10 −6 over the training and the batch size is 16. We performed the experiments on an NVIDIA Quadro P8000 using the Tensorflow framework.
For Bayesian MCMC inference we use an implementation of Transitional MCMC from [32] with 2048 samples per iteration.

1) Data:
The patient-specific tumor-free brain anatomies (WM, GM, CSF) in which we model the tumor were obtained from an MRI dataset of 15 glioblastoma patients by registration-based approach proposed in [3], [33]. The data have resolution 256x256x256 with isotropic voxels of 1mm side-length, but for simulation purposes, the calibration is performed in data downsampled to 128x128x128. The diffusion coefficient in white matter D w is considered to be 10 times greater than in grey matter D g . The following ranges are used for random uniform sampling of the model parameters: D w ∈ [0.001, 0.08] mm 2 /day, ρ ∈ [0.001, 0.03] 1/day, and T ∈ [50, 1000] day with a step size of 50 days. The tumor location coordinates x, y, and z are sampled within the brain volumes. Samples that have initial location {x, y, z} within the CSF do not trigger any tumor growth and were thus discarded. For samples with initial location closer than 32 pixels to any of the volume borders, we do extra padding when cropping to 64x64x64 size. In total, we have a set of 20k parameterssimulation pairs for training.
2) Training the surrogate: Conventional application of tumor modeling such as model personalization (via solving an inverse problem) implies estimating the model parameters by sampling over fixed, physiologically plausible ranges. Thus, we aim to employ the learnable surrogate which i) possesses an interpolation capacity for the parameters {D, ρ, T }, and ii) is capable to extrapolate for simulations in new geometries. To probe these properties, the test set of 2k samples was formed to have only the brain anatomies unseen by the network during training, while the parametric {D, ρ, T } triplets were sampled from the same ranges as for the training. a) Quantitative convergence: Fig. 5 shows the convergence of the loss L total for both training and test sets. Fig. 6 demonstrates the distribution of the mean absolute error within each class of the brain tissues evaluated on the whole test set. Even though we observe samples with the error in the order of 10 −1 , the majority of the distribution lies within the order of 10 −2 . In Fig. 7 we depict the histograms of the DICE score computed between the simulated and predicted tumor volumes which are thresholded at different levels of the tumor cell concentration (such thresholding is exactly the operation that we perform during the Bayesian inference to relate tumor concentration with MRI signals, Eq. 6). The distributions are centered close to the DICE 1.0. The smaller peak at DICE 0.0 is due to the fact that thresholding of some simulations results in volumes containing a few (up to 5) non-zero voxels, whereas thresholding of the corresponding network predictions outputs volumes of all-zero voxels (or vice versa). In our particular application of the tumor model personalization such volumes are significantly smaller that the binary MRI segmentation volumes to which we calibrate the model and thus do not affect the Bayesian inference outcome. b) Qualitative convergence: Fig. 4 illustrates visually the surrogate's predictions for simulations from the test set. The predicted distribution of the tumor concentration accurately emulates the ground truth simulations. The global constraints imposed by the CSF anatomy (in which the tumor is constrained to grow) are well captured by the network. Fig. 8 depicts a comparison between a network prediction and a simulated ground truth in 3D.
3) Performance speed-up: The computing time for a single simulation using the numerical solver is on average 120 seconds using an Intel Xeon with 8 CPU cores and 64GB of RAM. The time required by using the neural solver for processing a batch of size 8 during inference is equal to 2.4s, i.e. our surrogate is 400 times faster than the numerical solver per single simulation. Such a speed-up reduces the computing time for applications such as model personalization from weeks to minutes. We want to point out that a theoretical comparison between CPU-and GPU-based computations should be taken with care for a few reasons. First, the neural surrogate does not need to do any sequential computations and can batch process multiple frames, whereas the numerical solver is successively integrating the PDE equation in time and for a single sample. Second, the CPU-based solver acts on a grid of dimension 128x128x128 whereas the network predicts the tumor in a cropped region of interest of 64x64x64 voxels. Finally, the convolutional layers are very efficiently realized in most GPU-tailored deep learning frameworks. In this regard, a fair comparison could imply to contrast the GPU compute with a CPU 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 and GPU platforms that are widely available. Also, we emphasize that when both solvers are tested on the application such as Bayesian personalization, we treat the empirical imaging data as well as the simulated domain identically keeping the size and resolution same (even though the network predicts 64x64x64 tumor volumes centered around the initial location {x, y, z}, they are then embedded in the original 128x128x128 simulation domain for the personalization task). 4) Bayesian model personalization in patient data: We performed the Bayesian brain tumor model calibration on preoperative scans of 5 test patients using the proposed neural surrogate and numerical solver (the max-a-posterior (MAP) estimates of the tumor density are provided in Tab.1 in the appendix). Tumor concentrations modeled with the MAP estimates are shown in Fig. 9. We observe an agreement of the glioma profiles obtained by the two methods and the variability of the estimations is within the variability of the Bayesian calibration. We also note that the learnable surrogate is trained on a dataset with continuous uniform distribution of 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 sampled from a continuous interval. This implicitly suggests that the network's interpolation capacity is sufficient for parametric inference.
Lastly, even though we tested the surrogate on the simplistic Fisher-Kolmogorov model, a translation of the method to more complicated tumor growth descriptions should be straightforward (through conditioning with additional model parameters) at no extra inference time cost.

IV. DISCUSSION
Despite the mentioned advantages, we would like to address the limitations of the proposed method. As depicted in Fig.  7, even though the Dice distributions are peaked close to a perfect score, there is a room for improvement to reduce the deviation. As seen in Fig. 4, the anatomy encoder captures the global constraints imposed by the CSF anatomy (in which the tumor is constrained to grow). However, conditioning of the network on the WM and GM maps improves moderately the performance as compared to training without such conditioning (see Fig. 7 and 10). This is somewhat expected as the WM and GM geometries affect numerical simulations in a highly non-linear fashion as compared to the CSF, geometry of which is preserved pixel-wise in the output tumor volume.
We experimented with several architectural and hyperparametric configurations and selected the one providing the lowest error on the validation set ( Fig. 10 and 11 show examples of the ablation analysis). Given that the space of network designs is close to infinite, future research could benefit from using 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 mentioned limitations and resulting approximation errors are acceptable for clinical translation. As well as we do not know whether the 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-surgery analysis, etc. [1]- [4]). However, we find the proposed method to be a solid baseline in the search for optimal tumor model surrogate, which in turn can significantly speed up our search for a biophysical model descriptive enough for clinical trials.

V. CONCLUSION
We present the first learnable surrogate with anatomy encoder for tumor growth modeling that is capable of mapping the model parameters to the corresponding simulations while respecting patient-specific anatomy. Our method achieves close to real-time simulation allowing fast inference in Bayesian settings. The technique can be readily adopted to more complicated tumor growth models and similar 4D inverse modeling tasks.