SBML-PET-MPI: A parallel parameter estimation tool for SBML based models Zhike Zi
BIOSS Centre for Biological Signalling Studies University of Freiburg, 79104 Freiburg, Germany http://www.bioss.uni-freiburg.de/cms/sbml-pet-mpi.html http://sites.google.com/site/sbmlpetmpi/ E-mail:
December 2, 2010 SBML-PET-MPI version 1.1
Copyright © 2010, Zhike Zi
Table of Contents 1 INTRODUCTION .............................................................................................................................. 1 1.1 OVERVIEW OF SBML-PET-MPI ........................................................................................................... 1 1.2 METHODS IMPLEMENTED IN SBML-PET-MPI ........................................................................................ 1 1.2.1 Parameter estimation with global optimization algorithm................................................... 1 1.2.2 Parameter uncertainty analysis with profile likelihood exploit algorithm ............................ 1 1.2.3 Parameter uncertainty analysis with bootstrap method ...................................................... 2 2 QUICK START .................................................................................................................................. 3 2.1 SYSTEM AND PACKAGE REQUIREMENTS .................................................................................................. 3 2.2.1 Windows................................................................................................................................ 3 2.1.2 Linux ...................................................................................................................................... 3 2.2.3 Mac OS .................................................................................................................................. 3 2.2 INSTALLATION AND START INSTRUCTIONS ................................................................................................ 3 2.2.1 Windows................................................................................................................................ 3 2.2.2 Linux ...................................................................................................................................... 4 2.2.2 Mac OS .................................................................................................................................. 4 3 USING SBML-PET-MPI..................................................................................................................... 5 3.1 RUN SBML-PET-MPI ........................................................................................................................ 5 3.1.1 MPI daemon (mpd) launch and exit ...................................................................................... 5 3.1.2 Start SBML-PET-MPI .............................................................................................................. 5 3.1.3 Options for ODE Solver .......................................................................................................... 6 3.2 EXPERIMENTAL DATA FILE FOR PARAMETER ESTIMATION ........................................................................... 7 3.2.1 Annotations for data File....................................................................................................... 7 3.2.2 Information about parameters to be estimated ................................................................... 7 3.2.3 Information about experimental data ................................................................................... 7 3.2.4 Information about constraints .............................................................................................. 9 3.2.5 An important note about the SD values in the data file ........................................................ 9 3.2.6 About trigger of events with time variable ......................................................................... 10 3.3 RESULTS DISPLAY AND SAVE ................................................................................................................ 11 3.3.1 Results displayed in the terminal and saved as files ........................................................... 11 3.3.2 Plot of data fitting and parameter analysis with MATLAB .................................................. 12 3.4 SPECIFIC EXPLANATION FOR DATA FROM MULTIPLE EXPERIMENTAL CONDITIONS ......................................... 14 4 EXAMPLES .................................................................................................................................... 15 4.1 A SIMPLE MODEL FOR ENZYME SUBSTRATE REACTIONS ........................................................................... 15 4.2 THE EPO MODEL (WITH REAL EXPERIMENTAL DATA) .............................................................................. 18 4.3 E. COLI TRYPTOPHAN OPERON MODEL (WITH DATA FROM DIFFERENT CONDITIONS) ................................... 20 5 FAQ .............................................................................................................................................. 23 6 REFERENCES ................................................................................................................................. 27
1 INTRODUCTION 1.1 Overview of SBML-PET-MPI SBML-PET-MPI is a parallel parameter estimation tool for Systems Biology Markup Language (SBML) (Hucka et al., 2003) based models. The tool allows the user to perform parameter estimation, parameter uncertainty and identifiability analysis by collectively fitting multiple experimental data sets. SBML-PET-MPI can run on Windows, Linux and Mac OS systems. The features of SBML-PET-MPI include:
SBML-PET-MPI supports model import and export in SBML format, a widely accepted standard for the exchange of biological models. All estimated parameters will be saved in a new SBML file, which can be imported to other SBML supported software packages.
SBML-PET-MPI estimates parameter values by fitting multiple experimental data sets. Experimental data can be produced under different conditions.
SBML-PET-MPI can implement uncertainty analysis of the estimated parameters with profile likelihood exploit algorithm or bootstrap method with synthetic data sets.
SBML-PET-MPI supports events that describe the discontinuous state changes in the model.
SBML-PET-MPI supports normalized data or mathematical expression of the data, for example, the sum of several variables in the model. It also supports the standard deviation for the data and the noise or measurement error existed in the experiments.
SBML-PET-MPI supports common mathematical expressions for the qualitative and quantitative description of the model, for example, the constraints for parameters or their combinations.
This manual will guide the users through SBML-PET-MPI by introducing various features and functions of it in more details.
1.2 Methods Implemented in SBML-PET-MPI SBML-PET-MPI uses Message Passing Interface (MPI) standard to parallelize the optimization and other algorithms. The following algorithms were implemented in SBML-PET-MPI. 1.2.1 Parameter estimation with global optimization algorithm SBML-PET-MPI uses stochastic ranking evolution strategy (SRES) for global optimization of the parameters. SRES is an evolutionary optimization algorithm that uses stochastic ranking as the constraint handling technique (Ji and Xu, 2006; Runarsson and Yao, 2000; Zi and Klipp, 2006). SBML-PET-MPI uses a similar principle of MPI implementation for SRES algorithm as libSRES (Ji and Xu, 2006). 1.2.2 Parameter uncertainty analysis with profile likelihood exploit algorithm SBML-PET-MPI parallelizes the algorithm of profile likelihood exploit for the identifiability analysis and confidence intervals analysis of the estimated parameters. Levenberg–Marquardt algorithm is used for the identifiability analayis of the parameters within profile likelihood exploit algorithm. Detailed information about the algorithm should refer to the publication by Raue et al. (Raue et al., 2009).
1
1.2.3 Parameter uncertainty analysis with bootstrap method SBML-PET-MPI implement another algorithm to get the confidence limits of the estimated parameter values is based on a bootstrap method, which generates M new synthetic data sets D j ( j 1, 2, ..., M ) from the actual experimental data sets s
D0 (Press, 1992): Dsj (i) D0 (i) i , j 1, 2, ..., M where
i
(E1.1)
is taken from a normal distribution N(0,i ) ,
i
corresponds to the
standard deviation in the i-th experimental data.
SBML-PET-MPI can run optimizations with new synthetic data sets. The new optimized parameter sets for the synthetic data sets are saved and can be used for calculating the distribution and standard deviations of the estimated parameter values. It is worth noting that the bootstrap method are time-consuming because for each synthetic data sets, SBML-PET-MPI needs to run a global optimization. Therefore, the running time will be proportional to the number of synthetic data sets.
2
2 Quick Start SBML-PET-MPI is available at http://sites.google.com/site/sbmlpetmpi/ or http://www.bioss.uni-freiburg.de/cms/sbml-pet-mpi.html. It can run on Windows, Linux and Mac OS system.
2.1 System and Package Requirements 2.2.1 Windows Install Cygwin and select the packages for gcc in "Devel" category (gcc-core, gccg++, gcc-g77 etc.). In addition, please select Python in "Interpreter" category, which is required for MPICH2 library. Cygwin install information can get from http://www.cygwin.com/. MPICH2 library is required for penalization. MATLAB is optional for result plotting. 2.1.2 Linux Gcc library (including, c, c++, g77/gfortran), Python and MPICH2 library. MATLAB is optional for result plotting. We recommend recent versions of gcc. 2.2.3 Mac OS Xcode, http://developer.apple.com/technologies/tools/xcode.html, free developer tools library from Apple Inc. MATLAB is optional for result plotting. SBML-PET-MPI supports Mac OS X10.5 and above systems with Intel CPUs.
2.2 Installation and Start Instructions Before install SBML-PET-MPI, the user needs to install MPICH2 library, which is an implementation of the Message Passing Interface (MPI) standard. Please follow the instructions in the website of MPICH2 for installing it. MPIC2 website: http://www.mcs.anl.gov/research/projects/mpich2/index.php After first installation of MPICH2, please create .mpd.conf file in the following way. cd $HOME touch .mpd.conf chmod 600 .mpd.conf and then use an editor to insert a line like MPD_SECRETWORD=mr45-j9z into the file. (Of course use some other secret word than mr45-j9z.) MPICH2-1.2.1p1 might have some problem in the installation on Cygwin. In this case, please install its previous release MPICH2-1.2.1 from the following link: http://www.mcs.anl.gov/research/projects/mpich2/downloads/tarballs/1.2.1/mpich2-1.2.1.tar.gz Note: Single version of SBML-PET-MPI does NOT require MPICH2. 2.2.1 Windows Step 1: Download the SBML-PET-MPI ZIP file (SBML-PET-MPI-v1.0-Cygwin.zip) and unzip the zip file to the target directory. Step 2. Open a Cygwin terminal and change to the directory including the unzipped SBML-PET-MPI files. Start mpd in the host machine with command "mpd &" if it is not started yet. Then, run SBML-PET-MPI.
3
2.2.2 Linux Step 1: Download the SBML-PET-MPI ZIP file (SBML-PET-MPI-v1.0-Linux64.zip for 64 bit system or SBML-PET-MPI-v1.0-Linux32.zip for 32 bit system) and unzip the zip file to the target directory. Step 2: Add SBML-PET-MPI lib sub-directory to the LD_LIBRARY_PATH variable. Please insert the following line at $HOME/.bashrc file (or $HOME/.bash_profile file). $SBML-PET_PATH is the full path for the directory including the unzipped SBML-PET-MPI files. export LD_LIBRARY_PATH=$LD_LIBRARY_PATH:$SBML-PET_PATH/lib Step 3: Open a terminal and change to the directory including the unzipped SBML-PET-MPI files. Start mpd in the host machine with command "mpd &" if it is not started yet. Then, run SBML-PET-MPI. 2.2.2 Mac OS Step 1: Download the SBML-PET-MPI ZIP file (SBML-PET-MPI-v1.0-Mac64.zip or SBML-PET-MPI-v1.0-Mac32.zip) and unzip the zip file to the target directory. Step 2: Add SBML-PET-MPI lib sub-directory to the DYLD_LIBRARY_PATH variable. Please insert the following line at $HOME/.bashrc file (or $HOME/.bash_profile file). $SBML-PET_PATH is the full path for the directory including the unzipped SBML-PET-MPI files. export DYLD_LIBRARY_PATH=$DYLD_LIBRARY_PATH:$SBML-PET_PATH/lib Step 3: Open a terminal and change to the directory including the unzipped SBML-PET-MPI files. Then, run SBML-PET-MPI. Note for automatic results plotting in MATLAB: if you cannot start MATLAB from your terminal, please add MATLAB to the PATH environment in Linux or Mac by adding the following line in the .bashrc or .bash_profile file at your home directory: export PATH=$PATH:MATLAB_PATH MATLAB_PATH is the executable MATLAB path, for example, the default of MATLAB R2010b in Mac: export PATH=$PATH:/Applications/MATLAB_R2010b.app/bin
4
3 Using SBML-PET-MPI 3.1 Run SBML-PET-MPI 3.1.1 MPI daemon (mpd) launch and exit For Linux and Windows, if no mpd (MPI daemon) is started, first launch the daemons by type mpd & The above command only need to run once in a terminal. To exit previous launched daemons, type mpdallexit In Mac OS system, the mpd is automatically started. Single version of SBML-PET-MPI does NOT require mpd start. 3.1.2 Start SBML-PET-MPI After mpd is launched, run SBML-PET-MPI by type mpirun -n cpu_number ./SBML-PET-MPI model_file data_file [options] or mpiexec -n cpu_number ./SBML-PET-MPI model_file data_file [options] or ./SBML-PET-MPI model_file data_file [options] (for the computer with 1 processor) Please note that "cpu_number" should be no less than 2 for parameter estimation in SBML-PET-MPI. Options: -v -h [--help] --with-ple --with-ple=### --with-sda --with-sda=### --with-simulation --set-ODEsolver --with-gen=### --with-lambda=### --with-miu=###
get SBML-PET-MPI version number display SBML-PET-MPI help information parameter analysis with profile likelihood exploit default number of sampling points is 200 parameter analysis with profile likelihood exploit and set number of sampling points=### (positive integer) parameter optimizations with synthetic data sets default number of synthetic data sets is 100 parameter optimizations with synthetic data sets and set number of synthetic data sets=### (positive integer) run simulation for the model set ODE Solver options set the number of evolutionary generations set the number of offspring population size set the number of parent population size
Examples: mpirun -n 4 ./SBML-PET-MPI Examples/ES.xml Examples/ES_data.txt mpirun -n 4 ./SBML-PET-MPI Examples/ES.xml Examples/ES_data.txt --with-gen=1000 mpirun -n 10 ./SBML-PET-MPI Examples/ES.xml Examples/ES_data.txt --with-ple mpirun -n 4 ./SBML-PET-MPI Examples/ES.xml Examples/ES_data.txt --with-ple=200 mpirun -n 4 ./SBML-PET-MPI Examples/ES.xml Examples/ES_data.txt --with-sda=50 ./SBML-PET-MPI Examples/ES.xml --with-simulation
5
Once the SBML-PET-MPI is started, a welcome information will appear (Figure 3.1).
Figure 3.1 Start of SBML-PET-MPI Tip #1: Please prepare the model in SBML format in advance with Copasi (Hoops et al., 2006) or CellDesigner (Funahashi, 2003). If the model contains events, we recommend CellDesigner. 3.1.3 Options for ODE Solver Currently, SBML-PET-MPI uses the ODE Solvers from ODEPACK to solve the ODE systems (Hindmarsh, 1983). The users can change the options of the ODE solvers with the option of "--set-ODEsolver", SBML-PET-MPI will display following prompts: ************************************************* > Setting ODEPACK Parameters (Default values): Relative Tolerance Parameter (scalar): 1e-08 Relative Absolute Parameter (scalar): 1e-08 Maximum Number of Steps allowed during one call: 10000 ************************************************* Please input Relative Tolerance Parameter (scalar) The user can specify the following three parameters for the ODE solver.
Relative Tolerance Parameter (RTOL). RTOL and Absolute tolerance parameter (ATOL) are used to control the estimated local error for ODE solver. The default value of ATOL in SBML-PET-MPI is 1e-8.
Absolute Tolerance Parameter (ATOL). The estimated local error in Y(i) will be controlled so as to be roughly less (in magnitude) than EWT(i) = RTOL×ABS(Y(i)) + ATOL The default value of RTOL in SBML-PET-MPI is 1e-8.
Maximum Number of Steps (MXSTEP). It is the maximum number of (internally defined) steps allowed during one call to the solver. The default value is 10000.
6
3.2 Experimental Data File for Parameter Estimation The user can prepare the data file by modify the standard "data_template.txt".
input file
Attention: DO NOT MODIFY THE RESERVED WORDS AND PHRASES (UPPERCASE) IN THE DATA FILE. The reserved words and phrases in the data file are: PARAMETER ID MAXIMUM VALUE EXPERIMENTAL CONDITIONS TIME COURSES NUM_TIME_POINTS MATH CONSTRAINTS
MINIMUM VALUE TOTAL NUMBER AT CONDITION TIME COURSES NUM_EXP_DATA DATA CONSTRAINT
The following part will explain the data file in details. 3.2.1 Annotations for data File The user can add annotation by using "#" at the beginning of the annotation line. For example, # Annotation lines should begin with # #*********************************************************** # Input data file for SBML-PET-MPI #*********************************************************** 3.2.2 Information about parameters to be estimated Copy the information about parameters to be estimated to part after PART I in the data file. This part should include "PARAMETER ID", "MINIMUM VALUE" and "MAXIMUM VALUE" information. For example, #*********************************************************** # PART I: PAREMETERS TO BE ESTIMATED # Copy your information for parameters to be estimated to the below #*********************************************************** PARAMETER ID MINIMUM VALUE MAXIMUM VALUE E 0.01 100 k1 0.001 1e3 k2 0.001 1e3 k3 0.001 1e3 Tip #2: The species ID shown in "Parameter ID" corresponds to the initial concentration/amount of the species. 3.2.3 Information about experimental data The information about experimental data is should put in the PART II as the following. #*********************************************************** # PART II: Experimental Data # Copy your information for experimental data # to the below #*********************************************************** TOTAL NUMBER of EXPERIMENTAL CONDITIONS = 1 # The number of Time Courses at each Experimental Condition AT CONDITION 1, The TOTAL NUMBER of TIME COURSES is 2 DATA OF TIME COURSE 1 AT CONDITION 1 7
NUM_TIME_POINTS 8 NUM_EXP_DATA 3 Time (Epo + dEpoe) SD Epo_EpoR 0 1987.83 12.58 NaN NaN NaN 5 1726.43 13.02 278.31 12.16 56.8 20 1479.82 5.28 265.88 4.21 300.23 60 1391.35 36.11 163.57 15.96 503.89 120 1312.71 13.75 136.44 10.04 603.32 180 1392.98 21.15 81.1 7.33 580.3 240 1501.47 31.2 59.93 9.13 494.45 300 1594.16 24.7 33.28 2.7 403.09 Normalized 0 0
SD NaN 1.17 2.73 19.28 22.98 13.34 27.05 24.49 0
(Epo_EpoRi + dEpoi)
SD
(1) Specify the total number of experimental conditions in line: TOTAL NUMBER of EXPERIMENTAL CONDITIONS = 1 (2) Specify the time course number and condition number for each time course data in line: AT CONDITION 1, The TOTAL NUMBER of TIME COURSES is 1 (3) Specify the number of time points and number of data for each time course data in line: NUM_TIME_POINTS 10 NUM_EXP_DATA 2 (4) Input the mathematical expression for each data after the above line: Time
(Epo + dEpoe)
SD
Epo_EpoR
SD
(Epo_EpoRi + dEpoi)
SD
Here, SD is for Standard Deviation σ (measurement error) or noisy of data. If you have the experimental data for the sum of several variables, you can specify the measurement with the math expression of the variables. For example, the above case, the first measured data is the sum of Epo and dEpoe. Therefore, the Mathematical Expression of the first column data is "Epo+dEpoe" Note: the variable parameters (not fixed, constant="false") should not appear in the mathematical expression of the data. If it should appear in the mathematical expression, please remove the variable parameter and redefine it as a species with boundaryCondition="true" constant="false" in the SBML model. Tip #3: An alternative to specify the data for the sum of several variables (or other mathematical expression of several variables) is to add an assignment rule in the SBML model and specify the data as the new species (NOT new parameters!) with assignment rule. Set the new added species with boundaryCondition="true" constant="false". Then add new assignment rule for the new species. This trick can also speed up the optimization process. For the above example, one can add new species "totalEpo" and "totaEpoi", with the following assignment rule: totalEpo = Epo + dEpoe totalEpoi = Epo_EpoRi + dEpoi Then the mathematical expression of the data Time
(Epo + dEpoe)
SD
Epo_EpoR
SD
(Epo_EpoRi + dEpoi)
SD
can be replaced with the following line in the new model Time
totalEpo
SD
Epo_EpoR
SD
totalEpoi
SD
Tip #4: Please always use the id (NOT the name) of the species and parameters in the data file. 8
(5) Insert the time course data after the line for mathematical expression of data. 0 1987.83 12.58 NaN NaN NaN NaN 5 1726.43 13.02 278.31 12.16 56.8 1.17 20 1479.82 5.28 265.88 4.21 300.23 2.73 60 1391.35 36.11 163.57 15.96 503.89 19.28 ... ... ... ... ... ... ... The time course data can be copied from Excel. If the Standard Deviation (SD) of some data potions is unknown, specify it as "NaN". SBML-PET-MPI will replace the NaN SD value with default value of 1. In all the cases, the SD column must be provided. Tip #5: If a measurement of some data at a certain time points is not known or they are outliers, you can specify the data at this particular point as "NaN". SBML-PET-MPI will skip the fit of the NaN data point. (6) Insert normalized information for the data. The number 0 indicates not normalized data and 1 means that the data normalized to the maximum value of the data in this time course. For example, Normalized
0
1
0
….
3.2.4 Information about constraints The information about constraints is provided in the last part of the data file. # The number of Constraints at each Experimental Condition AT CONDITION 1, The TOTAL NUMBER of CONSTRAINTS is 1 MATH OF CONSTRAINT 1 AT CONDITION 1, k1 - 10*k2 The math expression g(X, p) corresponds to the constraints with the inequality of
g ( X , p) 0
(E3.1)
where X is the involved initial conditions of the molecule amount/concentration and p is the parameter set. For the example “k1-10*k2” corresponds to "k110*k2≤0" and that is k1≤10*k2. If no constrains are imposed, for each condition of the experiments, specify the number of constraints is 0, for example, # The number of Constraints at each Experimental Condition AT CONDITION 1, The TOTAL NUMBER of CONSTRAINTS is 0 AT CONDITION 2, The TOTAL NUMBER of CONSTRAINTS is 0 3.2.5 An important note about the SD values in the data file In SBML-PET-MPI, SD values in the data file are used as weight factor of the data points for evaluating the chi-square objective function (E3.2).
2
totalNumEx pData
n
ynpred yndata SD
2
(E3.2)
If the experimental data sets are not normalized AND no standard deviation information is available. We recommended the user to set some weight factor as SD value for all the data points in the corresponding data set. For example, the user can set the SD value with the maximum value of this experimental data set. This trick can avoid the over-fitting of a certain data set when the values of several data sets are in different scale.
9
For example, if the user get the following data sets. Time 0 1 2 3 4 5 Normalized
dataset1 1000 2000 5000 2500 1200 1000 0
SD NaN NaN NaN NaN NaN NaN
dataset2 0.005 0.002 0.001 0.003 0.004 0.005 0
SD NaN NaN NaN NaN NaN NaN
In this case, if the use provide the data information with the above SD information. According to equation E3.2, all the SD weight factors will be the same for dataset1 and dataset2 (SD reset with value of 1 for NaN in SBML-PETMPI). Since the values in dataset1 are much larger than those in dataset2, the fitting of dataset1 will overtake the fitting of datast2. Therefore, the user should specify SD information as weight factors to avoid the over-fitting of dataset1. The user can choose the maximum value of each data set as their weight factors (SD). Alternatively, the user can assume 5% or 10% of the data points as the weight factor (SD). The following modified data set information will increase the collective fitting of all the data sets. Time 0 1 2 3 4 5 Normalized
dataset1 1000 2000 5000 2500 1200 1000 0
SD 5000 5000 5000 5000 5000 5000
dataset2 0.005 0.002 0.001 0.003 0.004 0.005 0
SD 0.005 0.005 0.005 0.005 0.005 0.005
3.2.6 About trigger of events with time variable For time dependent events, the user should add the trigger time into the time course of the experimental data if the trigger time doesn't appear in the time course. For example, if the model has the event that at time point 60, the signal cue species L is washed out, which corresponds to the following SBML event definition.
10
Suppose the user has the following data set. Time 0 10 45 90 120 Normalized
data1 0.02 0.43 0.50 1.00 0.20 0
SD 0.02 0.05 0.10 0.10 0.05
Since the above data set doesn't include the event trigger time 60, the user can add a time point in the time course. Time 0 10 45 60 90 120 Normalized
data1 0.02 0.43 0.50 NaN 1.00 0.20 0
SD 0.02 0.05 0.10 NaN 0.10 0.05
For other trigger of events without time variable, SBML-PET-MPI will automatically detect the trigger time of the events.
3.3 Results display and save 3.3.1 Results displayed in the terminal and saved as files SBML-PET-MPI displays the result of optimization after each generation in the terminal as the following: -------------------- MESSAGE FROM SBML-PET-MPI ------------------->>> PARAMETER ESTIMATION with Experimental Data...... > generation: 1, best result from generation 1 > chisquare value = 2.874709e+04, penalty value for constraints = 0.000000 > best parameters = 7.208e-05 4.626e-05 1.698e-01 1.037e-04 8.938e+01 2.652e+02 1.941e+03 1.868e-06 1.498e+00 1.007e+03 > used time = 2.215 seconds > generation: 2, best result from generation 2 > chisquare value = 2.270294e+04, penalty value for constraints = 0.000000 > best parameters = 1.002e-04 2.092e-05 2.295e-01 7.647e-02 2.347e-02 7.488e-04 2.255e+03 4.425e-07 1.951e-06 9.985e+02 > used time = 2.496 seconds > generation: 3, best result from generation 3 > chisquare value = 1.221319e+04, penalty value for constraints = 0.000000 > best parameters = 4.048e-01 2.423e-04 1.325e-01 4.580e-04 2.467e+01 3.787e+01 1.940e+03 2.634e-06 2.700e+00 1.003e+03 > used time = 2.792 seconds ... ... When the running of SBML-PET-MPI is finished, the SBML models with the parameter set from the best fit, the optimization history and the summary of parameter estimation analysis are saved in "result" subdirectory of the installed SBML-PET-MPI directory (Figure 3.2).
11
Figure 3.2 Final results of parameter estimation The estimated parameters are stored in SBML files: "BestFitnessSBML.xml" (SBML format level 2 version 4) and "BestFitnessSBML_L2V1.xml" (SBML format level 2 version 1), which can be imported to other SBML software for analysis. Even SBML-PET-MPI is suddenly stopped or aborted, an updated SBML model with the latest best fitted parameters will be automatically saved as well. The summary of the estimated parameter values are saved in the file of "ParameterAnalysis.txt", which looks like the following: ******************************************************************** Summary of Estimated Parameters with Profile Likelihood Exploit ******************************************************************** chisquare = 2.11705 total number of experimental data, Nd = 38 total number of estimated parameters, Np = 4 chisquare/Nd = 0.055712 ID (Name) Value from Best Fit 95% confidence interval E 5.007936e+00 (4.970768e+00; 5.045104e+00) k1 1.977730e+00 (1.896612e+00; 2.066574e+00) k2 9.822647e-01 (9.477320e-01; 1.019867e+00) k3 9.985696e-01 (9.943764e-01; 1.002763e+00) 3.3.2 Plot of data fitting and parameter analysis with MATLAB SBML-PET-MPI provides the plot of data fitting and parameter analysis from profile likelihood with MATLAB. If MATLAB is installed in the host computer, the summary of the data fitting, optimization history and the profile likelihood of the parameters will be automatically plotted at the end of parameter estimation analysis (Figure 3.3 and Figure 3.4). The user can also copy the "result" subdirectory from the host to the local computer with MATLAB, then run plotResults (included in SBML-PET-MPI) in the MATLAB of the local computer. Please note that the "result" directory should put as a subdirectory of the directory where the file "plotResults.m" locates.
12
Figure 3.3 Summary of data fitting plots from SBML-PET-MPI
Figure 3.4 Profile likelihood plots of parameters from SBML-PET-MPI. The gray dotted line in denotes the threshold for 95% confidence intervals of the parameters calculated with profile likelihood exploit algorithm. The red circle point refers to the parameter value from the best fit. Tip #6: Do NOT modify the files saved in result directory. The change of
13
these files might affect the plot of data fitting and parameter analysis with MATLAB.
3.4 Specific Explanation for Data from Multiple Experimental Conditions If the experimental data are obtained from different experimental conditions, the user need modify some parameters defined in the model to distinguish the difference of the conditions. After SBML-PET-MPI started, a file named "ExpConditionsData.txt" will be produced at the "temp" sub-directory of SBML-PET-MPI. For different experimental conditions, the possible varied parameters values will be listed in this file. The user should modify the corresponding data for different conditions and save the file, then press ENTER key to continue. The following messages will display on the terminal: -------------------- MESSAGE FROM SBML-PET-MPI -------------------Please open the file temp/ExpConditionsData.txt and modify the data for different experimental conditions > If you have finished, press ENTER key to continue The value of the species listed in the file of "ExpConditionsData.txt" indicates the corresponding initial concentration/amount of the species. The user can also define some parameters to specify different experimental conditions, define them as global parameters rather than local parameters in the model.
14
4 Examples This section will present some examples showing the application and usage of SBML-PET-MPI. These examples cover different applications of SBML-PET-MPI.
4.1 A Simple Model for Enzyme Substrate Reactions The well-known Michaelis-Menten equations describe enzyme kinetics, which include four molecular species, namely the enzyme, E, the substrate, S, the product, P and the intermediate, ES. The reactions are listed below: (E4.1)
k1 k3 E S ES E P k2 The ODE system consists of 4 ordinary differential equations.
d[ E] k 2 [ ES ] k3 [ ES ] k1[ E ][ S ] dt
(E4.2)
d[S ] k2 [ ES ] k1[ E ][ S ] dt
(E4.3)
d [ ES ] k1[ E ][ S ] k2 [ ES ] k3[ ES ] dt
(E4.4)
d [ P] k3 [ ES ] dt
(E4.5)
To demonstrate the usage of SBML-PET-MPI, we assume that the true values of k1 = 2, k2 = 1, k3 = 1 and the initial concentration of [E] = 5 [S] = 10, [ES] = 0, [P] = 0. Then we run SBML-PET-MPI to estimate the values of k1, k2, k3 and the initial concentration of E and S within the following range Table 4.1 Range of parameters to be estimated Parameter ID
Minimum Value
Maximum Value
E
0.01
100
k1
0.001
1000
k2
0.001
1000
k3
0.001
1000
There are 4 parameters to be estimated. The initial concentration of E covers four magnitudes, while the value of kinetic parameters k1, k2 and k3 covers 6 magnitudes. The in silico experimental data is generated in the following way: we randomly sampled 20 experimental data sets from a normal distribution: the mean = the original simulation data and standard deviation = 5% corresponding data. The mean and corresponding standard deviations from these 20 artificial data sets are
15
used for parameter estimation. According the description in section "3.2 Experimental Data File for Parameter Estimation", the data file for this example is # data file for Michaelis-Menten Equations Model #*********************************************************** # PART I: PAREMETERS TO BE ESTIMATED #*********************************************************** PARAMETER ID MINIMUM VALUE MAXIMUM VALUE E 0.01 100 k1 0.001 1e3 k2 0.001 1e3 k3 0.001 1e3 #*********************************************************** # PART II: Experimental Data #*********************************************************** TOTAL NUMBER of EXPERIMENTAL CONDITIONS = 1 # The number of Time Courses at each Experimental Condition AT CONDITION 1, The TOTAL NUMBER of TIME COURSES is 1 DATA OF TIME COURSE 1 AT CONDITION 1 NUM_TIME_POINTS 19 NUM_EXP_DATA 2 Time S SD P SD 0 9.971476235 0.453141397 0 0 0.2 5.292235504 0.268280699 0.628160073 0.019133061 0.4 4.39726553 0.208274125 1.458557064 0.070669452 0.6 3.726573609 0.165386153 2.232857326 0.084280801 0.8 3.080806333 0.134350275 3.09237209 0.16668692 1 2.510472349 0.141611198 3.80900556 0.173941468 1.2 1.982882059 0.067026255 4.557345262 0.195609836 1.4 1.536801613 0.076873175 5.215831005 0.280492893 1.6 1.153292228 0.081686909 5.902796165 0.275531448 1.8 0.846651717 0.038701391 6.465823544 0.358515561 2 0.65483179 0.047542333 6.891697826 0.446185162 3 0.179484392 0.009625662 8.74544011 0.320153952 4 0.061421876 0.003239136 9.421806752 0.546102974 5 0.023203428 0.001171091 9.641262587 0.454142332 6 0.008947044 0.000525554 9.711707321 0.493301063 7 0.003764036 0.00012968 9.884865955 0.512322358 8 0.001521486 9.03E-05 9.941232427 0.46427798 9 0.000603876 2.71E-05 9.839203148 0.44095379 10 0.000251923 1.30E-05 10.06626382 0.314179785 Normalized 0 0 #*********************************************************** # PART III: Constraints Information #*********************************************************** # The number of Constraints at each Experimental Condition AT CONDITION 1, The TOTAL NUMBER of CONSTRAINTS is 0 Figure 4.1 Parameter estimation data file for the enzyme substrate model The fitting of the data and the profile likelihood of the parameters are shown in Figure 4.2 and Figure 4.3.
16
Figure 4.2 Summary of data fitting plots for the enzyme substrate model
Figure 4.3 Profile likelihood plots of parameters for the enzyme substrate model. The gray dotted line in denotes the threshold for 95% confidence intervals of the parameters calculated with profile likelihood exploit algorithm. The red circle point refers to the parameter value from the best fit.
17
4.2 The Epo Model (with Real Experimental Data) Becker et al. developed a mathematical model for the erythropoietin (Epo) and Epo receptor (EpoR) interaction (Becker et al., 2010). This model was calibrated with multiple real experimental data sets. Here we use the model and the experimental data sets in this work as a benchmark example for SBML-PET-MPI. Since the data sets in this work are generated for two models (the core model and the auxiliary model), we merged the original two models in one model ("Examples/Epo.xml") in order to simultaneously fitting the data sets for the two models. The SBML models and data sets are obtained from the following website: http://webber.physik.uni-freiburg.de/~jeti/Science_Becker_data_models/ In SBML-PET-MPI, we set the same estimated parameters with the same ranges and use the same data sets as the original work. Detailed information is described in the file of "Examples/Epo_data.txt". As shown in Table 4.2, SBML-PET-MPI can reproduce a similar parameter estimation result as original work. It also obtained similar parameter intervals information. Table 4.2 Values and statistical interval of the estimated parameters from SBML-PETMPI and those from the original study Parameter ID
Estimated value from SBML-PET-MPI
Estimated value in Becker et al.
95% confidence interval from SBML-PET-MPI
95% confidence interval in Becker et al.
kt
0.03270
0.03294
(0.03259;0.03278)
(0.0300; 0.0365)
kon
1.0523×10-4
1.0496×10-4
(1.049×10-4; 1.055×10-4)
(1.003×10-4;1.097×10-4)
ke
0.07465
0.07483
(0.0744; 0.0749)
(0.0723; 0.0776)
kex
0.00982
0.00994
(0.00979; 0.0985)
(0.0082; 0.0119)
kdi
0.003150
0.003179
(0.00314; 0.00316)
(0.00272; 0.00365)
kde
0.01635
0.01640
(0.01630; 0.01640)
(0.01557; 0.01726)
Epo
2030.1
2030.19
(2023.79; 2036.47)
(2024.98; 2035.41)
kon_SAv
2.287×10-6
2.294×10-6
(2.280×10-6; 2.294×10-6)
(2.162×10-6; 2.430×10-6)
kex_SAv
0.0108
0.0110
(0.0107; 0.0109)
(0.0041; 0.0186)
SAv
999.291
999.293
(996.168; 1002.413)
(999.173; 999.413)
The fitting of all the data sets is shown in Figure 4.4. In this case, the SBMLPET-MPI was run with 7 CPUs and the chi-square objective function converged within one minute.
18
Figure 4.4 Data fitting and convergence curve of Epo model To evaluate the performance of the parallelized parameter estimation analysis algorithms, we recorded the running time of SBML-PET-MPI for global optimization (parameter estimation) and profile likelihood exploit analysis (parameter identifiability analysis) of the Epo model. The speed up scalability is good when the number of processors is up to 10, but the communication time between the processors decreases the speed up performance for this small example (Figure 4.5). The speed up scalability is in general better for complex models with more estimated parameters.
Figure 4.5 The running time and speed up for the optimization (1000 generations) and parameter identifiability analysis of the Epo model with real experimental data sets using 1-20 CPUs. The running time might slightly vary at different execution of the program.
19
4.3 E. Coli Tryptophan Operon Model (With Data from Different Conditions) Sharad Bhartiya et al. developed a mathematical model to study the effect of external tryptophan on the trp operon (Bhartiya et al., 2003). The model accounts for the effect of feedback repression by tryptophan with a Hill equation. The model describes such a process: (1) In the first step, tryptophan molecule, T, and the aporepressor molecule, R, form the intermediate, RT. (2) Then, another tryptophan molecule binds with RT yielding the holorepressor, RT2. (3) The holorepressor binds with the free operator, O, and forms operon–holorepressor complex, ORT2, which represses tryptophan synthesis. A schematic representation of the trp operon system is provided in Figure 4.6.
Figure 4.6 Schematic representation of Dynamic Model of E. coli Tryptophan Operon. Tryptophan concentration is influenced by (a) enzyme synthesis (E) with kinetic constant k1, (b) enzyme catalyzed reaction for tryptophan synthesis from a nitrogenous substrate (NS) with kinetic constant kd, and (c) instantaneous uptake from the environment. Sufficient availability of tryptophan, T, leads to binding with the aporepressor molecule R, with a dissociation constant K1. The resulting holorepressor next binds with the free operon, O, with a dissociation constant K2, resulting in transcriptional and translational repression of enzyme synthesis. (Figure is modified from Fig. 1 in Bhartiya et al., 2003) The ODE system for this model includes two species and 11 parameters. Detailed information about this model is described in Bhartiya et al., 2003. The SBML file provided in SBML-PET-MPI is modified from the model in JWS Online models (http://jjj.biochem.sun.ac.za/index.html). The data produced with different extracellular trypophan concentrations (0.01 μM, 0.1 μM, 1.0 μM and 10 μM) are used to estimate 8 parameters values: ki1, Ot, eval, fval, Tomax, gval, kg, mu. The necessary files for this problem are located in directory of Examples with names of "Ecoli.xml" and "Ecoli_data.txt". Note: After SBML-PET-MPI is started, please modify the corresponding data in the file "ExpConditionsData.txt" at the "temp" sub-directory of SBML-PET-MPI. Set the values of parameter "Text" to be 0.01, 0.1, 1.0 and 10 at condition 1, condition 2, condition 3 and condition 4 respectively, which indicates that the extra-cellular trypophan concentrations are 0.01 0.1, 1.0 and 10μM for different conditions. Don’t change other parameters values listed in the file "ExpConditionsData.txt". After modifying the experimental conditions file, press 20
ENTER key to continue. The modified "ExpConditionsData.txt" has the following contents: # DATA FOR EXPERIMENTAL CONDITION 1 DATA ID VALUE Enz 0 Ts 0 s 0 Text 0.01 # DATA FOR EXPERIMENTAL CONDITION 2 DATA ID VALUE Enz 0 Ts 0 s 0 Text 0.1 # DATA FOR EXPERIMENTAL CONDITION 3 DATA ID VALUE Enz 0 Ts 0 s 0 Text 1 # DATA FOR EXPERIMENTAL CONDITION 4 DATA ID VALUE Enz 0 Ts 0 s 0 Text 10 The fitting of the data and the profile likelihood of the parameters are shown in Figure 4.7 and Figure 4.8.
Figure 4.7 Result of Best Fit for E. coli. Tryptophan Operon Model
21
Figure 4.8 Profile likelihood plots of parameters for E. coli. Tryptophan Operon Model. The gray dotted line in denotes the threshold for 95% confidence intervals of the parameters calculated with profile likelihood exploit algorithm. The red circle point refers to the parameter value from the best fit.
22
5 FAQ Q1: What is the objective function ( χ2) defined in SBML-PET-MPI? A: SBML-PET-MPI defines the objective/cost function (chi-square, χ2) with the sum of the weighted least squares between model simulation data and the experimental data. The SD values are the weight factors, which can be set as the measured standard deviation or the maximum value of the data sets. More details are described in section of "3.2.5 An important note about the SD values in the data file". Q2: How to prepare SBML file for my model? A: There are many tools for create the model in SBML format, for example, Copasi or CellDesigner. If the model contains events, we recommend CellDesigner for defining events. Q3: How to prepare the data file for parameter estimation? A: Read the section of "3.4 Prepare the Data File for Parameter Estimation" in manual documentation. Q4: How to define events in the SBML-PET-MPI? A: We recommend the user to use CellDesigner to edit the model with events. If time will appear in the events, use it as "time". The user does not need set a parameter/species for time, it will be automatically recognized by SBML-PET-MPI. For the time-dependent events, please pay attention to the notes for the data file preparation at section "3.2.6 About trigger of events with time variable". For other variables that are controlled by events, the user should define them as global parameters or species variables. Q5: How much of memory needed for running SBML-PET-MPI? A: SBML-PET-MPI use dynamic memory allocation and free strategy. Therefore, the usage of memory in SBML-PET-MPI depends on the number of processors (CPUs) started, the complexity of the model, the number of data points and parameters to be used. We recommend 100 M bytes of memory per processor. Q6: What does it mean if the estimated parameters are marked as nonidentifiable with profile likelihood exploit analysis? A: In this case, you should check the plot of profile likelihood of the parameters. If the plot of profile likelihood of the parameter is concave, then you can increase number of sampling points for profile likelihood exploit (--with-ple=###). If the plot of profile likelihood of the parameter is concave and it hits the bound of parameters (min or max), then you can try to increase the range of the corresponding parameter. If the plot of profile likelihood of the parameter is flat or convex, it means that the parameter value cannot be identified with the data sets. The user may either get more informative data set or measure (or fix) some estimated parameter values from other data sets. In addition, if the number of SRES generations is too small, this may also happen. In this case, please increase the number of SRES generations for global optimization. Q7: What about the speed up performance of SBML-PET-MPI and how fast is SBML-PET-MPI compared with other tools? A: SBML-PET-MPI has good speedup scalability with the increasing number of processors. SBML-PET-MPI speeds up the process of prameter estimation and identifiability analysis when the number of precoessors is more than 3. It would
23
be no apparent speed up between 1 processor and 2 processors due to the necesary model preocssing, results analysis designed in SBML-PET-MPI. Beacuse libSRES uses a different ODE solver (CVODE in C) from SBML-PET-MPI (ODEPACK in Fortran) and it doesn't need to decode the SBML model information. One should expect that after manually translating the specific model into ODE system, manually coding the estimated parameters and data information in libSRES, libSRES would be faster than SBML supported parameter estimation tools because our tool and other SBML tools need to process the SBML model, map the experimental data to the model simulation and update the model with the fitted parameters (write updated sbml model files during optimization). Therefore, it would be difficult to make a fair comparison of the absolute optimization time between compiled parameter estimation program with other SBML tools. Due to the above reason, it would be more reasonable to compare the speedup performance (not the absolute optimization time) of our tool to libSRES. We tested them with enzyme substrate reactions model at the same computer (same estimated parameters, generation, miu, lamada parameters). The speed up performance of both tools are similar (Figure 5.1). This is not a surprising result because our tool and libSRES used the same SRES algorithm and both are parallelized with MPI protocol.
Figure 5.1 Speed up performance of SBML-PET-MPI and libSRES Different SBML software tools might have different ODE solver, encoded in different programming language and use different optimization algorithms for parameter estimation. It would be difficult to compare the speed of optimization across different tools. Since both SBML-PET-MPI and COPASI are encoded in C. They have the same ODE solver from ODEPACK and also can use SRES for parameter estimation. We tested COPASI and SBML-PET-MPI with the enzyme substrate reactions model (details see Section 4.1). The estimated parameters and their range are set the same in COPASI as those in SBML-PET-MPI (details in Section 4.1). It takes about 14-15 minutes to finish the optimization with the following SRES algorithm setting: Number of Generation = 1000, Population Size = 238. With the same example running in the same computer system, we also implemented the optimization in SBML-PET-MPI with the same SRES algorithm setting as those set in COPASI. As shown in Table 5.1, SBML-PET-MPI has better performance than COPASI with 1 processor and it has good speed up scalability with the increasing number of processors. 24
Table 5.1 Performance of SBML-PET-MPI and its comparison to COPASI Number of Processors
Running time in COPASI 4.6
Running time in SBML-PET-MPI
1
890 seconds
370 seconds
2
NA
344 seconds
3
NA
188 seconds
4
NA
142 seconds
5
NA
120 seconds
6
NA
106 seconds
Note: The real running time might be slightly different at different execution of the programs
Q8: Can I run SBML-PET-MPI in computer with one processor? A: Yes, please download the single version of SBML-PET-MPI for running with one processor. In this case, there would be no speed up. The number of processors for SBML-PET-MPI should be more than 3 for speed up. Q9: What are the principles used to divide the global optimization, the profile likelihood method and the bootstrap method among multiple processors? A: The SRES global optimization process employs the classic (μ, λ)-ES evolution strategy algorithm, in which the selection is taken from the λ offspring only, whereas their μ parents are ignored even the fitness of the parents are better than that of the new generation. For each geneartion, SBML-PET-MPI distribute the evluations of the objective functions for λ offsprings into other n-1 processors with MPI protocol. The first processor is used to collecting the results from other processors and coordinate the overal algorithm. For the profile likelihood method, since the forward and backward profile likelihood can be independently calculated from the starting point of the best fitted parameter value. For m number of parameters, there are 2 × m profile likelihood calculations (forward and backward). Therefore, we distributed 2 × m profile likelihood calculations into n-1 processors. The first processor is used for result updating and the coordination of the whole algorithms. The bootstrap method with new synthetic data sets is based on the principle of parallel distribution as the SRES global optimization. The difference is that the bootstrap method requires many different runs of optimizations with new synthetic data sets. Q10: Why I get this message: "./SBML-PET-MPI: error while loading shared libraries: libsbml.so.4: cannot open shared object file: No such file or directory" A:The message comes out because the dynamic-link library (DLL) was not loaded in your system. The necessary DLL libraries are located at lib subdirectory of SBML-PET-MPI, to solve it, there are two solutions: (1) Permanent solution: you can export the LD_LIBRARY_PATH (for Linux system) or DYLD_LIBRARY_PATH (for Mac system) in the .bashrc or .bashr_profile file. Details please see the installation instruction notes in Page 4. (2) Temporary solution: After start a new terminal, first change to SBML-PET-MPI directory and then input the following command: export LD_LIBRARY_PATH=./lib
(for Linux system)
export DYLD_LIBRARY_PATH=./lib
(for Mac system)
25
Please note that solution 2 should be repeated for every new started terminal. If you stays in the same terminal, "export LD_LIBRARY_PATH=./lib" only need to do one time. Q11: How can I set the parameters for SRES global optimization algorithm? A: By default, SBML-PET-MPI uses an adaptive method for setting SRES algorithm parameters: lambda (offspring population size) and miu (parent population size), which depend on the number of parameters to be estimated: lambda = NumParameterToBeEstimated*7+210 miu = NumParameterToBeEstimated+30 The default maximum number of evolutionary generations is 2000. The uses can specify these parameter values by the following options: --with-gen=### --with-lambda=### --with-miu=###
set the number of evolutionary generations set the number of offspring population size set the number of parent population size
The ratio between parents and offspring population size is recommended to be about 1/3~1/10. According to the previous experiences reported, the ratio around 1/7 is good for most of the problems. Q12: More questions about SBML-PET-MPI? A: Please contact Zhike Zi by email:
.
26
6 References Becker, V., Schilling, M., Bachmann, J., Baumann, U., Raue, A., Maiwald, T., Timmer, J., and Klingmuller, U. (2010). Covering a broad dynamic range: information processing at the erythropoietin receptor. Science 328, 1404-1408. Bhartiya, S., Rawool, S., and Venkatesh, K.V. (2003). Dynamic model of Escherichia coli tryptophan operon shows an optimal structural design. Eur J Biochem 270, 26442651. Funahashi, A., Tanimura, N., Morohashi, M., and Kitano, H (2003). CellDesigner: a process diagram editor for gene-regulatory and biochemical networks. BIOSILICO 1, 159. Hindmarsh, A.C. (1983). ODEPACK, A Systematized Collection of ODE Solvers. In Scientific Computing, R.S. Stepleman, M. Carver, R. Peskin, W.F. Ames, and R. Vichnevetsky, eds. (Amsterdam, North-Holland), pp. 55-64. Hoops, S., Sahle, S., Gauges, R., Lee, C., Pahle, J., Simus, N., Singhal, M., Xu, L., Mendes, P., and Kummer, U. (2006). COPASI--a COmplex PAthway SImulator. Bioinformatics 22, 3067-3074. Hucka, M., Finney, A., Sauro, H.M., Bolouri, H., Doyle, J.C., Kitano, H., Arkin, A.P., Bornstein, B.J., Bray, D., Cornish-Bowden, A., et al. (2003). The systems biology markup language (SBML): a medium for representation and exchange of biochemical network models. Bioinformatics 19, 524-531. Ji, X., and Xu, Y. (2006). libSRES: a C library for stochastic ranking evolution strategy for parameter estimation. Bioinformatics 22, 124-126. Press, W.H. (1992). Numerical recipes in C : the art of scientific computing, 2nd ed. edn (Cambridge University Press). Raue, A., Kreutz, C., Maiwald, T., Bachmann, J., Schilling, M., Klingmuller, U., and Timmer, J. (2009). Structural and practical identifiability analysis of partially observed dynamical models by exploiting the profile likelihood. Bioinformatics 25, 1923-1929. Runarsson, T.P., and Yao, X. (2000). Stochastic ranking for constrained evolutionary optimization. Ieee Transactions on Evolutionary Computation 4, 284-294. Zi, Z., and Klipp, E. (2006). SBML-PET: a Systems Biology Markup Language-based parameter estimation tool. Bioinformatics 22, 2704-2705.
27
