How to fit an observable in a SimBiology model using sbiofitmixed function?

3 views (last 30 days)
I am trying to fit an observable of a Simbiology model with some available data. The observable is added to the model as the following:
o1=addobservable(IHY,'O1','log10(I./max(I))'); %Add observable
Where, I is a species in the model. It gives perfect results when I simulate it. However, when I am trying to fit the model to the data as the following:
index_Data = groupedData(index_data); %Create dataset
responseMap = 'O1 = index_R'; %Map the relation
estimatedParams = estimatedInfo({'a', 'b', 'c','d','e','f','g','I0','H0'}); %Declare which parameters to fit
[nlmeResults, simI, simP] = sbiofitmixed(IHY, index_Data, responseMap, ...
estimatedParams,[], 'nlmefitsa','UseParallel',true,'ProgressPlot',true);
It generates the following error:
Index in position 2 exceeds array bounds (must not exceed 3).
Error in (line 84)
yDataCell{i} = yDataCell{i}(:, responseIndices);
Error in nlmefitsa>calcf (line 1304)
f = fun(phiarg,X,V(grp,:));
Error in nlmefitsa>modelcaller (line 1206)
f = calcf(vect,fun,grp,phi,X,V);
Error in nlmefitsa>@(p,x,v,e)modelcaller(fvc,vect,modelfun,IdM,transphi(p),x,v,e) (line 596)
structural_model = @(p,x,v,e) modelcaller(fvc,vect,modelfun,IdM,transphi(p),x,v,e);
Error in nlmefitsa>starting_set (line 2130)
f = structural_model(phiM,XM,VM,errmod);
Error in nlmefitsa/onefit (line 601)
phiM = starting_set(nchains,chol_Gamma,mean_phi,XM,VM,IdM,...
Error in nlmefitsa (line 464)
[stop,betas(:,jrep),Gj,statsj,b_hat(idxRandPerm,uId,jrep)] = onefit(beta0(:,jbeta));
Error in sbiofitmixed (line 247)
[results.beta, results.psi, results.stats, results.b] = feval(estimationFunctionName, timeVectorStacked,
yVectorStacked, groupVectorStacked, groupList, yPredictFcn, beta0, nameValuePairs{:});
Error in observable_error_script (line 58)
[nlmeResults, simI, simP] = sbiofitmixed(IHY, index_Data, responseMap, ...
What am I doing wrong here? My data is such that I need to fit the observable, not a species. Note that my model does not have any dose. The dynamics is described by initial assignment and rate rules. The version information and matlab script to reproduce the problem is attached. Any help will be greatly appreciated.
Budhaditya  Chatterjee
Budhaditya Chatterjee on 28 Feb 2021
Thank you, Arthur for the reply. I have attached the script and the data. Also the version information is added. Thanks in advance.

Sign in to comment.

Accepted Answer

Arthur Goldsipe
Arthur Goldsipe on 1 Mar 2021
Edited: Arthur Goldsipe on 1 Mar 2021
HI Budhaditya Chatterjee,
This looks like a bug that someone else reported to me a few days ago. The bug is triggered whenever there is a simulation error, for example due to integration problems. The intended behavior is that the simulation errors are caught and NaN (not-a-number) values are reported for these simulation results.
Since I just found out about the bug, I can't yet tell you when it will be fixed. However, even if the bug was fixed, you would run into a different error due when running your code. So you should be able to avoid this bug once you've set up your fitting problem in a more appropriate way. I see two main issues with the problem as formulated.
First, your system of differential equations is quite sensitive to the parameters. I often encountered simulation errors during the fit, even for modest parameter changes. There are a couple of things you can do to avoid this problem. First, the parameters a through g can be constrained to have positive values, which you can do by estimating the log of these parameters (by setting the Transform property to 'log' in the appropriate elements of your estimatedInfo vector. This helps prevent ODE simulation errors. Second, I suggest updating the configuration set option AbsoluteToleranceScaling to false. The default value of true does not work well when fitting this sort of model.
Second, the observable you added to the model can evaluate to NaN, Inf, or complex values, if I becomes negative or is always 0. There are a number of ways you could address that. The simplest I could think of is to fit to I/max(I) by updating your model and data accordingly.
I tried these fixes, and I'm able to run the fit for quite some time without error. It looked like the fit would take some time to complete, so I didn't confirm it would run to completion. However, in case you find my changes helpful, I'm attaching them. You'll also notice that I made some minor changes to adapt the code to the way I would write it.
UPDATE: I left the fit running for a while, and it did eventually error with a simulation problem. So now I unfortunately think it may not be feasible to perform this fit.

More Answers (0)


More Answers in the  SimBiology Community





Community Treasure Hunt

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

Start Hunting!