Glioma Solver

Modeling Tumor Dynamics

In this tutorial we show how to:

1. Model tumor growth in patient-specific anatomy
2. Model tumor mass-effect and intracranial pressure
3. Use Bayesian inference to calibrate tumor model for patient-specific predictions
4. Tips & Tricks section offers additional useful information
5. Development section provides guidance on how to modify or extend the inference tools

1. Tumor Growth Modeling

This example simulates tumor growth in patient anatomy using a reaction-diffusion model. In this tutorial, we use synthetic brain anatomy, which can be downloaded here. The details of the model are explained in [Lipkova 2018] in References. For tissue segmentation tools please check the (References).

Input anatomy

The simulation is performed in the brain anatomy consisting of white matter (wm), gray matter (gm) and cerebrospinal fluid (csf). The GliomaSolver accepts the input anatomies in the form of a 3D matrix in the .dat-format. The folder GliomaSolver/tools/DataProcessing/source/ contains matlab scripts for converting medical nifty data from/to .dat. For each patient provide one folder with segmentation of wm, gm, and csf, each in a separate file with following name extension:

*WM.dat for wm
*GM.dat for gm
*CSF.dat for csf

The folder with the test data for this tutoruial contains multiple modalities and tissue segmentations. You can use the matlab script nii2dat to convert the tissue segmentation from nifty to .dat or use the already converted data in Atlas/anatomy_dat.

Executables

A) If you are using precompiled executables, you can decide between two options: brain128 or brain256, where the numbers 128 and 256 indicate the number of the simulation grid points per dimension.

B) If you installed the solver, you can set the grid resolution in the GliomaSolver/makefile/Makefile (see Documentation for details) and compile the solver:

cd GliomaSolver/makefile
make clean && make -j 4

to get the executable brain.

Simulation

The folder GliomaSolver/Simulations/TumorGrowth contains script run.sh, which allows to set-up the simulation input parameters and to execute the simulation. The input parameters are described below

The patient anatomy is mapped to the simulation grid represented by a unit cube. All spatial units are automatically rescaled. Due to the unit cube, the initial location of the tumor should be number between [0,1]. The tumor difusivity in the gray matter is assumed as Dg=0.1 Dw.

After setting up the input parameters, run the script as:

./run.sh

Output

If the input parameter vtk is set to 1, the solver will save the output in the form Data_000*.vtu, where * denotes the number of the output data in a chronological order, i.e. Data_0000.vtu is the initial condition and Data_0001.vtu is the solution at the time specified in dumpfreq. Figure 1 shows results of the tumor growth simulation with the adaptvie grid visualised with Paraview.

2. Tumor mass-effect and intracranial pressure model

This example simulates a tumor growth in the patient anatomy, accounting for the tumor-induced brain deformations and distribution of the intracranial pressure (ICP). In this example, we use synthetic brain anatomy which can be downloaded here. The manuscript with the details of the model is currently in a preparation, however, all the details are available upon request.

Input data

The input data should be provided in the .dat-format. See GliomaSolver/tools/DataProcessing/source/ for tools for converting nifty data to/from dat-format. For each patient provide one folder with segmentation of wm, gm, csf, and phase-field function (pff) of the brain mask, each in a separate file with following name extension:

*WM.dat for wm
*GM.dat for gm
*CSF.dat for csf
*PFF.dat for pff  

Compilation

Use the GliomaSolver/makefile/Makefile to set the grid resolution (see Documentation for details) and compile:

cd GliomaSolver/makefile
make clean && make helmholtz=hypre -j 4

to get the executable brain.

Simulation

The folder GliomaSolver/Simulations/BrainDeformations contains script run.sh, which allows to set-up the simulation input parameters and to execute the simulation. The solver for the brain deformation model does not support adaptive grid refinement, since in this case, the region of interest includes not only the tumor but also brain tissue and corresponding pressure distribution. However, the high performance of the software is ensured by additional MPI parallelisation. The simulation input parameters are described below

The simulation domain is treated as a unit cube. All spatial units are automatically rescaled to the simulation domain. Due to unit cube, the initial location of the tumor should be number between [0,1]. The tumor difusivity in the gray matter is assumed as Dg=0.1Dw.

Run the script as:

./run.sh

Output

If the input parameter vtk is set to 1, the solver will save output in the form Data_000*.vtu, where * denotes the number of the output data, i.e. Data_0000.vtu is the initial condition and Data_0001.vtu is the solution at the time specified in dumpfreq. For the content of the individual output channels see Documentation. Figure 2 shows the results of the tumor simulations.

3. Bayesian calibration for patient-specific predictions

This example shows how to calibrate the tumor growth model with respect to patient medical scans to obtain patient-specific predictions about the tumor cell density, which could be used for the radiotherapy planning (as presented in [Lipkova 2018] in References). The used medical data include structural MRI (T1Gd, FLAIR) and functional FET-PET scans. The inference procedure is automated and its usage is explained in the following subsections.

Input Data

For each patient provide one folder with the following input data in (nii or nii.gz format) and the following naming convetion:

• FLAIR.nii.gz = FLAIR MRI scan
• T1c.nii.gz = T1Gd MRI scan
• FET.nii.gz = PET-FET scan
• Tum_T1c.nii.gz = Binary segmentation of the T1Gd-enhanced tumor, including the necrotic core
• Tum_FLAIR.nii.gz = Binary segmentation of the FLAIR-enhanced tumor
• Tum_FET.nii.gz = FET-PET signal constrained to the tumor region only (obtained by subtracting the baseline signal)
• WM.nii.gz = Segmentation of the wm
• GM.nii.gz = Segmentation of the gm
• CSF.nii.gz = Segmentation of the csf

Figure 3 shows an overview of the input modalities and tumor observations. The test data for this tutorial can be obtained here and the set-up scripts are provided in GliomaSolver/Simulations/PatientInference/.

Executables

To speed-up the inference procedure, the model calibration is performed at the lower resolution (128^3 grid poitns), while the patient-specific prediction with the calibrated parameters is performed at high resolution (256^3). For this reason, one need to provide two brain executables, one for each resolution.

A) If you are using precompiled libraries, please copy brain128, brain256, and engine_tmcmc to the folder GliomaSolver/Simulations/PatientInference/. And rename brain128 to brain.

B) If you installed the solver, please compile the solver to get the brain executables for 128^3 and 256^3 resolution. To do so, first set the grid resolution in your GliomaSolver/makefile/Makefile as follows:

dim=3
bpd=16
bs=8
resjump=1

Then compile and copy the brain executable to the GliomaSolver/Simulations/PatientInference/ folder

 make clean && make -j 4
 cp brain ../Simulations/PatientInference

Afterward, create the high resolution executable by setting the bs=16, recompile and rename the executable brain to brain256 and copy it to the GliomaSolver/Simulations/PatientInference/.

Simulation

The folder GliomaSolver/Simulations/PatientInference/ should now contain following files:

• Input.txt = file for setting up the input parameters
• setupInference.sh = script for data preprocessing and setting up the inference enviroment
• brain = GliomaSolver executable at 128^3 resolution
• brain256 = GliomaSolver executable at 256^3 resolution
• engine_tmcmc = if using precompiled libraries

Open the folder Input.txt and set the following information:

DataPath=/home/jane/GliomaSolver/Simulations/Patient/InputData/
SolverPath=/home/jana/WORK/GliomaSolver/
Nsamples=2048

where DataPath= should contain the full path to the folder with the input nifty data. SolverPath= should contain the full path to the folder where you installed your GliomaSolver. Nsamples is the desired number of samples used for the Bayesian calibration. For the ideal performance, the number should scale with the number of processors. It is recommended to use at least 2000 samples.

After setting up the input file Input.txt call the setupInference script:

./setupInference.sh

This will:

1. Convert the data from nifty to dat and store them in the folder with _dat extension
2. Preprocess the data and map them to the simulation grid (in folder Preprocessing)
3. Create folder Inference that contains all necessary tools for the calibration for the given patient
4. Tell you how to run the inference

If you are using precompiled libraries, please copy the executable engine_tmcmc to the folder Inference. Now we are ready to run the inference simulation. To do so, please go to the folder Inference and use the script run_inference.sh (for UNIX computer) or run_inference_LRZ.sh (for SLURM cluster). The run_iference.sh takes as input parameter the number of cores on your computer. So if you have e.g. a machine with 32 cores, submit the simulation as follows:

./run_inference.sh 32

The simulation creates temporal files called tmpdir.*.* used for the computation (for details see Tips & Tricks sections). The script will compute the posterior probability distribution (pdf) of the model parameters w.r.t the patient input data, extract the results and store them in the folder Results . Please see the Tips & Tricks section below for information about a safe way for the simulation execution, restarts, and troubleshooting.

RESULTS

The results are stored in the folder Results, which is created at the same level as the Inference folder and contains the following results:

• MAP.nii = The most probable tumor cell density mapped to the patient input anatomy
• ReportFile.pdf = A pdf file reporting about
  • Simulation set-up
  • Table with the inferred parameters
  • Figure with the posterior pdf
  • Overview of the input data and the most probable tumor cell density.

4. Tips & Tricks

This section contains additional notes about data preprocessing, common mistakes and how to deal with them, as well as more details about the inference methods.

Segmentation corrections

The Inference tool predicts tumor cell density constrained by the patient-specific anatomy and the tumor observations available from the medical data. Errors in the tissue and tumor segmentations thus affect the model predictions. The most common and significant errors are:

• CSF segmentation is used as a boundary condition for the tumor growth, therefore:
  • If the hemisphere separations and/or ventricles are not segmented correctly (e.g. there are missing points), the model can infiltrate through such artificial pathways and lead to incorrect predictions. If needed, correct the CSF segmentation, mainly in the regions close to the tumor.
  • Please remove CSF segmentation from the region of the tumor segmentations, otherwise, they will act as boundary conditions.
• Tumor segmentations: should be connected. Please fill artificial gaps in the segmentations and remove isolated non-tumor voxels if present.
• Tum_FET should contain only the uptake values from the tumor. The best way is to subtract the baseline signal and then use a binding box around the tumor to exclude signal from non-tumor parts.

TMUX - safe simulations execution

It is highly recommended to use job handling system like tmux, which allows you to run the simulation even if you disconnect from the computer. For tmux installation on the Ubuntu execute the first line, for MAC the second one:

sudo apt-get install tmux
sudo port install tmux

To submit a simulation via tmux, go to the folder with your executable. Start the tmux environment and submit your simulation i.e.:

tmux
./run_inference.sh 32

To close the tmux press control+b and then d. To reopen the tmux, type in terminal tmux a. You can list all open tmux session, and stop them based on the id as follows:

tmux ls
$ 0: 1 windows (created Mon Feb 19 12:28:01 2018) [90x20]
$ 1: 1 windows (created Mon Feb 19 12:29:32 2018) [114x31]
tmux kill-session -t 0
tmux kill-session -t 1

For more info about tmux usage see here.

Bayesian Inference details

This section provides additional information for the Bayesian inference. These can be useful if you wish to understand or modify the inference tool or if you are facing some issues.

Posterior PDF

The Bayesian inference computes the samples from the posterior probability distribution (pdf) of the model parameters w.r.t the patient input data. The samples are obtained by the TMCMC algorithm, which generates a sequence of intermediate distributions (called generations), each stored in the file called curgen_db_*.txt, where the * stands for the number of the intermediate distributions. For details on the TMCMC please see [Ching 2007] in References.

Convergence

The intermediate distributions converge from prior pdf to posterior. The convergence is marked by the parameter p, (p=0 for prior and p=1 for posterior), which is computed automatically. To value of the parameter p is reported in the file Inference/runinfo.txt.

Troubleshooting & Restart

A) If your simulation crashed for some reason (e.g. power supply issues), restart it by calling the run_inference.sh script again. This will use the samples from the last completed curgen_db_*.txt file to continue the simulation

B) If your simulation gets stucked, i.e. the inference method didn’t converge yet, but it stops creating new temporal folders tmpdir.*.*. Then stop the simulation and resubmit it again. This can happen in some cases because of numerical instabilities and/or memory shortage. We are looking for possible improvements.

5. Development

The PatientInference set-up is based on the work presented in [Lipkova 2018]. This section describes individual parts of the Inference pipeline and explains how to proceed if you wish to modify the inference, e.g. to use other modalities or tumor model. You can still use the ./setupInference.sh and run_inference.sh scripts, but you need to do some modifications.

Preprocessing

The preprocessing part first converts all nifty data to dat-format. These are then mapped to the GliomaSolver grid using the Glioma_UQ_DataProcessing.cpp file, to ensure that the model predictions and input data are in the same reference space. This is perfomed by the script Prepreocessing/mapDataToMRAGrid.sh. If you introduce new modality, you will have to modify the Glioma_UQ_DataProcessing.cpp file so that all modalities are read-in and mapped to the GliomaSolver grid.

Likelihood

Folder Inference/Likelihood/ contains files used for the likelihood computation. If you introduce new data or model, you will have to incorporate them to the likelihood computation. Afterwards compile your new likelihood as:

cd Inference/Likelihood/makefile
make clean && make

which creates the executable called likelihood.

Prior Range

Next, you need to set the number of all unknown parameters that should be inferred and their prior range. This is done in file Inference/tmcmc_glioma.par, which looks like this:

# TMCMC and SN-TMCMC configuration file
# problem dimension
Nth             3
# Max number of generations
MaxStages       1000
# Population size
PopSize         2048
# Chain Lengths
MinChainLength  0
MaxChainLength  1

# Boundaries
Bdef             0          4
B0               -8.5172     -4.3428
B1               -5.9145     -1.6607
B2                5.1417     9.0766

Here you will need to modify the following items:

• problem dimension - the number of the unknown parameters
• Population size - the number of samples that should be used to populate the posterior pdf
• Boundaries - the prior range for all the unknown parameters. Here B0 contains the upper and lower bound for the first parameter and so on. Please provide the prior range for all your parameters. In some cases, it is convenient to work with a logarithm of the prior range, as done in this example.

Inference Files

Folder Inference/TMPTMP should contain all data and scripts needed for running the tumor model and computing the likelihood function, i.e. it should contain:

• brain = executable for the tumor model
• likelihood = executable to compute the likelihood function
• *.dat = all the available tumor observations used for the model calibration (e.g. tumor segmentations..).
• ROI.dat = (optional) 3D matrix marking the regions that should be used for the likelihood computation. It can be either brain mask or bounding box around the tumor.
• runAll.sh = execution script that runs the tumor model with the input parameters provided by the TMCMC and compute the corresponding likelihood function

The *.dat data, ROI.dat and runAll.sh are created by the Prepreocessing/mapDataToMRAGrid.sh script. Please copy all the necessary data into the Inference/TMPTMP and modify the runAll.sh script so it is compatible with your tumor model.

Inference function

Fie Inference/fitfun_glioma.c is used to collect the parameter values proposed by a sample of the TMCMC algorithm, perform the tumor simulation with the proposed parameters, compute the likelihood function and report the values back to the TMCMC method. You have to modify this file so it is compatible with your model and data, where the part that should be modified is marked as follows:

/*--------------USER INPUT STARTS HERE------------------- */
    ...
    param[]  = exp(x[])
    ...
    FILE *finp = fopen("InputParameters.txt", "w");
    for (i = 0; i < 3; i++) fprintf(finp, "%.16lf\n", param[i]);
    fclose(finp);
    ...
/*--------------USER INPUT STOPS HERE------------------- */    

The parameter values proposed by TMCMC are stored in the array x[]. If your prior range is in the log space, you will have to convert the array x to the real space (as indicated in the snippet above). Next, use the proposed parameter values to write the input scripts for the tumor model and the likelihood. The current inference writes the model parameters into the file InputParameters.txt (see the snippet above), tumor initial location to TumorIC.txt and the likelihood parameters to the LikelihoodInput.txt. Afterwards, the script runAll.sh is called.

Compilation

Finally, compile and run your simulation:

cd Inference
make clear; make clean; make
./run_inference 32