“Survival”: a simulation toolkit for radiobiological simulations for ion beam therapy

Andrea Attili (INFN)

  • XV Seminar on Software for Nuclear, Subnuclear and Applied Physics
  • May 29th, 2018 - Alghero, Italy

“Survival”

Software developed in C++ by the INFN (Istituto Nazionale di Fisica Nucleare) and University of Torino (UniTO, Physics Department) and provides different implementations of some radiobiological models to predict the cell survival after irradiation. The implemented models are (for the moment): LEMI, LEMII, LEMIII, MKM and MCt-MKM.

Reference:

  • Manganaro, L., Russo, G., Bourhaleb, F., Fausti, F., Giordanengo, S., Monaco, V., … Attili, A. (2018). “Survival”: a simulation toolkit introducing a modular approach for radiobiological evaluations in ion beam therapy. Physics in Medicine and Biology, 1–15. https://doi.org/10.1088/1361-6560/aab697

GitHub public repository:

Implemented models

  • LEM I: M. Scholz and G. Kraft, “Track structure and the calculation of biological effects of heavy charged particles”, Advances in Space Research 18, 5-14 (1996).
  • LEM II: T. Elsässer and M. Scholz, “Cluster effects within the local effect model”, Radiation Research 167, 319-329 (2007).
  • LEM III: T. Elsässer, M. Krämer and M. Scholz, “Accuracy of the local effect model for the prediction of biologic effects of carbon ion beams in vitro and in vivo”, International Journal of Radiation Oncology-Biology-Physics 71, 866-872 (2008).
  • LEM I,II,III rapid calculation GSI: M. Krämer and M. Scholz, “Rapid calculation of biological effects in ion radiotherapy”, Physics in medicine and biology 51, 1959-1970 (2006).
  • LEM I,II,III rapid calculation INFN: G. Russo (2011). PhD thesis
  • MKM + amorphous track structure: Y. Kase, T. Kanai, N. Matsufuji, Y. Furusawa, T. Elsasser, and M. Scholz, “Biophysical calculation of cell survival probabilities using amorphous track structure models for heavy-ion irradiation”, Physics in Medicine and Biology 53, 37-59 (2008).
  • MCt-MKM: Manganaro, L., Russo, G., Cirio, R., Dalmasso, F., Giordanengo, S., Monaco, V., … Attili, A. (2017). A Monte Carlo approach to the microdosimetric kinetic model to account for dose rate time structure effects in ion beam therapy with application in treatment planning simulations. Medical Physics, 44(4), 1577–1589.

Prerequisites

  • gcc or clang (C++11 compilers) + make utility
  • GSL (GNU Scientific Libraries).
  • OpenMP (Open Multi-Processing).
  • Git version control system and source code management.

Installation of the presequisites (for Ubuntu)

sudo apt-get install build-essential
sudo apt-get install libgsl-dev
sudo apt-get install libpomp-dev
sudo apt-get install git

Installation of the presequisites (for RedHat/CentOS)

sudo yum install gcc-c++ make
sudo yum install gsl-devel
sudo yum install libgomp-devel
sudo yum install git

furhter infos at https://github.com/batuff/Survival/blob/master/INSTALLATION.md

Download & compilation

Download

  • download via Git (to be installed if it is not present in the system)
git clone https://github.com/batuff/Survival
  • or directly wia web browser (“Clone or download” button at the GitHub page)

Compilation

  • with g++
cd Survival
make
./survival
  • or with clang (e.g. for Mac OSX)
cd Survival
make -f Makefile_clang
./survival

Documentation

Files useful for these exercises

Configuration: system variables configuration script

Edit the file setenv.sh to set the following system variables:

  • $install_folder: the folder where Survival has been downloaded and compiled (to be added to the variable $PATH)
  • $DATA: the folder where additional data can be found (e.g. SRIM and some tracks data). Usually found in $install_folder/data
ls data
Be_Srim.dat
C_Srim.dat
F_Srim.dat
H_Srim.dat
Ne_Srim.dat
O_Srim.dat
...
  • For the Mac OSX (High Sierra) there is a folder (ext_lib) to be added to $DYLD_LIBRARY_PATH containing precompiled libraries for GSL and OpenMP.

SRIM data

Contains information about the stopping power of different ions in water. It is used by Survival to perform conversions between dE/dx <–> Energy of these particles.

   Ion       dE/dx      dE/dx 
  Energy     Elec.      Nuclear
 (MeV/u)    (MeV/um)    (MeV/um)
---------- ---------- ----------
9.090e-04   5.345E-02  6.865E-02
1.000e-03   5.606E-02  6.648E-02
1.090e-03   5.856E-02  6.446E-02
1.181e-03   6.095E-02  6.258E-02
1.272e-03   6.325E-02  6.082E-02
1.363e-03   6.547E-02  5.917E-02
1.454e-03   6.762E-02  5.762E-02
1.545e-03   6.970E-02  5.617E-02
1.636e-03   7.172E-02  5.480E-02
1.818e-03   7.560E-02  5.228E-02
...         ...        ...

Usage

To execute the program the user has to call from the command line:

source setenv.sh
survival -SIMULATION_OPTION CHOSEN_VALUES ...

The user has the possibility to set a number of physical and biological parameters by using the syntax:

-PARAMETER_NAME PARAMETER_VALUE.

Running survival without arguments displays all the possible input arguments:

survival
Model Parameters:
   (MKM and tMKM)
     survival -MKM_alpha0 [double]
     survival -MKM_beta0  [double]
     survival -MKM_rNucleus [double]
     survival -MKM_rDomain  [double]
...

Arguments -- Outputs

  • -projectName It's a string representing the prefix to give at any file and directories that will be created in the simulation. The default value is “NewProject”.

  • -output The user has the possibility to choose between three kinds of output (and all possible combination between them):

    • LQ_pars a file will be created, named “projectname_LQparameters_model.csv”, containing the information about the parameters chosen for the simulation and the values of the simulated LQ α and β parameters (a new line for each energy evaluated).
    • meanValues a file will be created, named “projectname_survival_model.csv”, containing the information about the parameters chosen for the simulation and the values of doses delivered and survival observed (a new line for each energy or dose evaluated).
    • cellValues This kind of output is supported only by the MonteCarlo calculusType. (See next slides)
    • all to enable all output types.

Arguments -- Outputs

Arguments -- Models

  • -model It's a string representing the model to be used in the simulation. Five basic models are supported (see the reference paper):

    • LEMI
    • LEMII
    • LEMIII
    • MKM
    • tMKM_Manganaro2017

The default value for this option is “MKM”

Arguments -- Calculus

  • -calculusType A string identifying the type of calculus to be done. Some possibilities are available:

    • rapidLEM_Scholz2006 It's an implementation of the method described in: M. Krämer and M. Scholz, “Rapid calculation of biological effects in ion radiotherapy”, Physics in medicine and Biology 51, 1959-1970 (2006). It's compatible only with LEMI, LEMII and LEMIII models.
    • rapidLEM_Russo2011 A new rapid method for LEM, proved to be more accurate, described in: G. Russo, “Develpment of a radiobiological database for carbon ion Treatment Planning Systems - Modelling and simulating the irradiation process”, Ph.D. Thesis, Università degli studi di Torino (2011). It is compatible only with LEMI, **LEMII* and LEMIII models.
    • rapidMKM_Kase2008 A fast implementation of the MKM calculation as described in: Kase, Y., Kanai, T., Matsufuji, N., Furusawa, Y., Elsässer, T., & Scholz, M. (2008). Biophysical calculation of cell survival probabilities using amorphous track structure models for heavy-ion irradiation. Physics in Medicine and Biology, 53(1), 37–59 It's compatible only with the MKM.

Arguments -- Calculus (2)

  • -calculusType A string identifying the type of calculus to be done. Some possibilities are available:

    • rapidMKM_Kase2008_corrected_beta An extension of the rapidMKM_Kase2008 method in which a non-Poissonian correction factor is added also for the \( \beta \) parameter.
    • rapidMKM_Attili2013 A fast original implementation of the MKM model, combining the methods described in: Hawkins_2003 Hawkins, R. B. (2003). A microdosimetric-kinetic model for the effect of non-Poisson distribution of lethal lesions on the variation of RBE with LET. Radiation Research, 160(1), 61–69, and Kase, Y., Kanai, T., Matsufuji, N., Furusawa, Y., Elsässer, T., & Scholz, M. (2008). Biophysical calculation of cell survival probabilities using amorphous track structure models for heavy-ion irradiation. Physics in Medicine and Biology, 53(1), 37–59.
    • rapidMKM_Attili2013_corrected_beta An extension of the rapidMKM_Attili2013 method in which a non-Poissonian correction factor is added also for the \( \beta \) parameter. See Survival::Calculus::rapidMKM_Attili2017_corrected_beta() for details.
    • MonteCarlo Compatible with all models implemented, performs a Monte Carlo simulation of the irradiation process to get the LQ parameters. Compatible with all models.

Arguments -- radiobiological parameters (MKM)

  • -cellType A string identifying the name of the cell lline used in the calculation. The default value is “Cell1”.

Note: The cell line in reality is completely determined by the model parameters chosen. This is only a tag to indicate the cell but it isn't used in the simulation.

  • MKM Parameters

    • -MKM_alpha0 A double representing the linear quadratic α parameter characteristic for X-rays, expressed in Gy-1. The simulated α parameter will tend to this value for low LET.
    • -MKM_beta0 A double representing the linear quadratic β parameter characteristic for X-rays, expressed in Gy-2. The simulated β parameter will tend to this value for low LET.
    • -MKM_rNucleus The radius of the cell expressed in μm.
    • -MKM_rDomain The radius of domains wich constitute the MKM nucleus expressed in μm.
    • -MKM_timeConst The time constant associated to the repair kinetics of the cell.

Arguments -- radiobiological parameters (LEM)

  • LEMI, LEMII and LEMIII Parameters

    • -LEM_alpha0 A double representing the linear quadratic α parameter characteristic for X-rays, expressed in Gy-1.
    • -LEM_beta0 A double representing the linear quadratic β parameter characteristic for X-rays, expressed in Gy-2.
    • -LEM_rNucleus The radius of the cell expressed in μm.
    • -LEM_Dt The transition dose beyond which the standard linear quadratic parametrization became pure linear expressed in Gy.

Arguments -- source definition

  • -doses The sequence of doses to be delivered. The default is 1 to 6 Gy (step 1), in order to construct a survival curve.

  • Monoenergetic irradiation

    • -ion A string identifying the chemical symbol of the element, without mass number specifications. Ions from proton to neon are supported.
    • -energies A sequence of kinetic energies of the primary ions to be evaluated, expressed in MeV. This and the -lets option are mutually exclusive, but one of the two has to be specified.
    • -lets A sequence of LET of the primary ions to be evaluated, expressed in MeV. This and the -energies option are mutually exclusive, but one of the two has to be specified.
  • “Mixed fields” irradiation

    • -spectrum_file Useful for complex irradiation evaluation. This option provide the possibility to irradiate the cellular population with different ions and different energies in a mixed field, described in an external file. The user has to specify the name of the file containing the spectrum complete with its relative or absolute path.

Arguments -- source definition (Spectrum file format)

The spectrum file is a simple ascii file with 6 columns:

  1. Ion Type (1H, 4He, etc..)

  2. Charge

  3. A

  4. Kinetic energy (MeV)

  5. LET (keV/μm)

  6. Weight

Arguments -- source definition (temporal structure)

  • -nFraction Supported only by the tMKM model, it require an integer representing the number of fraction in which to divide each nominal dose to be delivered. The default is 1 (single fraction).

  • -timeSpacing Supported only by the tMKM model, it indicate the time spacing between consecutive fractions, expressed in hours. The default is 0 (no time spacing).

  • -fracDeliveryTime Supported only by the tMKM model, it indicate the fraction delivery time, expressed in hours. The default is 0 (istantaneous delivering).

Arguments -- Monte Carlo options

  • -precision Supported only by the MonteCarlo calculusType, it's a double identifying the precision to be reached in the calculation. Two possibilities are provided, as the user can indicate:

    • A positive integer: it will be taken as the number of cell to irradiate for each level of nominal dose imposed.
    • A double between 0 and 1: it indicates the relative error on the simulated survival to be reached before stopping the calculation.
  • -parallelismType A positive integer indicating the level of parallelism to be used in the simulation. The user has the possibility to specify the number of threads to dedicate at the calculation, in particular:

    • 1: Only 1 thread = parallelism disabled.
    • N>1: The number of threads to be used.
    • 0: Let the program to chose a number of thread corresponding to the number of core of the computer executing the program. This is set as default value.

Bash scripts for the exercises

Example of script (run-example1-a.sh):

#!/bin/bash

# Evaluate a simple survival curve with LEMI, for HSG irradiated with a monoenergetic beam

source Survival/setenv.sh

projectName="Example1-a"

output="meanValues"

model="LEMI"
calculusType="rapidLEM_Scholz2006"

precision=0.015
parallelismType=0

doses="0.5 1 2 3 4 5 6"

cellType="HSG"
LEM_alpha0=0.313
LEM_beta0=0.0615
LEM_rNucleus=5
LEM_Dt=30

ion="C"

trackMode="histogram"

lets="40"

survival -projectName $projectName \
        -output $output \
        -model $model \
        -calculusType $calculusType \
        -precision $precision \
        -parallelismType $parallelismType \
        -doses $doses \
        -cellType $cellType \
        -LEM_alpha0 $LEM_alpha0 \
        -LEM_beta0 $LEM_beta0 \
        -LEM_rNucleus $LEM_rNucleus \
        -LEM_Dt $LEM_Dt \
        -ion $ion \
        -trackMode $trackMode \
        -lets $lets \

Excercise 1a: Evaluation of a survival curve [HSG, LEM I]


  • Model: LEM I
  • LEM I parameters:
    • \( \alpha_0 = 0.313 \) Gy\( ^{-1} \),
    • \( \beta_0 = 0.0615 \) Gy\( ^{-2} \),
    • \( R_N = 5 \) \( \mu\text{m} \),
    • \( D_t = 30 \) Gy.
  • calculusType: rapidLEM_Scholz2006
  • Primary Ion: C
  • LET: \( \ 40 \) \( \ \text{keV}/\mu\text{m} \)
  • Doses: 1, 2, 3, 4, 5, 6 Gy
  • output: meanValues

Example file: run-example1-a.sh

Survival curve output format ("meanValues" output)

file: Example1-a_survival.csv

model calculusType cell alpha_0 beta_0 r_nucleus D_t radiation_type particle meanEnergy meanEnergy_sigma
LEMI rapidLEM_Scholz2006 HSG 0.313 0.0615 5 30 Monoenergetic C 687.241 0
LEMI rapidLEM_Scholz2006 HSG 0.313 0.0615 5 30 Monoenergetic C 687.241 0
LEMI rapidLEM_Scholz2006 HSG 0.313 0.0615 5 30 Monoenergetic C 687.241 0
meanLET meanLET_sigma LETd LETd_sigma nFractions timeSpacing fracDeliveryTime
0.04 0 0.04 0 1 0 0
0.04 0 0.04 0 1 0 0
0.04 0 0.04 0 1 0 0
nominal_dose mean_dose mean_survival mean_doseUncertainty mean_survivalUncertainty
0.001 0.001 0.9989036 0 0
1.000 1.000 0.3197383 0 0
2.000 2.000 0.0937445 0 0

Excercise 1b: Evaluation of (α, β) vs. LET [HSG, LEM I]


  • Model: LEM I
  • LEM I parameters:
    • \( \alpha_0 = 0.313 \) Gy\( ^{-1} \),
    • \( \beta_0 = 0.0615 \) Gy\( ^{-2} \),
    • \( R_N = 5 \) \( \mu\text{m} \),
    • \( D_t = 30 \) Gy.
  • calculusType: rapidLEM_Scholz2006
  • Primary Ion: C
  • LET: 10, 12, 15, … , 331, 406, 500 \( \ \text{keV}/\mu\text{m} \)
  • Doses: 1, 2, 3, 4, 5, 6 Gy
  • output: LQ_pars


Example file: run-example1-b.sh

Survival curve output format ("meanValues" output)

file: Example1-b_LQparameters_LEMI.csv

model calculusType cell alpha_0 beta_0 r_nucleus D_t radiation_type particle
LEMI rapidLEM_Scholz2006 HSG 0.313 0.0615 5 30 Monoenergetic C
LEMI rapidLEM_Scholz2006 HSG 0.313 0.0615 5 30 Monoenergetic C
LEMI rapidLEM_Scholz2006 HSG 0.313 0.0615 5 30 Monoenergetic C
meanEnergy meanEnergy_sigma meanLET meanLET_sigma LETd LETd_sigma nFractions timeSpacing fracDeliveryTime
5553.84 0 0.0100000 0 0.0100000 0 1 0 0
3678.81 0 0.0122863 0 0.0122863 0 1 0 0
2606.32 0 0.0150952 0 0.0150952 0 1 0 0
alpha beta alpha_sigma beta_sigma
0.578043 0.0563536 0 0
0.650056 0.0548854 0 0
0.721866 0.0533395 0 0

Excercise 1c: Evaluation of (α, β) vs. LET [HSG, MKM]


  • Model: MKM
  • MKM parameters:
    • \( \alpha_0 = 0.313 \) Gy\( ^{-1} \),
    • \( \beta_0 = 0.0615 \) Gy\( ^{-2} \),
    • \( R_N = 4.1 \) \( \mu\text{m} \),
    • \( R_d = 0.34 \) \( \mu\text{m} \).
  • calculusType: rapidMKM_Kase2008_corrected_beta
  • Primary Ion: C
  • LET: 10, 12, 15, … , 331, 406, 500 \( \ \text{keV}/\mu\text{m} \)
  • Doses: 1, 2, 3, 4, 5, 6 Gy
  • output: LQ_pars


Example file: run-example1-c.sh

Excercise 1d: Comparison between LEM I, MKM and exp. data

Excercise 2a-b: Comparison between Monte Carlo and fast approx. [HSG, MKM]


  • Model: MKM
  • MKM parameters:
    • \( \alpha_0 = 0.313 \) Gy\( ^{-1} \),
    • \( \beta_0 = 0.0615 \) Gy\( ^{-2} \),
    • \( R_N = 4.1 \) \( \mu\text{m} \),
    • \( R_d = 0.34 \) \( \mu\text{m} \).
  • calculusType: rapidMKM_Kase2008_corrected_beta
  • Primary Ion: C
  • LET: 40 \( \ \text{keV}/\mu\text{m} \)
  • Doses: 1, 2, 3, 4, 5, 6 Gy
  • output: LQ_pars meanValues
    Example file: run-example2-a.sh


  • Model: MKM
  • MKM parameters:
    • \( \alpha_0 = 0.313 \) Gy\( ^{-1} \),
    • \( \beta_0 = 0.0615 \) Gy\( ^{-2} \),
    • \( R_N = 4.6 \) \( \mu\text{m} \),
    • \( R_d = 0.34 \) \( \mu\text{m} \).
  • calculusType: MonteCarlo
  • Primary Ion: C
  • LET: 40 \( \ \text{keV}/\mu\text{m} \)
  • Doses: 1, 2, 3, 4, 5, 6 Gy
  • output: LQ_pars meanValues cellValues
  • precision: 0.15
  • parallelismType: 4
    Example file: run-example2-b.sh

Monte Carlo data output format ("cellValues" output)

Monte Carlo data output format ("cellValues" output)

file: Example2-b_survival.csv

model calculusType cell alpha_0 beta_0 r_nucleus r_domain radiation_type particle meanEnergy meanEnergy_sigma
MKM MonteCarlo HSG 0.313 0.0615 4 0.38 Monoenergetic C 979.415 0
MKM MonteCarlo HSG 0.313 0.0615 4 0.38 Monoenergetic C 979.415 0
MKM MonteCarlo HSG 0.313 0.0615 4 0.38 Monoenergetic C 979.415 0
meanLET meanLET_sigma LETd LETd_sigma nFractions timeSpacing fracDeliveryTime
0.0303 0 0.0303 0 1 0 0
0.0303 0 0.0303 0 1 0 0
0.0303 0 0.0303 0 1 0 0
nominal_dose mean_dose mean_survival mean_doseUncertainty mean_survivalUncertainty
1 0.9930091 0.4917271 0.1114103 0.0530305
2 1.9526080 0.2167977 0.1055546 0.0253573
3 3.2252050 0.0687924 0.1453497 0.0100640

Monte Carlo data output format ("cellValues" output)

file: Example2-b_survival_data/ 000_MonteCarlo_parameters.csv


PARAMETER_NAME VALUE
1 Model MKM
2 Parametrization LQ_noDt
3 Calculus MonteCarlo
4 Cell HSG
5 alpha_0 0.313
6 beta_0 0.061
7 r_nucleus 4.000
8 r_domain 0.380




9 Radiation_type Monoenergetic
10 Particle C
11 ParticleA 12
12 ParticleZ 6
13 Track_Type KieferChatterjee
14 Track_Mode histogram
15 Number_of_Fractions 1
16 Time_Spacing 0.000
17 Fraction_Delivery_Time 0.000

Monte Carlo data output format ("cellValues" output)

file: Example2-b_survival_data/ 000_MonteCarlo_survival_data/ dose_survival_30.30(keV-um)_1.00(Gy).dat

dose (Gy) survival
0.7146526 0.6139393
1.0488400 0.4238110
1.0031290 0.4842538
1.3100990 0.3492380
1.3288060 0.3731056
0.9499197 0.4629264
0.9246260 0.5170134
0.2349154 0.8918953
1.5776330 0.2439996
0.8374691 0.5570889

Excercise 2a-b: HSG irradiated with C @ keV/μm