reproducing sensitivities from MATLAB GUI

Hello,

I am trying to do some basic sensitivity analysis in MATLAB (command line) using the exported traces of a sensitivity analysis task from the Simbiology Desktop GUI (please see attached screenshot for summary of what's exported from GUI task to workspace).

Upon export, I have my sensitivities (dx(t)/dk(t)) in xs and my states (x(t)) in x. Now, if I run:

abs(trapz(t,xs(:,1,1)))

for instance, I can easily reproduce the non-normalized value for the sensitivity of species 1 to parameter 1 calculated from the GUI (not shown), and so on for every species/parameter combo. However, when I try to compute half-normalized sensitivities using:

abs(trapz(t,xs(:,1,1)./x(:,k)))

where k is the index for species 1 (for some reason these are not automatically aligned with xs and so k is not necessarily equal to 1), I'm unable to reproduce. I'm not sure if I'm missing something really simple here, and would appreciate the guidance. Thanks in advance!

Answer 1

[Sietse Braakman]

A typical use case for this question is when calculating a sensitivity on species C, where species C is a linear combination of species A and B (e.g. C = A+B) and is calculated using a repeated assignment. Currently (R2019a), SimBiology does not support calculating sensitivities on species that are the subject of a repeated assignment. Instead, you can use differentation rules to calculate sensitivity dC/dp from dA/dp and dB/dp (e.g. sum rule for dC/dp = dA/dp + dB/dp).

The order in which the sensitivities are calculated depends on the order in which you select the inputs and outputs for the sensitivity analysis. This holds for both sensitivity analyses set up in the task editor as well as those set up programmatically.

Automatic solution: This might or might not apply to your case. SimBiology allows for the sensitivities to be normalized, either just with respect to the numerator (i.e. the response or output of your model, this is called 'Half' normalization) or with respect to both the numerator and the denominator (i.e. the parameter or input to your model, this is called 'Full' dedimensionalization). See also this link. Set the normalization name-value pair on the configuration set to 'Half':

cs = getconfigset(model);
cs.SensitivityAnalysisOptions.Normalization = 'Half';

You can also do this in the task editor:

Manual solutions: In your case, it might just be a matter of rearranging the output order in your sensitivity analysis to match the order of the logged species during simulation (those in x).

Another option is to write a piece of code that matches the indices of the numerators in xsnames with those in names and perform the normalization using that. The function contains can help here.

indices = contains(xsnames,names{1});

will give you all the indices where names{1} is contained in the string of each entry in xsnames. For example

>> xsnames{1}
ans = 'd[Plasma.Drug]/d[K_plu]'
>> names{1}
ans = 'Plasma.Drug'
>> contains(xsnames{1},names{1})
ans =
  logical
   1

[Jeremy Huard]

To add up to the answer given by Sietse to compute the half-normalized time-dependent sensitivities, if you wish to then compute the time-independent sensitivities, you will need to compute the AUC of the absolute values of the time-dependent sentitivities.

So basically trapz(t,abs(xs(:,1,1))) instead of abs(trapz(t,xs(:,1,1))).