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.
Like most projects in D2D this example is based on
- 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
These definition, data and setup files are located in /Examples/Raia_CancerResearch2011, and their respective subfolders.
D2D works hierarchical in the sense that specifications in the model definition are considered as defaults but may be overwritten for particular data sets by specifications in the respective data definition file. 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 data and/or 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 dynamic interactions by providing reaction schemes. A model definition file is divided in several sections that are described in detail at the Setting up models page. Most prominently the differential equation system is defined in the REACTIONS section 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"
...
Each line represents a biochemical reaction. Reactants and products are specified on the left and right of the reaction arrow -> and the reaction rate is specified next. Here, the membrane associated receptor Rec either forms a complex IL13_Rec with the ligand IL13 (first line) or is internalized (second line). The internalized receptor Rec_i is again recruited to the membrane (third line). Each reaction rate is added to the right-hand-side of the ODEs for the products and subtracted for the reactants. The D2D software parses the specified reaction rates and considers each term that is not a dynamic variable nor an input as parameter, which can be estimated later 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 the observation functions are specified to make the default connection between dynamic variables and experimental data. In particular the observation functions are defined in the OBSERVABLES section 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 is finished once, the full functionality if the d2d framework can be applied for the analyses. All information about model and data are now available in the MATLAB environment as a global variable ar. Below some basic functions 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 user has several possibilities to specify which kind of plots should be shown by the default plotting command arPlot. After the setup, one figure for each data-set is shown (one in total for this example).
One way to also visualize 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 the 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.
All commands executed in the d2d framework via a graphical user interface can of course also be perfomed at the MATLAB command line which enables the creation of analysis workflows.
arFit starts the selected numerical optimization algorithm to calibrate the model to the experimental data. Maximum likelihood estimation is the default fitting strategy and is used here. If the the option flag in the global variable ar.config.showFitting = true model calibration is shown for each iteration by arPlot. In this example, fitted parameters values were already loaded from a workspace in the setup file. Therefore arFit does improve the likelihood further.
To check whether the current local optimum is also globally optimal a sequence of n fits for random initial guesses is executed using the command arFitLHS(n). For deterministic optimization the initial parameter guesses are generated by Latin hypercube sampling (LHS). arFitLHS replace previous values if a better fit is found.The resulting likelihoods of the n fits are displayed in a sorted manner by arPlotChi2s.
This figure shows the result for 100 different fits and indicates 1) the improvement of the initial likelihood value by several orders of magnitude (gray crosses vs. red circles), and 2) the plateaus after optimization demonstrating that local minima are reliably detected. The more often the same local minima are found, the better the chance that all existing local minima including the global one were found.
A convenient way to assess whether a parameter of a non-linear model has a finite confidence region (i.e. is identifiable) is by calculating the profile likelihood. This concept generalizes the concept of standard errors and Fisher-Information for the non-linear setting. After initializing the profile likelihood calculation by arPLEInit, the profile for all parameters are calculated automatically by ple. The calculation for a specific parameter, e.g. for the internalization rate Rec_intern might be performed by ple('Rec_intern',m). The second argument m can be used to specify the maximum number of steps (default: 50) along the profile. For this example, 200 points in lower an upper direction were calculated via ple('Rec_intern',200).
In the upper panel the profile likelihood is plotted. The parameter is identifiable because the profile (blue line) crosses the statistical threshold (red) at both sides of the optimum. The lower panel shows along the profile, i.e. for specific values of Rec_intern how all other parameters have to be adjusted to describe the data. This dependencies are of particular interest for non-identifiabilities as demonstrated in the following exemplary profile.
The second parameter examined for identifiability is the initial concentration of the IL13 receptor init_Rec_i. Here, the default settings only shows a part of the profile of interest. The following commands quickly provide the whole profile likelihood.
#!matlab
k = arPrint('init_Rec_i'); % Find index of init_Rec_i
global ar; global pleGlobals; % Get access to global variables
ar.ub(k) = 5; % Set upper boundary to 10^5 for optimization
arPLEInit; % Initialization, e.g. adopting ar.ub for the profile likelihood calculation
ple(k,[],[],.1,.1) % Calculate profile with fixed stepsize of 0.1
The resulting figure shows the profile likelihood of init_Rec_i.
The following conclusions can be drawn from this plot:
- The parameter has a lower boundary at around 2 (on the log10-scale)
- The profile shows only a minor increase for values > 4
- The profile increases at around 4.7 because another parameter can not be decreased further to compensate depicted by the circles. To check this statement, the profile is recalculated with larger boundaries for the receptor recycling
#!matlab
j = arPrint('Rec_recycle'); % Find index of Rec_recycle
ar.lb(j) = -8; % Set lower boundary to 10^-8 for optimization
pleGlobals.lb(j) = 5; % Set lower boundary to 10^-8 for profiling
ple(k,[],[],.1,.1) % Re-calculate profile of init_Rec_i
This demonstrates that init_Rec_i can not be restricted to high values, hence is practically non-identifiable.