ACCOMPANYING WEBSITE
Simon D. Angus & Monika Joanna Piotrowska (2014)
A matter of timing: identifying significant multi-dose radiotherapy improvements by numerical simulation and genetic algorithm search

Contents

Overview

This website provides accompanying material for the study listed above. Specifically, we provide here access to the full code to run the model and produce all data that stands behind the study, with examples for the two main modes of code use. In addition, all data either used to calibrate and validate the model or produced by model simulation runs can be downloaded. (Where multiple random seeds were run (as explained in the paper), we have provided in the data download only a single seed to reduce the size of the download.)

The main purpose of providing these materials is to allow exploration of the codes used to accomplish the study, and to explore the data sets behind the study. We have not provided these materials with full replication in mind (though this is entirely possible with the materials here).

We are happy to provide further details and limited support to users who wish to explore the model or data further in MATLAB.

The interested user should consult the documentation provided in every code file provided in the code base.

Download Code & Data

Please fill out the following details and click the 'download' button and you will be taken to the download page where you can get links for the code or data as you would like. Note: * indicates required field.

Outline of the code

Below we present a listing of the non primitive functions included in the model code folder. In addition, the model utilises the parallel and statistics toolboxes.

AngusPiotrowska2014_code-online/
  setup.m              % run first to setup paths
  tumor_evo.m          % the main calling function used from the prompt
  tumor_run.m          % the core routine at the heart of the model
  cleanup.m            % run at conclusion of session to remove paths
   lib-tumor/           % codes focussed on the underlying model with irradiation
     InitializeRunStats.m
     ReplenishSubstrate2.m
     Diffusion2.m
     IrradiateTumor_evo.m
     DelayTime.m
     SetLocalStates.m
     KillSites.m
     MetaboliseNutrients.m
     SetGlobalStates.m
     RepairCells.m
     ProbDeathInRepair.m
     CellDivision2.m
     BirthCells.m
     GetStats2.m
     GetDiams2a.m
     FinishCheck.m
   lib-opt/             % codes focussed on the GA Search procedure
     GetBenchmarks.m
     GetInitialPop.m
       LoadLib.m
       EvolvePopulation.m
         tournament.m
         topdecile.m
         topfive.m
         localone.m
         mutate.m
           randpermNC.m
         mutate_local.m
         crossover.m
         sol2prot.m
         checklib.m
     InitializeEvoRun.m
     ImportTumorState.m
     CompileStats.m
     ReportProg.m
     AddToLib.m
       UpdateFitness.m
         cells2fin_cells.m
     SaveLib.m
       vis_lib.m
     max_sol_len.m
   lib-ca/              % Cellular automata methods
     NumericalDiffusion2b.m
     moore.m
   lib-other/           % generic, enabling functions
     choose.m
     roulette.m
     ols.m
     normgen.m
     normdist.m
     evalfile.m

Outline of the Data

Below we present a listing of the data contained in the download, by folder, referring to each study conducted, in order, in the paper.

Calibration & Validation Studies to Obtain Optimal \tau
-------------------------------------------------------
  Opt.Calib.Sugie2006
  Opt.Calib.Otsuka2011
  Opt.Calib.CombineSO
  Opt.Valid.KalHahn
Benchmark I and II in-depth Studies
-----------------------------------
  Opt.Benchmark.BMI
  Opt.Benchmark.BMII
GA Search (BMI and BMII contexts)
---------------------------------
  Opt.Search.BMII
  Opt.Search.BMI
In-depth Study of Top 3 Candidates from GA Search
-------------------------------------------------
  Opt.Best.BMI
  Opt.Best.BMII
Periodic Candidate Study
------------------------
  Opt.Periodic.BMI
  Opt.Periodic.BMII
Full Data Studies of Top Periodic & GA Candidates
-------------------------------------------------
  Opt.FullData
'Tremble' (robustness) Experiments Around Top Candidates
--------------------------------------------------------
  Opt.Tremble.BMI
  Opt.Tremble.BMII

EXAMPLE: Running a 'static' (full-data, 1 generation) trial

Since we have asked only for a single generation (.n_gn=1), the model will switch into 'single-run' mode and will not conduct GA Search. Further, since we have asked for full data to be taken (.keep_full_stats=1) full, rich data will be captured by the run.

First, navigate to the main code folder and run setup,

cd ../../AngusPiotrowska2014_code-online/
setup

to set up the paths required for the model. Next, move to the example/ folder,

cd example/

Now, we are in a position to run the model with a given runfile that typically is located in the current working directory (in this case example/). See the tumor_evo documentation for the main file tumor_evo to familiarise yourself with it and its primary arguments.

The key input is the runfile_example_static.txt runfile which tells tumor_evo how to behave for the given experiment.

type runfile_example_static.txt
% -------------------------------- %
% RUNFILE -- STATIC STUDY EXAMPLE
%
%    S. Angus AND M. Piotrowska
%
% Note:
%   by setting .n_gn = 1, we do not engage the GA Search.
%   by setting .keep_full_stats = 1, full data will be taken.
%
% S. Angus AND M. Piotrowska
% -------------------------------- %
P_EVO.resol_hr = 1;
P_EVO.days_pr   = 4/24;	% days to run protocol for
P_EVO.days_post = 4/24; % days to keep growing tumor post-protocol
P_EVO.n_gn = 1;			% num. generations
P_EVO.n_tr = 10;			% num. trials (unique tumors) to test a solution on per generation
P_EVO.n_pop = 1;		% number of solutions in each population
P_EVO.tot_dose = 20;	% max total dose over treatment time
% -------------------------------- %
P_EVO.keep_full_stats = 1;
P_EVO.n_workers = 10;
P_EVO.lib_fname = 'GA_LIBna.mat';
P_EVO.BM_fname = 'bm1_parallel.mat';
% -------------------------------- %
P_EVO.first_dt = 6/3600;		% ensure first wait-time is this
P_EVO.prot_lib.R  = [10 20];		% library of possible fraction dose levels (Gy)
P_EVO.prot_lib.dt = [6/3600 0.25 0.5 1 2 4 6];  % library of possible wait-times between fractions (hr)
P_EVO.evo.pr_mut = 0.05;	% prob. of mutation to solution component
P_EVO.tum_lib_s = 'vnw';	% encoding of tumor lib to use
% --- initial solution strings --- %
P_EVO.len_v = max_sol_len(P_EVO);
P_EVO.ini.pop = [ ...
	[[1 2] zeros(1,P_EVO.len_v-2)]
    ];
% -------------------------------- %
% particulars of tumor dataset
P_EVO.L = 51;	% grid side length
P_EVO.tstep = 6;		% (s) frequency to update model
% -------------------------------- %

We pass a value for $\tau$ (the multi-fraction calibration parameter) and a random seed which is typically the replicate number (so that we can easily re-run the experiment with exactly the same random number stream),

tau = 0.8;      % the optimal fitting \tau value
seed = 1;       % random seed for this replicate

Now, we can run the main model. This time, the runfile will run a single generation of the model only, and so, the GA Search process will not be invoked. And, we've asked (in the runfile) for full data to be taken from the model run for detailed analysis.

tumor_evo('runfile_example_static.txt',tau,seed)

After the run has finished, it will have produced the file, runfile_example_static.txt_tau0.8_s1-out.mat in the current working directory. We can inspect this as normal,

load runfile_example_static.txt_tau0.8_s1-out.mat

The results data structure has the form,

disp(res)
          P_EVO: [1x1 struct]
         PARAMS: [1x1 struct]
            sol: [1x1 struct]
    full_params: [1x1 struct]

With the summary data held in res.sol.stats,

disp(res.sol.stats)
    finish_time: [10x1 double]
         Ncells: [10x8 int32]
      CellState: [1x1 struct]
           NecR: [10x8 single]
           TumR: [10x8 single]
    PotProlifPc: [10x8 single]
       ProlifPc: [10x8 single]

EXAMPLE: Running a genetic algorithm (GA) 'search' (multi-generation, limited stats)

NB: Prior to running the GA Search mode of the model you must ensure there is no '.lock' file in the current working directory associated with the runfile '<runfile_name>.lock' (i.e. which was generated from a previous run), or the GA will get stuck, assuming that another search node presently is using the library .mat file.

To run the GA Search conduct the same setup as before if needed, setup, cd example, tau = 0.8; seed = 1;

Again, the runfile carries the key information required to set up the run,

type runfile_example_search.txt
% -------------------------------- %
% RUNFILE -- GA_Example.mat
% S. Angus AND M. Piotrowska
% -------------------------------- %
P_EVO.resol_hr = 1;
P_EVO.days_pr   = 1/24;	% days to run protocol for
P_EVO.days_post = 1/24;    % days to keep growing tumor post-protocol
P_EVO.n_gn = 3;		% num. generations
P_EVO.n_tr = 10;			% num. trials (unique tumors) to test a solution on per generation
P_EVO.tot_dose = 10;	% max total dose over treatment time
P_EVO.n_pop = 2;		% number of solutions in each population
% -------------------------------- %
P_EVO.keep_full_stats = 0;
P_EVO.n_workers = 10;
P_EVO.lib_fname = 'GA_LIB_Example.mat';
P_EVO.BM_fname = 'bm1_parallel.mat';
% -------------------------------- %
P_EVO.first_dt = 0.01;		% ensure first wait-time is this
P_EVO.prot_lib.R  = [1.25];		% library of possible fraction dose levels (Gy)
P_EVO.prot_lib.dt = [0.01 0.02 0.03];  % library of possible wait-times between fractions (hr)
P_EVO.evo.pr_mut = 0.10;	% prob. of mutation to solution component
P_EVO.evo.lcl1_steps = 4;	% # steps up and down vector to search locally
P_EVO.evo.rank_type = 'RAND';	% first candidates arise randomly
P_EVO.tum_lib_s = 'vnw';	% encoding of tumor lib to use
% --- initial solution strings --- %
P_EVO.len_v = max_sol_len(P_EVO);
P_EVO.ini.pop = [ ];
% -------------------------------- %
% particulars of tumor dataset
P_EVO.L = 51;	% grid side length
P_EVO.tstep = 6;		% (s) frequency to update model
% -------------------------------- %

And then run the GA Search with,

tumor_evo('runfile_example_search.txt',0.8,1)

During the simulation, a .lock file will appear and disappear as the GA performance library is augmented by the successive generations. Note that it is possible to run the model on multiple nodes, each with multiple CPUs (cores) so long as the nodes access the same file structure by cluster-based linking. Dropbox or similar may be appropriate, but file transfer rates may lag the rate of library update and manipulation. Some experimentation may be required if running in this mode.

As each generation finishes (on a processing node), the GA performance library will be updated. In this case, the file is called GA_LIB_Example.mat.

In addition, each library update, a performance output file will be created which provides a convenient summary of the current top- performing (by p-value) candidates in the library. A copy of the current performance log file follows,

type performance-log.txt

Update @ Oct.14,2014 16:25:41 (Melbourne time)

LIB contains 6 candidates. Showing top 6 (by p-val) ...
 -------------------------------------------
  id   avg_s   min_s   max_s   p-val    type
 -------------------------------------------
   1  238.00    9.46   22.23    0.00    RAND
   2  238.00    9.62   22.07    0.00    RAND
   3  238.00    9.62   22.23    0.00    TMNT
   4  238.00    9.46   22.23    0.00    TMNT
   5  238.00    9.47   22.23    0.00    TPDC
   6  238.00    9.92   22.15    0.00    TPDC
 -------------------------------------------

Last 6 computes took (avg,min,max): 0, 0, 0 (min).

To explore the library, we can do as normal,

load GA_LIB_Example.mat
disp(LIB)
1x6 struct array with fields:

    ttl_days
    T
    R
    s_ft
    s_n
    s
    avg_s
    p_val
    inh
    x
    cells_v_t
    comp_time
    Params

Each location in the Library contains an individual candidate that has been evaluated by the model. We can inspect one individually by using the index, e.g.,

disp(LIB(1))
     ttl_days: 0.0833
            T: [0.0100 0.0200 0.0400 0.0600 0.0700 0.0800 0.1100 0.1400]
            R: [1.2500 1.2500 1.2500 1.2500 1.2500 1.2500 1.2500 1.2500]
         s_ft: [10x1 double]
          s_n: [10x1 double]
            s: [10x1 double]
        avg_s: 238
        p_val: 0
          inh: [0 0]
            x: [3 1 1 1 2 1 2 1 1 1 1 1 3 1 3 1]
    cells_v_t: [10x2 int32]
    comp_time: 2.9310e-04
       Params: [1x1 struct]