How to set up estimation method in sbiofit

Hi All,

I'm trying to estimate parametrs in my model using parameter estimation task in SimBiology toolbox.In the documentation available here it is mentioned that

"The default estimation method that sbiofit uses will change depending on which toolboxes are available."

From the list of methods available here, I want to select fmincon .

Could someone explain how to set up fmincon in the following syntax?

sbiofit(model,gData,responseMap,estimatedParam,'fmincon')

I want to define an objective function and pass the following arguments to fmincon

fmincon(@objective,p0,A,b,Aeq,beq,lb,ub,nlcon)

I have a confusion here,

p0 is the initial value of parameters, A,b,Aeq,beq are linear equality and inequality constraints, which is (A=[ ],b=[ ],Aeq=[ ],beq=[ ]) in my case.lb and ub are the bounds of paramaters and ncoln is the nonlinear constraints that are the differential equations taht will be generated from my model present in xml file.

I've already specified p0,lb,ub in estimatedParam array with properties Name, InitialValue and Bounds

estimatedParams     = estimatedInfo(paramsToEstimate,'InitialValue',ParameterInitialValues,'Bounds',ParameterBounds)

It's not to me how to define

fmincon(@objective,p0,A,b,Aeq,beq,lb,ub,nlcon)

since p0,lb,ub is already defined in estimatedParam array. However, I'd like to know

how to pass the other arguments objective ,A,b,Aeq,beq, nlcon in fmincon.

nlcon function is defined like this

Any suggestions will be highly appreciated.

[Jeremy Huard]

The estimatedInfo object will let you define the parameters initial values and their lower and upper bounds.

But can you explain what you would like to achieve with nlcon?

[Deepa Maheshvare]

Hi, Thanks you very much for the response. I want to pass ode's as constraints to the minimization problem as explained here. I'd like to implement it like this. The lb and ub in the inputs of the nonlincon(ncoln) function defined here are the bounds of the variables(Please note: lb and ub in nlcon are not the bounds of parameters that are already specified in estimatedParam) in my model. I'm trying to restrict the values of the solution of the variables at all times to be in a certain range[lb ub] . The range[lb ub] is the only experimental input that I have. [lb ub] is the range of the values of transient variables obtained from different experimental observations . i.e if there are 20 variables in my model, there are 20 differential equations generated form the sbml/xml file and the [lb ub] is available for all 20 variables.These bounds are steady state measurements from experiments.

In addition to this, I have the initial values of the parameters to be estimated and the [lb ub] of the parameters.

I don't have time varying experimental input to assign to gData,responseMap in sbiofit.

I'm not sure how to implemet this in SimBiology Toolbox. If there is a way to generate the ode's from the modelObject generated after reading the sml file using sbmlimport , I would be able to implement the above easily.

Any suggestions will be highly appreciated

Answer 1

[Jeremy Huard]

Hi Deepa,

let me summarize my understanding of your problem and please correct me if I'm wrong:

• you only have steady state data for several experimental conditions
• there is no other data to fit the model to
• you also have parameter bounds

Is this correct?

If yes, I suggest to create a dataset with these steady state data as explained here. If this is experimental data, you will know when this data was measured so you can build a table with time points and measurements.

Then you can use sbiofit to fit your model to those data points. Again, the App will be of great help.

[Deepa Maheshvare]

Hi Jeremy, Yes, you are right. I had a look at the dataset table . I'm afraid I don't have the time (independent ) data at which the steady state values have been observed. There are around 500 values and it would be impossible for me to get the steady state times. That's why I prefer to use the ncoln function in fmincon. Also,I am not sure now to define the objective function and define the odes as constraint.

[Jeremy Huard]

So, if this steady state data is the only data you have and you can't use it as such, I wonder how your objective function would look like. You need to give fmincon an objective function to minimize, so what would it be?

Also, you observe a steady state after some time but you don't really know what happens before. The model concentrations might have a complicated behaviour (maybe oschillations) that get dampened over time and reach a steady state after a while, which is what you observe. So your constraints might too strong.

I see two ways to tackle your problem:

1. either you create a dataset and define a (reasonable) large time as being your the steady state time. In that case you can use sbiofit which will be more convenient and will give you the possibility to call sbioparameterci to get confidence intervals of parameters for example
2. or you use fmincon or lsqnonlin/lsqcurvefit without nonlinear constraints but with an objective function that calls the function sbiosteadystate and compute the differences with your observed data.

[Deepa Maheshvare]

This is how my objective function is defined.It is the squared difference between the values computed from the model at a sufficiently long time tSpan, which would allow the system to reach steady state,(this I consider to be the steady state value from the model) and the experimental observation summed over all species. This is essentially the same approach that has been suggested by you in the second point where sbiosteadystate function can be used. I will go ahead with the second suggestion. Thanks a lot. If there is any way to implement non linear constraint, that would be my first choice.

[Jeremy Huard]

What you describe is actually solution 1. You could create a data set with a single time point at tspan. I believe this is the easiest way to get access to loglikelihood, confidence intervals, get the results in a variant for instance. But both approaches are valid.

[Deepa Maheshvare]

Hi Jeremy,

Thanks a lot for the response. I have a few doubts in the first approach.

When I do this ,"You could create a data set with a single time point at tspan" . When I do this ,

would it be incorrect, in the sense, not all species will take the same amount of time to reach steady state. For instance, I have computed transient times of species of all species in the model. The time scale ranges from 0.01 seconds to 1000seconds. When the time scale of equilibration is different for each species, I am not sure how "a data set with a single time point at tspan" can be justified.

I'd be happy to know your views on this.

[Jeremy Huard]

This is what you are currently doing by calling ode15s for [0,tspan] and taking the last data points.

And this is fine in my opinion. All you know from your data is that it was taken at steady state but you don't know exactly when. So all you need to do is take tspan large for everything in your system to be at steady state.

The idea with the data set with a single time point is to reproduce what you actually have, which is a single time point.

[Deepa Maheshvare]

Thanks a lot for the advice. I'll try it out