Calculating Senstivities with SimBiology GUI and from 'sbiosimuate' - How do I deal with NaNs and Infs

Hi

I am trying to replicate the results of Sensitivity Analysis I get from Simbiology GUI in my script file. The code is written below:

% set the outputs
set(csObj.SensitivityAnalysisOptions, 'Outputs',[s1 s2 s3]);

% set the inputs
set(csObj.SensitivityAnalysisOptions,'Inputs', [k1 k2 k3]);

% Set Normalization
set(csObj.SensitivityAnalysisOptions,'Normalization','Full');      

% Enable Sensitivity
set(csObj.SolverOptions, 'SensitivityAnalysis',true);

% accelerate
sbioaccelerate(m);

% Simulate
simDataObj = sbiosimulate(m);

% Get the data
[T, R, snames, ifacs] = getsensmatrix(simDataObj);

The Sensitivity matrix I get ('R') has NaNs and Infs in it, so I convert them to zero and sum across the time axis.

R(isnan(R)) = 0;
R(isinf(R)) = 0;
R = abs(R);

R1 = sum(R,1);          % Sum across time
R1 = squeeze(R1);       % size of R1 now is states * parameters

And compare the results with the plot obtained from GUI

[Value, Index] = sort(R1(1,:),'descend')

I see that my results are different from those in the GUI, both in terms of values of Sensitivity and order of sensitive parameter.

Can anyone help as to what I am doing wrong here in the normalization of 'R' matrix?

Thanks S N

Answer 1

[Ingrid Tigges]

The results between the GUI based calculation and your calculation in the command line differ because of a difference in the calculation. Your command line version calculates the sum of sensitivity across time while in the GUI version the area under the curve is used. Thus if you replace your sum with the area under the curve the results will be the same.

A little more background information: The Solver computes the sensitivity of each output with respect to each input at all the time points that it steps through. So In the raw form sensitivity is a 3 dimensional matrix(time x n_Output x n_Input). This is the Matrix that getsensmatrix returns. If you want to create a time plot (time against sensitivity) you can use it as is. But if you want to create a bar chart to answer the question, to which input factor is an output most sensitive to, then you need to sum up the trajectory in some fashion to come up with one number for each output/input pair. We do this by numerically integrating the sensitivity trajectory in time:

[t, R, states, inputs] = getsensmatrix(tobj, x, y); [~, in, out] = size®; result = zeros(in, out); for i = 1:in for j = 1:out index = ~isnan(R(:,i,j)) & ~isinf(R(:,i,j)); result(i,j) = trapz(t(index), abs(R(index,i,j))); end end

[S N]

Thanks Ingrid, that solves my problem.

S N