Raia_CancerResearch2011
Reference: Raia et al. Dynamic Mathematical Modeling of IL13-induced Signaling in Hodgkin and Primary Mediastinal B-cell Lymphoma Allows Prediction of Therapeutic Targets. Cancer Research 71, pp. 693-704, 2011.
Folder: /Examples/Raia_CancerResearch2011.
Facts: The model contains 205 data points, 39 free parameters and 4 experimental conditions.
This page provides a guided tour through the D2D software using the dynamical model of JAK2-STAT5 signal transduction published in Raia et al., Cancer Research 2011. All required definition files are located in /Examples/Raia_CancerResearch2011, and their respective subfolders.
Like most projects in D2D this example is based on the following.
- a model definition file specifying inputs, reactions, observations, and error model
- data definition files depicting experimental setups different from those defined previously
- data files containing the measured data
- a setup file to initialize the software, load model as well as data and create the associated C files
D2D works hierarchical in the sense that specifications in the model definition are considered as defaults but may be overridden in a particular data definition. If on the other hand a data definition file is empty except the mandatory keywords, the defaults from the model definition are used. The latter case is the situation in this example.
The following tasks are routinely performed after model and data have been loaded.
- Plotting to visualize the dynamics generated from the current set of parameters
- Fitting to adjust the parameters to describe the data as good as possible
- Uncertainty analysis to assess parameter identifiability and calculate confidence intervals
The file il13_jak2_stat5.def in the /models subfolder specifies the ODE model, i.e. names and units of the dynamic variables and their interplay by providing reaction schemes. A model definition file is divided in several sections that are described in detail on Setting up models. Most prominently the differential equation system is defined in the REACTIONS part as conversion scheme.
#!matlab
...
REACTIONS
Rec -> IL13_Rec CUSTOM "il13_level * Kon_IL13Rec * 2.265 * Rec"
Rec -> Rec_i CUSTOM "Rec_intern * Rec"
Rec_i -> Rec CUSTOM "Rec_recycle * Rec_i"
...
Here, the membrane associated receptor either forms a complex with the ligand IL13 (first line) or is internalized (second line). The internalized receptor is again recruited to the membrane (third line). Each reaction rate is relative to the involved dynamic variables. The D2D software parses the reaction rate specification and considers each term that is not a dynamic variable as parameter, which later can be estimated from the data. As il13_level is set to the input dose for each experiment the introduced parameters from the reactions above are Kon_IL13Rec, Rec_intern and Rec_recycle.
Later in the same file il13_jak2_stat5.def in the /models subfolder the observation functions are specified to make the connection between dynamic variables and experimental data. In particular the observation functions are defined in the OBSERVABLES part and the associated error model after the keyword ERRORS.
#!matlab
...
OBSERVABLES
RecSurf_obs C au conc. 0 0 "Rec + IL13_Rec + p_IL13_Rec"
IL13_cell_obs C au conc. 0 0 "scale_IL13_cell_obs * IL13_cell"
...
ERRORS
RecSurf_obs "RecSurf_obs1 + RecSurf_obs2*RecSurf_obs"
IL13_cell_obs "IL13_cell_obs1 + IL13_cell_obs2*IL13_cell_obs"
...
In these exemplary lines the total amount of receptor on the surface is defined as the observable RecSurf_obs. As the IL13 concentration can only be observed on a relative scale, the observational parameter scale_IL13_cell_obs is introduced. For each observable it is required to specify an error model, in this case a relative error RecSurf_obs2*RecSurf_obs with an offset RecSurf_obs1 is implemented.
The file MedB1_real_data.def in the /data subfolder is empty except for keywords hence the observables defined in the model definition above are not altered.
The associated file MedB1_real_data.xls contains the experimental data.
#!
time il13_level pIL4Ra_obs pJAK2_obs IL13_cell_obs RecSurf_obs ...
0.0 0 NaN NaN 0.616 1.289752779
5.0 0 NaN NaN NaN NaN
10.0 0 NaN NaN NaN 1.57854438
...
Note that the headers of this file must match either the independent variable (here: time), a model parameter (e.g. il13_level), or the name of an observable (e.g. RecSurf_obs). Columns with unrecognized headers are ignored.
Executing the script Setup.m loads model and data in MATLAB and automatically generates the C-files.
After the setup script finished successfully full functionality can be used. All information about model and data are now available in the global variable ar. Below some basic function are explained briefly, for further information see [Extended Functionalities](Extended Functionalities) and [Function Reference](Function Reference).
Calling arPrint shows the list of all model parameters, their current values, upper bound (ub) and lower bound (lb), ... in the MATLAB command window.
#!matlab
Parameters: # = free, C = constant, D = dynamic, I = initial value, E = error model
name lb value ub 10^value fitted prior
# 1|D | CD274mRNA_production | -5 -1.7 +3 | 1 +0.021 | 1 | uniform(-5,3)
# 2|D | DecoyR_binding | -5 -2.4 +3 | 1 +0.0039 | 1 | uniform(-5,3)
# 3|D | JAK2_p_inhibition | -5 -1.1 +3 | 1 +0.078 | 1 | uniform(-5,3)
...
The command arPlot generates one figure for each data-set (one in total for this example).
By default the observations are displayed. A way to visualize also the underlying ODE solutions is to invoke arPlotter which opens a small graphical user interface that allows selecting data sets and available quantities to display by arPlot command.
After activating the checkbox PlotX and hitting the button Update Plots a new figure is opened showing the dynamics of the internal variables.
arFit runs the currently selected numerical optimization algorithm to calibrate the model to the experimental data. Maximum likelihood estimation is used here. If the the option flag in the global variable ar.config.showFitting = true the process of the model calibration is shown during the estimation (see arPlot). In this example the optimization routine is already converged to a local optimum hence calling arFit again does improve the likelihood further.
To check whether the current local optimum is also globally optimal a sequence of n fits is executed using the command arFitLHS(n). For deterministic optimization the initial parameter guesses are generated by Latin hypercube sampling (LHS). The results of the n fits are displayed by arPlotChi2s sorted by goodness of fit and the parameter values of the best fit replace the previous values.
In this figure one can nicely see 1) the improvement of the initial likelihood value by several orders of magnitude, and 2) the plateaus after optimization demonstrating the robustness of the results.
ple, plePlot, arPLETrajectories, Export figures