MATLAB Answers

How can I specify some bounds for parameter estimation in SimBiology?

7 views (last 30 days)
I am trying to implement parameter estimation using the function SBIOPARAMESTIM from Simbiology but I was not able to find information on how to specify some bounds for the parameter values I aim to estimate. I would like to know if it is possible to specify some bounds (lower & upper) for parameter estimation under Simbiology.
 
Specifying parameter bounds is required to the numerical identification of our models, which have a large number of parameters and could not be calibrated without some 'a-priori' information on their values.

Accepted Answer

MathWorks Support Team
MathWorks Support Team on 17 Mar 2016
Edited: MathWorks Support Team on 17 Mar 2016
For MATLAB R2014b and later, specifying parameter bounds during parameter estimation is directly supported by using "sbiofit" (from the command line) and the Fit Task (from the SimBiology Desktop):
 
1) Load your data
load(fullfile(matlabroot,'examples','simbio','data15.mat'))
 
2) Create a SimBiology model or use an existing model
pkmd = PKModelDesign;
pkc1 = addCompartment(pkmd,'Central');
pkc1.DosingType = 'Bolus';
pkc1.EliminationType = 'linear-clearance';
pkc1.HasResponseVariable = true;
model = construct(pkmd);
configset = getconfigset(model);
configset.CompileOptions.UnitConversion = true;
3) Define dosing
 
dose = sbiodose('dose');
dose.TargetName = 'Drug_Central';
dose.StartTime = 0;
dose.Amount = 10;
dose.AmountUnits = 'milligram';
dose.TimeUnits = 'hour';
4) Specify the parameters to estimate
paramsToEstimate = {'log(Central)','log(Cl_Central)'};
 
5) Specify the boundaries of the parameters
estimatedParams = estimatedInfo(paramsToEstimate,'InitialValue',[1 1],'Bounds',[1 5;0.5 2]);
 
6) Estimate the parameters
fitConst = sbiofit(model,gData,responseMap,estimatedParams,dose);
 
You can find the complete example in the attached Fitting_with_Bounds_Example.m.
For MATLAB R2014a and older the above approach above cannot be used. For these releases use the following approach:
The basic steps are as follows:
1) Load your .sbproj file
sbioloadproject('Tumor_Growth.sbproj', 'm1')
2) Load your data
ds = dataset('xlsfile', 'Data.xlsx') ;
3) Create vectors containing the input arguments for FMINCON. Among them you can specify upper and lower boundaries for the parameters as well as the initial values
p0 = [0.1, 0.1, 0.01] ; % initial estimate
lb = [1e-3 1e-3 1e-3] ; % lower bound
ub = [10 10 10] ; % upper bound
A = [] ; %
b = [] ; %
Aeq = [] ; %
beq = [] ; %
nlc = [] ; %
4) Define the optimization options
opt1 = optimoptions(@fmincon , ...
'Display' , 'iter' , ... %Iterations should be displayed
'Algorithm' , 'interior-point' , ... %Optimization Algorithm that should be used
'FinDiffRelStep' , m1Exp.SimulationOptions.RelativeTolerance ,.... %Step size factor
'TolFun' , 1e-3 ,... %Termination tolerance on the functionvalue, a positive scalar
'TolX' , 1e-6) ; %Termination tolerance on x,a positive scalar
5) Create your regression function
a) Set values of parameters to be estimated to value in current iteration, p
b) Set dose, amount and time
c) Call SIMULATE to simulate the model for this parameter set
d) Remove NaNs from the output
e) Calculate the sum of the squared errors
function f = objectiveFcn(p, yObs, m1Exp, ds)
% p - parameter
% yObs - Observed response
% yPred - predicted response
% f - value of the objective function
%%Simulation
m1Exp.InitialValues = p ; % Set values of parameters to be estimated to value in current iteration, p
d1 = getdose(m1Exp) ; % get dose
d1.Amount = ds.Dose(~isnan(ds.Dose)) ; % set amount
d1.Time = ds.Time(~isnan(ds.Dose)) ; % set time
m1Exp.SimulationOptions.OutputTimes = ds.Time ; % Set sampling time
[~, yPred] = simulate(m1Exp, d1) ; % Simulate
idx = ~isnan(yObs) ; % Identify non-nan values in yObs
%%Objective function: sum of squared errors
f = sum((yPred(idx) - yObs(idx)).^2) ;
end
6) Call your fitting algorithm
fh = @(p) objectiveFcn(p, ds.Tumor, m1Exp, ds) ;
pFit_FMC = fmincon(fh, p0, A, b, Aeq, beq, lb, ub, nlc, opt1 ) ;
You can also use another fitting algorithm that supports boundaries, for example lsqnonlin.

  0 Comments

Sign in to comment.

More Answers (0)

Tags

No tags entered yet.

Products


Release

R2016a

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!