Simulation of Tumor Cell Population Growth Dynamics

Research Group at Shepherd - Dr. Qing Wang (Mathematics, PI), Dr. Zhijun Wang (Computer Science).

Project Mentor - Dr. David Klinke, Department of Chemical Engineering, Department of Microbiology, Immunology & Cell Biology, West Virginia University.

Experimental Data by Dr. Jonathan Bramson, Department of Pathology and Molecular Medicine, McMaster University, Hamilton, Ontario, Canada

Goal of Simulation

Find the close-to-optimal, hopefully the optimal, parameter set in a high-dimensional parameter space using a computationally tractable method. Please see handouts for the model. Test set comes from experimental data.

Parameter Set Optimization

Genetic Algorithm. Encode with the parameter set. Evolve using crossover and mutation. For each generation, the ODE set is solved for each individual and the fitness function value (variance here) is calculated. The top half is passed to the next generation. Note the best individual in each generation is free from mutation. Crossover point(s) and mutation rate can be dynamically adjustable parameters so the algorithm can move more efficiently in the large parameter space and hopefully not get trapped in not-so-good local minimums.

Simulated Annealing + AMCMC. If applicable to the model, good estimates of initial parameter values may be calculated using simulated annealing. Then this initial parameter set is used as the starting point for a Markov chain. Please refer to the methodology article from Dr. Klinke [1].

Initial Value ODE Solvers

Runge-Kutta (order 4). Adams-Moulton formulas (non-stiff problems). Backward differentiation formulas (stiff problems). Any ODE solver should be implemented in the most efficient way due to the amount of calls in either the genetic or the AMCMC method. The ODE solver is the computational bottleneck of the simulation.

Please note: All variable and parameter values are equally weighted during the initial phase. The relative weights can be changed easily.

Design and Implementation Plan

Parameter ranges,  initial values, and simulation parameters can be specified on the interface. Final converged or time-up parameter values and the numerical solution to ODEs are saved in a text file and used for analysis and drawing.

Please note that functional curves for variables can be drawn using the saved data file. There are quite a few variables and they have different units so they can not be drawn in the same graph on the simulation interface as before.

The simulators in the following table are hosted and tested from this page. Please note the simulation process will consume client machine's resources. Internet Explorer with .Net framework 4 and/or Visual Basic Power Pack are required to run the application.

  07/04/2014 - New Holiday Edition
Equation Set Solver Version Notes
Blue 2 with M Sample impulse parameter files - Impulses_LV. Impulses_M.
Please put above files on Windows Desktop. Parameters can be modified.
  • Added effect of chemotherapy. Please refer to handouts for equations and variable/parameter meanings.
  • Impulses for chemotherapy drug M are read from a file. The sample file can be modified.
  • Added new variable M and its equation with an initial value of 0. Parameter kd8 was added.
  • Equations for TE1, TE3, CMHCI+, and CMHCI- were modified to reflect the effects of M. Parameters d, g1, g2, g3, and g4 were added.
  Second Generation Model - Possible Future Work
Simulator Equation Set Solver Version Notes
Violet 1.   Violet 1  Linear Violet 2
  • New second generation model. Please refer to handouts for equations and variable/parameter meanings.
  • Linear versions are used to fine tune the parameter values only. Logarithm based scales should be used for initial searches in a big parameter space.
  • Fitness remains the same as the indigo versions below.
  New Fitness with Separated Components - Possible Future Work
Simulator Equation Set Solver Version Notes
Indigo 1.  Indigo 1 Linear Indigo 2
  • New first generation model. Please refer to handouts for equations and variable/parameter meanings.
  • Linear versions are used to fine tune the parameter values only. Logarithm based scales should be used for initial searches in a big parameter space.
  • Fitness is separated into two components - 
    1. The normalized error squared (Divided by max of experimental data squared).
    2. The normalized differential slope squared (Divided by max of experimental data slope squared).
  11/22/2013 - First Generation Model Update
Simulator Equation Set Solver Version Notes
Blue 1.   Blue 1  Linear Blue 2 (Sample impulse parameter file)
  • New first generation model. Please refer to handouts for equations and variable/parameter meanings.
  • Updated 11/22/2013. Added impulse control to LV. Sample impulse parameter file can be downloaded and modified to an impulse pattern in the format of day - amount. Total number of days is now a parameter from user input.
  • Updated 10/27/2013.

    • Variance v = a * v1 + b * v2. Parameters a and b are positive and determined by user. gp100 data is not included.

    • v1 =  sum of difference of slopes squared of the ODE solution and corresponding experimental data (Normalized by the max of the slope of experimental data squared for each variable).  

    •  v2 = sum of error squared of the ODE solution and corresponding experimental data (normalized by the max of the experimental data squared for each variable).

    • TE3 is excluded in v1 as TE3 only has one day of data (implemented in the previous version).

  • Updated 10/20/2013.
  • Linear versions are used to fine tune the parameter values only. Logarithm based scales should be used for initial searches in a big parameter space.
  • Fitness remains the same as versions below.
  6/6/2013 - New Fitness with Components, New Version of Equations (1.5 & 1.6), etc.
Simulator Equation Set Solver Version Notes
Green 1 
Green 1 Linear		 Green 2
  • Equations (1.1) - (1.9) not changed from the last version (Yellow 3 & 4).
  • From w/o Volt in eq. 1.5 versions (Originally from Orange 5 and 6).
Green 3
Green 3 Linear		 Green 4

Following changes were made to Green 1 & Green 2.

  • (1.5)  CMHCI+ / Volt is changed to CMHCI+ / (r1+CMHCI+).
  • (1.6)  Term - r2 CMHCI-2 added.
  • Updated 7/22/2013. kp2 – kd4 is changed to 0.5. Tn(0) is changed 0.0714. Cmhcm(0) is changed to 2E6.
  • Updated 7/18/2013. Linear parameter ranges were added to fine tune parameter values in small ranges.
  • Updated 7/16/2013. Due to the limitation of experimental data, slope variance component is not calculated for TE3. LV parameter kd2 is fixed to 0.37 in Green 3 & 4.
  • Updated 6/27/2013. Calculated TE3 will not convert using the formula (parameter β1 will not be used). Will fit LV so parameter kd2 will be "unfixed".
  • Updated 6/20/2013. #2 is not calculated if set to 0.
  • Removed term kp * TE1 from the right side of (1.2).
  • Option "Decouple Eqs 1.1-1.4" added. When this option is checked, please ignore all other equations and their parameters since they have no effect on 1.1-1.4.
  • For versions Green 1 - 4 above, please put your experimental data file "exp_data.txt" on the desktop, or else the built-in default experimental data file will be used.
  • LVpower / ( γ + LVpower) for all 3 terms in 1.1 and 1.2.
  • Fitness has the following two components with adjustable relative weights ("Slope Portion" indicates the portion of fitness component #2 below) -
    1. The normalized variance (Divide by max of experimental data squared. Apply (x-min)/(max-min). Average over number of data points. In [0, 1]).
    2. The normalized differential slope squared (Divide by max of experimental data slope squared. Apply (x-min)/(max-min). Average over number of data points. In [0, 1]).
  12/20/2012 - New Terms (last update - 2/17/2013 with new fitness/variance)
Simulator Equation Set Solver Version Notes
Yellow 1 Yellow 2 Equations (1.1) - (1.4).
Set fitness/variance thresholds for TE1 and TE2 respectively.
TE1 or TE2 has
1. The previous fitness definition using normalized ratios (always ≥ 1).
2. The normalized (divided by the maximum of experimental data squared and the number of data points) variance  (always ≥ 0).
Yellow 3 Yellow 4 Equations (1.1) - (1.9). Not changed.
  • Added new term kp * TE1 to the right side of (1.2).
  • For versions Yellow 1 - 4 above, please put your experimental data file "exp_data.txt" on the desktop, or else the built-in default experimental data file will be used.
  • LVpower / ( γ + LVpower) for all 3 terms in 1.1 and 1.2.
  • From w/o volt in eq. 1.5 versions (Modified from Orange 5 and 6).
  12/7/2012 - Fitting the Curves — Complete Model Edition
Simulator for 1.1-1.9 (from w/o volt in eq. 1.5 versions) Equation Set Solver Version Notes
Orange 1 Orange 2 LV2 / ( γ + LV2) for all 3 terms in 1.1 and 1.2.
Orange 3 Orange 4 LV3 / ( γ + LV3) for all 3 terms in 1.1 and 1.2.
Orange 5 Orange 6 LVpower / ( γ + LVpower) for all 3 terms in 1.1 and 1.2.
For versions Orange 1 - 6 above, please put your experimental data file "exp_data.txt" on the desktop, or else the default experimental data file will be used.
  11/23/2012 - Fitting the Curves Edition
Simulator for 1.1-1.4 (from w/o volt in eq. 1.5 versions) Equation Set Solver Version Notes
Red 1 Red 2 LV for all 3 terms in 1.1 and 1.2.
Red 3 Red 4 LV / ( γ + LV) for all 3 terms in 1.1 and 1.2.
Red 5 Red 6 LV2 / ( γ + LV2) for all 3 terms in 1.1 and 1.2.
Red 7 Red 8 LV1/2 for all 3 terms in 1.1 and 1.2.
Red 9 Red 0 LV3 / ( γ + LV3) for all 3 terms in 1.1 and 1.2.
For versions Red 1 - 0 above, please put your experimental data file "exp_data.txt" on the desktop, or else the default experimental data file will be used.
Simulators - 11/19/2012 Version Notes
         w/o volt in eq. 1.8 ← Violet 6 Equation solver only. LV / ( γ + LV) ---> LV ^ 1/2 for all equations. Others same as the last row. Updated 11/20/2012.
       w/o volt in eq. 1.8 ← Violet 5 Special version (1.1-1.7 only).
       w/o volt in eq. 1.8 ← Violet 4 Equation solver only. LV / ( γ + LV) ---> LV for all equations. Others same as the last row. Updated 11/20/2012.
       w/o volt in eq. 1.8 ← Violet 3 Special version (1.1-1.7 only).
       w/o volt in eq. 1.8 ← Violet 2 Equation solver only. LV ---> LV / ( γ + LV) in (1.1). Equation for TE1 is decomposed into 4 new equations. Terms with parameters α and a21 are removed from equations. TE1 is changed to TE1d in (1.7).
       w/o volt in eq. 1.8 ← Violet 1 Special version (1.1-1.7 only).
 
Simulators - 11/13/2012 Version Notes
  w/ volt in eq. 1.5 ← Blue 7       w/o volt in eq. 1.5 ← Blue 8 Equation solvers only. LV ---> LV ^ 1/2 for all 3 terms in (1.1) and (1.2). Updated 11/20/2012.
w/ volt in eq. 1.5 ← Blue 5       w/o volt in eq. 1.5 ← Blue 6 Special versions.
w/ volt in eq. 1.5 ← Blue 3       w/o volt in eq. 1.5 ← Blue 4 Equation solvers only. LV ---> LV2 / ( γ + LV2) for all 3 terms in (1.1) and (1.2).
w/ volt in eq. 1.5 ← Blue 1       w/o volt in eq. 1.5 ← Blue 2 Special versions.
Simulators - 11/8/2012 Version Notes
  Version w/ volt in eq. 1.5.          Version w/o volt in eq. 1.5. Equation solvers only.
LV ---> LV / ( γ + LV) in the second term of (1.2).
Other features same as in the table below.
This is an incremental version.
Version S w/ volt in eq. 1.5.       Version S w/o volt in eq. 1.5. Special versions.
Version 3 w/ volt in eq. 1.5.       Version 3 w/o volt in eq. 1.5 (O1). Discussions on 10/15/2012.
Simulators - 10/26/2012 Version Notes
  Version w/ volt in eq. 1.5.          Version w/o volt in eq. 1.5. Equation solvers only. LV ---> LV / ( γ + LV) in (1.1) and the first term of (1.2).
Other features same as in the table below.
Version S w/ volt in eq. 1.5.       Version S w/o volt in eq. 1.5. Special versions.
Version 3 w/ volt in eq. 1.5.       Version 3 w/o volt in eq. 1.5. Discussions on 10/15/2012.

Simulators Version Notes
  Version w/ volt in eq. 1.5.          Version w/o volt in eq. 1.5. Equation solvers only. ↑↑↑ Modified
versions
Version S w/ volt in eq. 1.5.       Version S w/o volt in eq. 1.5. Special versions - Simplified version with equations 1.1 - 1.4 only. a23 and a32 are 0s.
Note: The first four equations are decoupled from the others. Please ignore parameters and variables that are not in 1.1 - 1.4.
Version 3 w/ volt in eq. 1.5.       Version 3 w/o volt in eq. 1.5. Discussions on 10/15/2012. New default ranges for beta1 and beta2. Version #s added in output files. Each simulation run will have its own folder with output files partitioned as requested. Limit TE3, IFNg, and TNFa exp. data to first 14 days. Step size is changeable.
Partitioned files are written in real time to prevent interruptions such as power outrage so they are not sorted. If a simulation finishes, the whole output file is sorted by variance.
  Version 2 w/ volt in eq. 1.5.       Version 2 w/o volt in eq. 1.5. Changes implemented from discussions on 10/3/2012. Major changes - kd1, kd2, and c1 become constants for each simulation. gp100 exp. data (Te3, C(t), IFNg, TNFa) added. Variance component option added.
Version w/ volt in eq. 1.5.          Version w/o volt in eq. 1.5. Major changes - New 0s. Variable step size. INFg excluded from fitness. Use C(t). New kp2 and kd4.
  Version 1.      Version ex1.      Version ex2. Archives. Last update - 9/20/2012. See hard copy version notes.

[1] David Klinke, "An Empirical Bayesian Approach for Model-based Inference of Cellular Signaling Networks", BMC Bioinformatics, November 2009.

Send me Suggestions about Design and Coding of the Simulator

Reference Simulators from Previous Projects