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

Hi,

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

%Run SAEM

[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 SimBiology.fit.internal.getYPredictFcn/calcYPredict (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

[Arthur Goldsipe]

Hi,

My initial guess is that this sounds like a bug related to observables. It's difficult to tell when I can't reproduce the problem myself. Could you share the files necessary to reprodue the problem? (You can attach files to your question.) Can you also confirm what version of MATLAB you are using?

Thanks,

Arthur

[Budhaditya Chatterjee]

Thank you, Arthur for the reply. I have attached the script and the data. Also the version information is added. Thanks in advance.

Answer 1

[Arthur Goldsipe]

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.

-Arthur

[Budhaditya Chatterjee]

Thank you very much, Arthur.