Formatting for 'addobservables'

I've built a SIR epidemic model on SimBiology and want to calculate its sobol indices.

I have 2 parameters, alpha and gamma, which describe infection rate and recovery rate, and 3 species, 'S', 'I' and 'R' to represent the population compartments.

I want to determine peak infections after the simulation (what is the highest value for the 'I' species) and set that as the observable which I will the input into the sbiosobol function to determine effect of alpha and gamma on the value of peak infections.

This is the code:

m1=sbiomodel('SIR');

N=10000;

I0=10;

S=addspecies(m1,'S','InitialAmount',N-I0);

I=addspecies(m1,'I','InitialAmount',I0);

R=addspecies(m1,'R','InitialAmount',0);

r1=addrule(m1,'S = -beta*S*I/N','RuleType','rate');

r2=addrule(m1,'I = beta*S*I/N-gamma*I','RuleType','rate');

r3=addrule(m1,'R = gamma*I','RuleType','rate');

p1=addparameter(m1,'beta','Value',0.2);

p2=addparameter(m1,'gamma','Value',0.1);

p3=addparameter(m1,'N','Value',N);

T=200;

conf=getconfigset(m1,'active');

set(conf,'Stoptime',T,'SolverType','ode45');

set(conf.SolverOptions,'OutputTimes',[1:T]);

[t,pop,names]=sbiosimulate(m1);

imax=addobservable(m1,'imax','max(pop(:,2))');

sobol=sbiosobol(m1,{'beta','gamma'},'imax');

However, I am getting an error for the imax=... line which reads: "Name 'pop' in observable 'imax' does not uniquely refer to any species, parameters, compartments, or observables according to SimBiology precedence rules."

Thanks!!

Answer 1

[Florian Augustin]

Hi Shaun,

expressions of SimBiology observables can only reference names of components on the model, or functions on the Matlab path; referencing variables in the Matlab Workspace is not supported.

If I understand you use case correctly, you can reference "I" instead of "pop(:,2)" in the observable's expression:

addobservable(m1, "imax", "max(I)"); % <- reference I here instead of pops(:,2).

% Run Sobol analysis:

sobolResults = sbiosobol(m1, ["beta", "gamma"], "imax", "ShowWaitbar", true);

% Plot results:

bar(sobolResults);

I hope this helps.

-Florian

[Shaun Tan]

oh yes it works, thanks!

[Shaun Tan]

Hi Florian,

Do all syntax work in SimBiology?

I'm trying to perform analysis on another observable but I keep getting this error: "Unable to perform global sensitivity analysis. All samples contain incomplete or failed model simulations. Adjust simulation settings or parameter bounds to values that result in successful simulations."

Now I'm trying to find the day on which this infection peak occurs, and the value can be calculated successfully using:

imax1=max(pop(:,2));

iday1=find(pop(:,2)==imax1);

However when I express it in the addobservable format shown below it throws up the error, I think it's got to do with the find() expression?

iday=addobservable(m1,'iday','find((I)==max(I))');

Thanks,

Shaun

[Florian Augustin]

Hi Shaun,

You are very much on the right track with the observable.

If I understand the use-case correctly, the observable should return the time at which the max. I is attained. You can do this by indexing into the time variable, which is a reserved token in SimBiology that gives you access to simulation times of the ODE solver.

When evaluating expressions of observables, SimBiology first tries to determine whether the result is a scalar value or a vector. I therefore recommend to write observables' expressions in a way that they either always return a scalar value or always a vector of values, for any value of I. In your case, find(I == max(I)), returns a scalar if the max. is unique, but it can also evaluate to a vector if the max. is not unique.

If you know that the max. is unique in all your model simulations, you can force the index to be unique by choosing the first occurance as follows:

addobservable(m1, "iday", "time(find(I == max(I), 1))");

If you are unsure and your model could attain the max. value at several different time points and you want to be alerted if this happens, you can write a helper function (an example, getTimeOfMaxValue.m, is attached) to guard against this case:

addobservable(m1, "iday", "getTimeOfMaxValue(time, I)");

Best,

Florian

[Shaun Tan]

Hi Florian,

Thank you very much for your explanations and suggested improvements, the sensitivity analysis is up and running now. One last question concerning the concept of global sensitivity analysis:

If I change the input parameters, shouldn't the sobol indices of the output remain the same since it's global? If the indices change wouldn't that be local sensitivity analysis instead? Right now running the analysis for a particular output while changing the input parameters leads to different sobol indices.

Thanks again!

Shaun

[Florian Augustin]

Hi Shaun,

You are right, changing the (distributions of the) input parameters will generally change the Sobol indices.

Sobol indices are considered a global analysis because they reflect the sensitivity of a model's response over a whole range of input parameters (as specified by the input parameters' distributions). This is in contrast to local sensitivity analyses that only compute the sensitivity of the model response at a single point in in the parameter space.

A more specific way to describe Sobol indices would probably be "variance-based sensitivity analysis".

Sobol indices are computed with respect to the distributions of the input parameters as specified in the analysis. First order Sobol indices can be thought of as a pie chart of the total variance: the total pie is the model response's variance, and each slice (first order index) is the fraction of the variance that can be attributed to variations of a single input parameter alone. If there are interactions between parameters, then the first order Sobol indices won't explain the total variance of the model's response and we get one extra slice: "variance that is unexplained by first order Sobol indices". This is because Sobol indices are not a one-at-a-time sensitivity measure, meaning they account for interactions between input parameters.

When the distribution of the input parameters is changed, the relative effect of the sensitivity inputs changes, and as a result the Sobol indices will be different. Similarly, when input parameters are swapped out / added / removed from the analysis, then the relative effect of variations of the sensitivity inputs (and thuse the Sobol indices) changes as well. Additionally, the analysis now computes interactions between different input parameters as in the original analysis, which also changes the Sobol indices.

Best,

-Florian

[Shaun Tan]

Thanks for the explanations Florian, you've been a great help!!