Global fitting using Kinetic scripts
From QuB
A practical tutorial by Lorin S. Milescu
Contents |
Overview
QuB provides several sophisticated optimizers for analyzing single channel data (MIL, MPL, MIP), and one for macroscopic currents (Mac). However, these algorithms are limited to only one type of data, i.e., single or macroscopic. It may be impossible to efficiently extract all the necessary information from either type. Moreover, the hard constraints that can be imposed on the kinetic parameters through the use of the Singular Value Decomposition, as provided by QuB, may not be enough. For example, this constraining mechanism can be used to enforce microscopic reversibility, or equal and identical binding sites, but is not suitable for enforcing soft constraints, such as a maximum pOpen, or a Kd value. Note that the hard constraints reduce the number of free parameters, whereas the soft constraints narrow the parameter space. Additionally, the dedicated QuB optimizers are limited to optimizing only the k0 and the k1 factors of rate constant expressions, and the so-called extra kinetic parameters, which can be used to implement allosteric factors [add page]. In practice, one may need to optimize in terms of derived quantities, such as half-activation voltage, etc.
QuB provides another, very flexible mechanism for general optimization problems. The idea is simple: let the user define the model, the parameters, the data to be fitted, and the cost function using Kinetic scripts, and then run a generic optimizer, such as simplex. Using this approach, one can optimize kinetic models from an arbitrary combination of single and macroscopic data, use a variety of free parameters, and enforce additional constraints. While very versatile, this general optimization procedure involves writing new scripts (and scriptlets), or modifying existing ones, which could be quite challenging.
In the following, I will discuss the necessary steps for setting up a general optimization problem in QuB.
Preparing the interface layout
Having the right interface layout is helpful. Depending on the particular problem, the following windows should be made visible in the interface: the Data, the Model, the Graphic, the Script, the Scriptlets, the Global Fitting and the Report. The Data and the Model windows are strictly necessary only if one of the single channel or macroscopic optimizers is called to calculate the cost function (or a component of it), in which case the data and the model files have to be opened. If too many windows are required, they can be split between two or more interfaces. For example, one interface may contain the Scripts, the Scriptlets, the Model and the Report windows, while another interface can display the Graphic. Even better, if two monitors are connected to the computer, the Graphic can be displayed as a floating window in the second monitor.
Preparing the data
Data used in the fitting can be loaded as QuB data files, or as traces loaded in one or several of the available pages in the Graphic window. Loading a trace can be done interactively, by using the Open trace button, or programmatically, by running a scriptlet. The Kinetic commands for loading a trace are the following:
//Declare an object called "graph", of type "graphic"
graphic graph;
//Load a single data trace from the "mytrace.txt" file, on trace 0 and page 1
graph.loadtrace(name = "mytrace.txt", trace = 0, page = 1);
//or load an entire list of traces; in this case, the trace parameter is disregarded
graph.loadtrace(name = "mytraces.lst", trace = 0, page = 1);
One way or another, just make sure that all the data that are to be used in the fitting are loaded in QuB.
Preparing the model
Conventionally, QuB allows only Markov models with Eyring rates. For example, the expression of a voltage-dependent rate constant is k = k0 x exp(k1 x V), where k0 and k0 are the rate at zero transmembrane potential (s-1) and the voltage dependency factor (mV-1), respectively. With the Kinetic scripting language, the user has full freedom to define and use other rate constant expressions, such as k = k0 x exp(k1 x (V - Vh)).
Regardless of the rate constant expressions used, having a visual model opened in QuB is not necessary, but it may be helpful.
Preparing the optimizer
The optimizer is interfaced in the Global Fitting window. The user can choose between simplex and dfpmin. Note that the dfpmin optimizer requires the gradients of the cost function, whereas simplex doesn't. Since analytical expressions for the gradients may be unavailable or too difficult to calculate, dfpmin uses numerical approximations. However, for some problems, even the numerical derivatives cannot be calculated accurately, or the cost function is simply not differentiable, in which case simplex is the only choice.
The free parameters must be specified in the main Kinetic script, using a specialized optimizer object, as illustrated in the following example:
//Declare optimizer object
optimizer MyOpt;
//Declare parameters to be optimized
par dbl Par1 = 1.0, Par2 = 0.1;
//Add parameters to the optimizer object
MyOpt.addpar(par = Par1, min = -10.0, max = 10);
MyOpt.addpar(par = Par2, min = -1.0, max = 1.0);
In the Global Fitting window, enter "MyOpt" in the "Params name" edit field. This way, the global optimizer knows
what to optimize.
In addition to the name of the optimizer object, the user must specify the name of several scripts (actually, scriptlets): "On start", "On end", "Iteration", "Penalty" and "Set par". Each one of these scripts is called by the optimizer at different times during the fitting,and their usage is up to the user. For example, I use "On start" to initialize the free parameters. This is particularly important if certain transformations between the parameters of the model and the free parameters are necessary. The best example is the k0 parameter (the rate constant at zero potential). Since k0 must be >= 0, its logarithm is optimized instead. Thus, the "On start" script contains the following code:
lnk0 = ln(k0);
Of course, lnk0 is a parameter that must be declared in the main script, along with k0:
par dbl k0, lnk0;
The "Set par" script can be used to do the inverse transformation, i.e., to convert the free parameters back into model parameters:
k0 = exp(k0);
The "Iteration" script can be used to display the progress of the cost function and of the optimized parameters during the fit.
The "Penalty" script can be used to calculate the cost function. The mechanism by which the value of the cost function is passed back to the Global optimizer is described in the next section.
Preparing the main script and the scriptlets
For the Global optimizer to work, one must have a Kinetic script loaded in the Scripts window, and this script must be running. In fact, unless there is a specific reason to have it running - such as real-time fitting of action potentials - the script can and should be paused. The reason is that the main script must contain, at the minimum, declarations of the optimizer object and of the model and free parameters, and possibly some other code. Clearly, if the main script is not running, these objects and variables are "dead". The main script might look like this:
Main script:
//Declare optimizer object
optimizer optNa;
//Declare parameters
par dbl am_k0, am_k1, bm_k0, bm_k1, ah_k0, ah_k1, bh_k0, bh_k1, am_lnk0, bm_lnk0, ah_lnk0, bh_lnk0;
//Declare ion channel object
ionchannel hh icNa;
//Set up the ion channel model
//. . .
//Add parameters to the optimizer object
optNa.addpar(par = am_lnk0, min = ln(1e-6), max = ln(1e6));
optNa.addpar(par = am_k1, min = -1.0, max = 1.0);
optNa.addpar(par = bm_lnk0, min = ln(1e-6), max = ln(1e6));
optNa.addpar(par = bm_k1, min = -1.0, max = 1.0);
optNa.addpar(par = ah_lnk0, min = ln(1e-6), max = ln(1e6));
optNa.addpar(par = ah_k1, min = -1.0, max = 1.0);
optNa.addpar(par = bh_lnk0, min = ln(1e-6), max = ln(1e6));
optNa.addpar(par = bh_k1, min = -1.0, max = 1.0);
//Use a "stop" statement to pause the script
stop;
beep;
A step-by-step example
To run this example, you will have to download several files from: http://lambda.med.buffalo.edu/milescu
Contact
For questions send me an email at: Lorin_Milescu at hms dot harvard dot edu
