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

Question

Kendric Aaron 2023-8-28

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2013747-questions-regarding-sensitivity-analysis-for-simbiology-equation-value-range

评论： Jeremy Huard 2023-10-30

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

What is the specific mathematical equation used for computing the sensitivity coefficient of a parameter?
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!

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Jeremy Huard 2023-8-29

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2013747-questions-regarding-sensitivity-analysis-for-simbiology-equation-value-range#answer_1295856

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:

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.
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

12 个评论
显示 10更早的评论隐藏 10更早的评论

Jeremy Huard 2023-8-31

编辑：Jeremy Huard 2023-9-12

在 MATLAB Online 中打开

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.

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.

Jeremy Huard 2023-9-4

编辑：Jeremy Huard 2023-9-13

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 2023-10-24

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 2023-10-30

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

请先登录，再进行评论。

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

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

回答（1 个）

12 个评论
显示 10更早的评论隐藏 10更早的评论

社区

另请参阅

类别

标签

产品

版本

Community Treasure Hunt

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

0 个评论 显示 -2更早的评论隐藏 -2更早的评论

回答（1 个）

12 个评论 显示 10更早的评论隐藏 10更早的评论

社区

另请参阅

类别

标签

产品

版本

Community Treasure Hunt

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

12 个评论
显示 10更早的评论隐藏 10更早的评论