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 of 100 initial parameter sets 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.
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. It generalizes the concept of inverse Fisher-Information to the non-linear setting. After initializing by arPLEInit the profile of the parameter k is calculated by ple(k,m). The second argument m specifies the maximum number of steps (default: 50). For this example, the first profile likelihood to look at is for the parameter Rec_intern calculated by the command ple(k6,200).
In the upper panel the profile likelihood is plotted. The parameter is identifiable because the profile (blue) crosses the statistical threshold (red) at both sides of the optimum. In the lower panel the variations of all other parameters along the profile are shown. This figure is especially informative for non-identifiabilites as demonstrated in the following exemplary profile.
The second parameter investigated for identifiability is the initial concentration of the IL13 receptor init_Rec_i. Here, the default settings are not optimal. The commands provided below lead quickly to a satisfying result.
#!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
pleGlobals.ub(k) = 5; % Set upper boundary to 10^5 for profiling
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
- The profile shows only a minor increase for x-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.