How do I do parameters Estimation in SimBiology using Optimization toolbox?

Hi

I am trying to perform parameter optimization in SimBiology using optimization function from the optimization toolbox. As an test, I tried to do parameter estimation shown in the 'sbioparamestim' example here - http://www.mathworks.com/help/simbio/ref/sbioparamestim.html to check if I get the same answers. My code is as follows, the driver file (Estimation_Driver.m)

%%Driver to test sbioparamestim using optimization toolbox

sbioloadproject('gprotein_norules');

p_array = sbioselect(m1,'Type','parameter');

for index = 1:size(p_array,1)
    k0(index,1) = p_array(index).Value;
end

[error] = Objective_func(k0,m1);

f = @(k)Objective_func(k,m1);

options = optimset('Display','iter');
[k,res_error] = lsqnonlin(f,k0);

And the file calculating the error is as follows (Objective_func.m)

function [error] = Objective_func(k,m)

Gt = 10000;
tspan  = [0 10 30 60 110 210 300 450 600]';
Ga_frac = [0 0.35 0.4 0.36 0.39 0.33 0.24 0.17 0.2]';
xtarget = Ga_frac * Gt;

% Assign the parameters
for index = 1:length(k)
    m.Parameters(index).Value = k(index,1);
end

[soln_obj] = sbiosimulate(m);

[t1,Ga,~] = selectbyname(soln_obj,'Ga');

Ga_sim = interp1(t1,Ga,tspan);

error = xtarget - Ga_sim;

My question is - why are my answers (below) different than the answers in the example?

    ans =

      0.0100
     -0.0000
      0.0004
      4.0000
      0.0040
      1.0000
      0.0000
      0.1100

and the residual is

1.2176e+06

which are different from the values on the 'sbioparamestim' help page, and even the termination message are different :(my message - below):

Optimization running.
Objective function value: 1217576.5248185194
Local minimum possible.

lsqnonlin stopped because the size of the current step is less than
the default value of the step size tolerance.

and the message when I run using 'sbioparamestim':

Local minimum possible.

lsqnonlin stopped because the size of the current step is less than
the default value of the step size tolerance.

Stopping criteria details:

Optimization stopped because the norm of the current step, 5.676554e-11,
is less than options.TolX = 1.000000e-06.

Optimization Metric                                 Options
norm(step) =   5.68e-11                      TolX =   1e-06 (default)

Am I not using the same conditions as in 'sbioparamestim' (I didn't change any of the default conditions in 'lsqnonlin'?

More importantly, is this the correct way to use an optimization solver if I don't want to use 'sbioparamestim' (I find it easier to pose my problem this way, for simulating experiments which are much more complex)?

Any help here would be useful. Thanks!

SN

[Ingrid Tigges]

Which release of MATLAB are you using? I had a look into the example but it only contains the lines

sbioloadproject gprotein_norules m1;
Gt = 10000;
tspan  = [0 10 30 60 110 210 300 450 600]';
Ga_frac = [0 0.35 0.4 0.36 0.39 0.33 0.24 0.17 0.2]';
xtarget = Ga_frac * Gt;
p_array = sbioselect(m1,'Type','parameter');
Ga = sbioselect(m1,'Type','species','Name','Ga');
[k, result] = sbioparamestim(m1, tspan, xtarget, Ga, p_array)

And not all the other code you provided.

[S N]

I am using

MATLAB - R2012b

SimBiology - Version 4.2 (R2012b).

The code you wrote is the what you will find on the help page, whereas the code I wrote in my question is what I wrote to use do parameter estimation without using 'sbioparamestim'.

Answer 1

[Arthur Goldsipe]

Hi,

Calling lsqnonlin via sbioparamestim uses slightly different default options than calling lsqnonlin directly. The reference page and the function help list those differences. Most important are the different values for TolFun, FinDiffRelStep, and TypicalX.

As you can see, without these defaults the optimization can converge to a poor minimum. This can happen when using gradient-based optimization methods to estimate parameters of ordinary differential equations (ODEs) that are solved numerically. The accuracy of the ODE solution (as determined by the solver tolerances) limit the finite differencing step size that the optimizer can use while estimating the gradients.

If you wish to call optimization methods directly, I would recommend either tightening the solver tolerances and/or changing the option FinDiffRelStep. Your objective function is reasonable but there are a couple of things you can do to speed it up dramatically. First, you can get rid of the interpolation step and instead use the OutputTimes configset option (see the OutputTimes reference page for details). In addition, you can dramatically speed up the call to sbiosimulate by using sbioaccelerate, perhaps making the simulations run in 1/20th the time. (This may require you to install a compiler and run "mex -setup".) I would re-structure the diver to load the model, make any necessary configuration changes (like changes to OutputTimes or tolerances), call sbioaccelerate, and then call the optimizer. This is the sort of work that sbioparamestim keeps you from having to do.

-Arthur

[S N]

Thanks Arthur for the answer and the tips.

I am going to try with relaxed tolerances. Currently I am simulating ~40 different experiments simultaneously (different parameters and Initial Conditions in each one of them), most of which calculate the time at which AUC (Area Under the Curve) reaches a certain value, hence I am not sure 'OutputTimes' would be of that much help.

But, I am going to try 'sbioaccelerate' and see if it makes any difference.

Thanks S N

Answer 2

[Shashank Prasanna]

Optimization is tricky business. You will have to see why you solver stopped and then relax the tolerances till you find a better result. You may end up getting stuck in local minima depending on your initial conditions as well.

This message should give you the hint:

Local minimum possible.

lsqnonlin stopped because the size of the current step is less than
the default value of the step size tolerance.

Try relaxing the tolerances for TolX and try again.

I recommend going through this nicely written page which has good steps to follow when you hit a wall this way. http://www.mathworks.com/help/optim/ug/when-the-solver-fails.html

[Sietse Braakman]

Here is a general guidance on parameter estimation with SimBiology:

https://www.mathworks.com/matlabcentral/discussions/simbiology/134-guidance-and-troubleshooting-for-parameter-estimation-in-simbiology