Questions Regarding Sensitivity Analysis for Simbiology (Equation, Value Range)

I have some questions regarding how sensitivity analysis is performed in SimBiology.

1. What is the specific mathematical equation used for computing the sensitivity coefficient of a parameter?
2. Normally what value range for the parameter (e.g. ± 5% of original value) is used for the computation? Is there any way to tweak this?

I cannot seem to find any reference for #2 specifically.

Thank you!

Answer 1

[Jeremy Huard]

Hi @Kendric Aaron,

I will assume that your question relates to Local Sensitivity Analysis. But please let me know if this is not the case.

Both questions make me think of an approximation of the sensitivities using finite differences. But let me start off by saying that SimBiology does not use the finite differences approach because it does not yield very accurate results. Instead, it uses a technique called complex step differentiation.

Let me now answer your specific questions:

1. The equation of the sensitivities being approximated depends on which normalization has been selected. Typically, it is recommended to use Full Dedimensionalization to make sure all sensitivities are dimensionless and be able to compare them. In this case, it will be . That being said, this equation is approximated by the complex step differentiation technique and the final equation used will differ. More details in #2.
2. You can find a nice short description of the algorithm and the final equation in Cleve Moler's blog, whose pioneer work in the 60's laid the groundwork for today's algorithm. Complex step differentiation also requires a step size parameter, which is set to a very small value because complex step differentiation is not only more accurate than finite differences but also more robust with decreasing step size.

Please note that complex step differentiation requires the ODEs to be complex analytic, that is, to be infinitely differentiable in the complex plane.

So, models with

• Nonconstant compartments
• Algebraic rules
• Events
• or nonanalytic functions

are not supported.

In this case, you can run a parameter scan to compute sensitivities with finite differences. For example, by varying parameter values by 5% as you suggested.

I hope this helps.

Best regards,

Jérémy

[Kendric Aaron]

Oh I see, I think the parameter scan might be what I am looking for. In that case, I would like to know how I can calculate the area underneath the curves generated by the parameter scan so that I can derive the sensitivity coefficients.

Thank you!

[Jeremy Huard]

Hi @Kendric Aaron,

you can ask SimBiology to compute an AUC after each simulation by adding an Observable to your model. Observables are computed post-simulation and will use the full timecourse to compute a scalar or time-dependent quantity. Plese note that you can also add an observable to simulation results instead of the model itself.

That said, in your case you need to not only to simulate the model for each parameter variation but also compare to the baseline result without any parameter perturbation.

This will be easier to achieve with a script. Here is a generic example that you should be able to adapt to your model:

% load model
s = sbioloadproject("tmdd_with_TO.sbproj");
model = s.m1;

% add observable to compute AUC of Target Occupancy
addobservable(model, "AUC", "trapz(time, TO)", Units="hour");

% simulation setup
parameterNames = ["ksyn", "kdeg", "kel", "kon", "km"];
output         = "AUC";
doseName       = "Daily Dose";

dose = getdose(model, doseName);
simfun = createSimFunction(model, parameterNames, output, dose, AutoAccelerate=false)

Warning: Dosing information includes numeric values for time, amount, or rate that are ignored.

Warning: Reported from Dimensional Analysis:
Cannot perform dimensional analysis for observable 'AUC' because of the function 'trapz' in the observable. Because UnitConversion is on, correct simulation results will depend on this expression being dimensionally correct. Additionally, SimBiology simulates the model in a unit system determined at runtime. The units are determined by the units used in the model and the model's configset. Unless the inputs and outputs to the function are dimensionless, results may change due to configset option changes, changes to the model, or version changes in SimBiology. It is recommended that input and output arguments to functions be dimensionless to ensure correct results.

Warning: Reported from Dimensional Analysis:
Cannot perform dimensional analysis for observable 'AUC' because of the function 'trapz' in the observable. Because UnitConversion is on, correct simulation results will depend on this expression being dimensionally correct. Additionally, SimBiology simulates the model in a unit system determined at runtime. The units are determined by the units used in the model and the model's configset. Unless the inputs and outputs to the function are dimensionless, results may change due to configset option changes, changes to the model, or version changes in SimBiology. It is recommended that input and output arguments to functions be dimensionless to ensure correct results.

simfun = 
SimFunction

Parameters:

      Name      Value         Type                  Units         
    ________    ______    _____________    _______________________

    {'ksyn'}         1    {'parameter'}    {'nanomole/hour'      }
    {'kdeg'}    0.0934    {'parameter'}    {'1/hour'             }
    {'kel' }     0.523    {'parameter'}    {'1/hour'             }
    {'kon' }    0.0485    {'parameter'}    {'liter/nanomole/hour'}
    {'km'  }    0.0458    {'parameter'}    {'1/hour'             }

Observables: 

     Name           Type          Units  
    _______    ______________    ________

    {'AUC'}    {'observable'}    {'hour'}

Dosed: 

      TargetName                 TargetDimension          
    _______________    ___________________________________

    {'Plasma.Drug'}    {'Amount (e.g., mole or molecule)'}


TimeUnits: day

accelerate(simfun);

% generate parameter values for simulation by adding 5% to each parameter
baseValues = arrayfun(@(x) sbioselect(model, Name=x).Value, parameterNames);

nPars = numel(parameterNames);
percentChange = 5/100;  % percentChange=5%
multiplier = [ones(1,nPars); ones(nPars) + percentChange*eye(nPars)];
samples = baseValues .* multiplier; % one row per simulation, one column per parameter
samplesT = array2table(samples, VariableNames=parameterNames, ...
    RowNames=["baseline","baseline + Δp" + (1:nPars)])

samplesT = 6×5 table
                      ksyn     kdeg        kel        kon         km   
                      ____    _______    _______    ________    _______

    baseline             1     0.0934      0.523      0.0485     0.0458
    baseline + Δp1    1.05     0.0934      0.523      0.0485     0.0458
    baseline + Δp2       1    0.09807      0.523      0.0485     0.0458
    baseline + Δp3       1     0.0934    0.54915      0.0485     0.0458
    baseline + Δp4       1     0.0934      0.523    0.050925     0.0458
    baseline + Δp5       1     0.0934      0.523      0.0485    0.04809

% simulate
stopTime = 8; % in simfun.TimeUnits (here, =day)
doseTable = getTable(dose); 
if iscell(doseTable)
    % transpose to get 1xN cell array in case of multiple doses to make
    % sure all are applied for each simulation
    doseTable = doseTable(:)';
end
results = simfun(samplesT.Variables, stopTime, doseTable);

% compute normalized sensitivities
resultBase = results(1).ScalarObservables.(output);  % first simulation
resultsPerturb = vertcat(results(2:end).ScalarObservables).(output); % all other simulations
sensitivities = abs(resultsPerturb - resultBase) ./ (percentChange*baseValues');
normalizedSensitivities = sensitivities .* baseValues' / resultBase;

% reorder in descending order
[normalizedSensitivities, idx] = sort(normalizedSensitivities, "descend");
parameterNames = parameterNames(idx);

% create bar plot
figure;
pNames = categorical(parameterNames);
pNames = reordercats(pNames, parameterNames);
bh = bar(pNames, normalizedSensitivities);

xtips = bh.XEndPoints;
ytips = bh.YEndPoints;
labels = string(bh.YData);
text(xtips, ytips, labels, HorizontalAlignment="center",...
    VerticalAlignment="bottom", Color=bh.FaceColor);

bh.Parent.YLim(1) = -bh.Parent.YLim(2)/20;
bh.Parent.YGrid = "on";
bh.EdgeColor = "none";
bh.BaseLine.Visible = "off";
ylabel("normalized sensitivity");
title("$$\frac{p}{AUC(p)}\frac{|AUC(p+\Delta p) - AUC(p)|}{\Delta p}$$", ...
    Interpreter="latex");

I hope this helps.

Best regards,

Jérémy

EDIT: Corrected sensitivities calculation with normalization.

[Kendric Aaron]

Just to be clear, the script needs to be executed from the main MatLab console? And I just need to copy-paste the script while changing the relevant information?

Also, for the getdose command, my issue is that my model has multiple active doses, so in my case, I guess I have to type all the active doses?

Thank you!

[Jeremy Huard]

Yes, you can create a new script or Live Script in the MATLAB Editor, copy/paste that code and run it.

The model used in this example is shipped with SimBiology, so the script should run without issues.

Concerning Active doses:

I would like to start by saying that it is not recommended to set the Active property to true. This property will be removed in a future release. It is recommended to specify doses explicitely when running simulations instead.

Now, if you already have multiple doses in your model with Active=true, you can get them with sbioselect:

doses = sbioselect(model, 'Active', true,'where','Type','==',["repeatDose","scheduleDose"])

Or you can get several doses by specifying their names:

doseNames = ["Daily Dose","Single Dose"];

doses = sbioselect(model, 'Name', doseNames)

Please note that sbioselect will not necessarily return the dose objects in the same order as the list of names. But this is not important for this example.

Best regards,

Jérémy

[Kendric Aaron]

For the loadproject command, I just have some question regarding the file name.

Currently, my file name is something like "Test Sensitivity.sbproj"

When I try to run the script (after making the necessary modifications), I get an error message stating that the program is unable to open the project. Do I need to change the file name to use underscores instead of spaces?

Also, it states something regarding file path error, even though the current active directory should be correct. I do not know what the problem might be.

Best wishes,

Kendric

[Jeremy Huard]

Instead of the command form I used in the example you can use the function form, which will support spaces in the name:

sbioloadproject("Test Sensitivity.sbproj")

[Kendric Aaron]

Hello,

Thank you for the information.

I have a question regarding the dose.TargetName.

Currently, the target name in my model is Liver.LiverBFR. However, this seems to cause an error when I attempt to run the script as it is unable to resolve the target name. How do I address this issue?

Thank you!

[Jeremy Huard]

Hi,

it is difficult to find the issue without the actual model.

My first advice would be to open the project in SimBiology Model Builder and check the dose to make sure that there is no red indicator next to the Target Name. A red indicator would mean that the target name is indeed not resolved and you could pick the right species instead.

Best,

Jérémy

[Kendric Aaron]

I have some question regarding the following part of the script:

simfun = createSimFunction(m1, parameterNames, output, dose.TargetName, AutoAccelerate=false)

accelerate(simfun);

As my model has multiple active doses,should I type "doses.TargetName" instead of "dose.TargetName"?

When I use "doses.TargetName", I receive an error:

Unrecognized method, property, or field 'TargetName' for class 'SimBiology.RepeatDose'.

What could be the reasons for this error, and how do I resolve it?

Thanks!

[Jeremy Huard]

As stated in the doc page of createSimFunction, if you have multiple doses to pass you should specifiy a list of names as cell array or string array or pass the list of doses directly.

So you could use:

simfun = createSimFunction(model, parameterNames, output, doses, AutoAccelerate=false)

or

targetNames = get(doses, "TargetName");

simfun = createSimFunction(model, parameterNames, output, targetNames, AutoAccelerate=false)

When it comes to the simulations, you will need to make sure that all doses are applied to each simulation. This can be achieved with:

doseTables = getTable(doses)'; % please keep the ' (=transpose) to make sure you

% get a cell array of size 1xN

results = simfun(samplesT.Variables, stopTime, doseTables);

Best,

Jérémy

[Kendric Aaron]

I was able to get the code to run properly and made some modification to suit my specific needs.

I would just like to ask, is there a preferred way to cite the forum post for official purposes (e.g. publication)?

Best wishes,

Kendric

[Jeremy Huard]

Hi @Kendric Aaron,

thank you for the update. I am glad to hear that it worked for you.

I am not aware of a specific way to cite forum posts. But I would assume that citing it like any other URL would suffice. Bibtex has a dedicated URL entry in case you use Latex.

As title you could use 'SimBiology Community Website'.

Best regards,

Jérémy