Main Content

Perform global sensitivity analysis by computing first- and total-order Sobol indices (requires Statistics and Machine Learning Toolbox)

performs global sensitivity analysis [1] on a SimBiology model
`sobolResults`

= sbiosobol(`modelObj`

,`params`

,`observables`

)`modelObj`

by decomposing the variances of
`observables`

with respect to the sensitivity inputs
`params`

.

uses additional options specified by one or more name-value pair arguments.`sobolResults`

= sbiosobol(`modelObj`

,`params`

,`observables`

,`Name,Value`

)

Load the Tumor Growth Model.

`sbioloadproject tumor_growth_vpop_sa.sbproj`

Get a variant with the estimated parameters and the dose to apply to the model.

```
v = getvariant(m1);
d = getdose(m1,'interval_dose');
```

Get the active configset and set the tumor weight as the response.

```
cs = getconfigset(m1);
cs.RuntimeOptions.StatesToLog = 'tumor_weight';
```

Simulate the model and plot the tumor growth profile.

sbioplot(sbiosimulate(m1,cs,v,d));

Perform global sensitivity analysis (GSA) on the model to find the model parameters that the tumor growth is sensitive to.

First, retrieve model parameters of interest that are involved in the pharmacodynamics of the tumor growth. Define the model response as the tumor weight.

modelParamNames = {'L0','L1','w0','k1','k2'}; outputName = 'tumor_weight';

Then perform GSA by computing the first- and total-order Sobol indices using `sbiosobol`

. Set `'ShowWaitBar'`

to `true`

to show the simulation progress. By default, the function uses 1000 parameter samples to compute the Sobol indices [1].

sobolResults = sbiosobol(m1,modelParamNames,outputName,'Variants',v,'Doses',d,'ShowWaitBar',true)

sobolResults = Sobol with properties: Time: [444x1 double] SobolIndices: [5x1 struct] Variance: [444x1 table] Observables: {'[Tumor Growth Model].tumor_weight'} ParameterSamples: [1000x5 table] SimulationInfo: [1x1 struct]

You can change the number of samples by specifying the `'NumberSamples'`

name-value pair argument. The function requires a total of `(number of input parameters + 2) * NumberSamples`

model simulations.

Show the mean model response, the simulation results, and a shaded region covering 90% of the simulation results.

plotData(sobolResults);

You can adjust the quantile region to a different percentage by specifying `'Alphas' `

for the lower and upper quantiles of all model responses. For instance, an alpha value of 0.1 plots a shaded region between the `100 * alpha`

and `100 * (1 - alpha)`

quantiles of all simulated model responses.

`plotData(sobolResults,'Alphas',0.1);`

Plot the time course of the first- and total-order Sobol indices.

```
h = plot(sobolResults);
% Resize the figure.
h.Position(:) = [100 100 1280 800];
```

The first-order Sobol index of an input parameter gives the fraction of the overall response variance that can be attributed to variations in the input parameter alone. The total-order index gives the fraction of the overall response variance that can be attributed to any joint parameter variations that include variations of the input parameter.

From the Sobol indices plots, parameters `L1`

and `w0`

seem to be the most sensitive parameters to the tumor weight before the dose was applied at t = 7. But after the dose is applied, `k1`

and `k2`

become more sensitive parameters and contribute most to the after-dosing stage of the tumor weight. The total variance plot also shows a larger variance for the after-dose stage at t > 35 than for the before-dose stage of the tumor growth, indicating that `k1`

and `k2`

might be more important parameters to investigate further. The fraction of unexplained variance shows some variance at around t = 33, but the total variance plot shows little variance at t = 33, meaning the unexplained variance could be insignificant. The fraction of unexplained variance is calculated as 1 - (sum of all the first-order Sobol indices), and the total variance is calculated using `var(response)`

, where `response`

is the model response at every time point.

You can also display the magnitudes of the sensitivities in a bar plot.

bar(sobolResults)

You can specify more samples to increase the accuracy of the Sobol indices, but the simulation can take longer to finish. Use `addsamples`

to add more samples. For example, if you specify 1500 samples, the function performs `1500 * (2 + number of input parameters)`

simulations.

gsaMoreSamples = addsamples(gsaResults,1500)

The SimulationInfo property of the result object contains various information for computing the Sobol indices. For instance, the model simulation data (SimData) for each simulation using a set of parameter samples is stored in the `SimData`

field of the property. This field is an array of `SimData`

objects.

sobolResults.SimulationInfo.SimData

SimBiology SimData Array : 1000-by-7 Index: Name: ModelName: DataCount: 1 - Tumor Growth Model 1 2 - Tumor Growth Model 1 3 - Tumor Growth Model 1 ... 7000 - Tumor Growth Model 1

You can find out if any model simulation failed during the computation by checking the `ValidSample`

field of `SimulationInfo`

. In this example, the field shows no failed simulation runs.

all(sobolResults.SimulationInfo.ValidSample)

`ans = `*1x7 logical array*
1 1 1 1 1 1 1

`SimulationInfo.ValidSample`

is a table of logical values. It has the same size as `SimulationInfo.SimData`

. If `ValidSample`

indicates that any simulations failed, you can get more information about those simulation runs and the samples used for those runs by extracting information from the corresponding column of `SimulationInfo.SimDat`

a. Suppose that the fourth column contains one or more failed simulation runs. Get the simulation data and sample values used for that simulation using `getSimulationResults`

.

[samplesUsed,sd,validruns] = getSimulationResults(sobolResults,4);

You can add custom expressions as observables and compute Sobol indices for the added observables. For example, you can compute the Sobol indices for the maximum tumor weight by defining a custom expression as follows.

% Suppress an information warning that is issued during simulation. warnSettings = warning('off', 'SimBiology:sbservices:SB_DIMANALYSISNOTDONE_MATLABFCN_UCON'); % Add the observable expression. sobolObs = addobservable(sobolResults,'Maximum tumor_weight','max(tumor_weight)','Units','gram');

Plot the computed simulation results showing the 90% quantile region.

h2 = plotData(sobolObs); h2.Position(:) = [100 100 1280 800];

You can also remove the observable by specifying its name.

`gsaNoObs = removeobservable(sobolObs,'Maximum tumor_weight');`

Restore the warning settings.

warning(warnSettings);

`modelObj`

— SimBiology modelSimBiology model object

SimBiology model, specified as a SimBiology `model object`

.

`params`

— Names of model parameters, species, or compartmentscharacter vector | string | string vector | cell array of character vectors

Names of model parameters, species, or compartments, specified as a character vector, string, string vector, or cell array of character vectors.

**Example: **`["k1","k2"]`

**Data Types: **`char`

| `string`

| `cell`

`observables`

— Model responsescharacter vector | string | string vector | cell array of character vectors

Model responses, specified as a character vector, string, string vector, or cell
array of character vectors. Specify the names of species, parameters, compartments, or
`observables`

.

**Example: **`["tumor_growth"]`

**Data Types: **`char`

| `string`

| `cell`

Specify optional
comma-separated pairs of `Name,Value`

arguments. `Name`

is
the argument name and `Value`

is the corresponding value.
`Name`

must appear inside quotes. You can specify several name and value
pair arguments in any order as
`Name1,Value1,...,NameN,ValueN`

.

```
sobolResults =
sbiosobol(modelObj,params,observables,'ShowWaitbar',true)
```

specifies to show a
simulation progress bar.`'Bounds'`

— Parameter boundsnumeric matrix

Parameter bounds, specified as the comma-separated pair consisting of
`'Bounds'`

and a numeric matrix
with two columns. The first column contains the lower bounds and
the second column contains the upper bounds. The number of rows
must be equal to the number of parameters in
`params`

.

If a parameter has a nonzero value, the default bounds are ±10% of the value. If the
parameter value is zero, the default bounds are `[0 1]`

.

**Example: **`[0.5 5]`

**Data Types: **`double`

`'Doses'`

— Doses to use during simulations`ScheduleDose`

object | `RepeatDose`

object | vector of dose objectsDoses to use during model simulations, specified as the comma-separated pair consisting of
`'Doses'`

and a `ScheduleDose`

or
`RepeatDose`

object or a vector of dose objects.

`'Variants'`

— Variants to apply before simulationsvariant object | vector of variant objects

Variants to apply before model simulations, specified as the comma-separated pair consisting
of `'Variants'`

and a variant object or vector of variant
objects.

When there are multiple variants and if there are duplicate specifications for a property's value, the last occurrence for the property value in the array of variants is used during simulation.

`'NumberSamples'`

— Number of samples to compute Sobol indices`1000`

(default) | positive integerNumber of samples to compute Sobol indices, specified as the comma-separated pair
consisting of `'NumberSamples'`

and a positive integer. The function
requires `(number of input `

model simulations to compute the first- and total-order
Sobol indices.`params`

+ 2) *
NumberSamples

**Data Types: **`double`

`'SamplingMethod'`

— Method to generate parameter samples`'Sobol'`

(default) | character vector | stringMethod to generate parameter samples, specified as the comma-separated pair consisting of
`'SamplingMethod'`

and a character vector or string. Valid options
are:

`'Sobol'`

— Use the low-discrepancy Sobol sequence to generate samples.`'Halton'`

— Use the low-discrepancy Halton sequence to generate samples.`'lhs'`

— Use the low-discrepancy Latin hypercube samples.`'random'`

— Use uniformly distributed random samples.

**Data Types: **`char`

| `string`

`'StopTime'`

— Simulation stop timenonnegative scalar

Simulation stop time, specified as the comma-separated pair consisting of
`'StopTime'`

and a nonnegative scalar. If you specify neither
`StopTime`

nor `OutputTimes`

, the function uses
the stop time from the active configuration set of the model. You cannot specify both
`StopTime`

and `OutputTimes`

.

**Data Types: **`double`

`'OutputTimes'`

— Simulation output timesnumeric vector

Simulation output times, specified as the comma-separated pair consisting of
`'OutputTimes'`

and a numeric vector. The function computes the
Sobol indices at these output time points. You cannot specify both
`StopTime`

and `OutputTimes`

. By default, the
function uses the output times of the first model simulation.

**Example: **`[0 1 2 3.5 4 5 5.5]`

**Data Types: **`double`

`'UseParallel'`

— Flag to run model simulations in parallel`false`

(default) | `true`

Flag to run model simulations in parallel, specified as the comma-separated pair consisting of
`'UseParallel'`

and `true`

or
`false`

.

**Data Types: **`logical`

`'Accelerate'`

— Flag to turn on model acceleration`true`

(default) | `false`

Flag to turn on model acceleration, specified as the comma-separated pair consisting of
`'Accelerate'`

and `true`

or
`false`

.

**Data Types: **`logical`

`'InterpolationMethod'`

— Method for interpolation of model simulations`"interp1q"`

(default) | character vector | stringMethod for interpolation of model simulations to new output times, specified as the
comma-separated pair consisting of `'InterpolationMethod'`

and a
character vector or string. The valid options follow.

**Data Types: **`char`

| `string`

`'ShowWaitbar'`

— Flag to show progress of model simulations`false`

(default) | `true`

Flag to show the progress of model simulations by displaying a progress bar,
specified as the comma-separated pair consisting of `'ShowWaitbar'`

and `true`

or `false`

. By default, no wait bar is
displayed.

**Data Types: **`logical`

`sobolResults`

— Results containing Sobol indices`SimBiology.gsa.Sobol`

objectResults containing the first- and total-order Sobol indices, returned as a `SimBiology.gsa.Sobol`

object. The object also contains the parameter sample
values and model simulation data used to compute the Sobol indices.

The results object can contain a significant amount of
simulation data (SimData). The size of the object exceeds ```
(1 + number of observables) *
number of output time points * (2 + number of parameters) * number of samples * 8
```

bytes. For example, if you have one observable, 500 output time points, 8 parameters, and
100,000 samples, the object size is `(1 + 1) * 500 * (2 + 8) * 100000 * 8`

bytes = 8 GB. If you need to save such large objects, use this syntax:

`save(fileName,variableName,'-v7.3');`

`sbiosobol`

implements the Saltelli method [1] to compute Sobol
indices.

Consider a SimBiology model response *Y* expressed as a mathematical model $$Y=f\left({X}_{1},{X}_{2},{X}_{3},\mathrm{...},{X}_{k}\right)$$, where *X _{i}* is a model parameter
and

`i = 1,…,k`

.The first-order Sobol index (*S _{i}*) gives the fraction
of the overall response variance

`V(`*Y*)

that can be
attributed to variations in $${S}_{i}=\frac{{V}_{{X}_{i}}\left({E}_{{X}_{\sim i}}\left(Y|{X}_{i}\right)\right)}{V\left(Y\right)}$$

The total-order Sobol index (*S _{Ti}*) gives the fraction
of the overall response variance

`V(`*Y*)

that can be
attributed to any joint parameter variations that include variations of
$${S}_{Ti}=1-\frac{{V}_{{X}_{\sim i}}\left({E}_{{X}_{i}}\left(Y|{X}_{\sim i}\right)\right)}{V\left(Y\right)}=\frac{{E}_{{X}_{\sim i}}\left({V}_{{X}_{i}}\left(Y|{X}_{\sim i}\right)\right)}{V\left(Y\right)}$$

To compute individual values for *Y* corresponding to samples of parameters
*X _{1}*,

$$A=\left(\begin{array}{cccc}{X}_{11}& {X}_{12}& \mathrm{...}& {X}_{1k}\\ {X}_{21}& {X}_{22}& \mathrm{...}& {X}_{2k}\\ \mathrm{...}& \mathrm{...}& \mathrm{...}& \mathrm{...}\\ {X}_{n1}& {X}_{n2}& \mathrm{...}& {X}_{nk}\end{array}\right)$$

$$B=\left(\begin{array}{cccc}{X}_{11}^{\text{'}}& {X}_{12}^{\text{'}}& \mathrm{...}& {X}_{1k}^{\text{'}}\\ {X}_{21}^{\text{'}}& {X}_{22}^{\text{'}}& \mathrm{...}& {X}_{2k}^{\text{'}}\\ \mathrm{...}& \mathrm{...}& \mathrm{...}& \mathrm{...}\\ {X}_{n1}^{\text{'}}& {X}_{n2}^{\text{'}}& \mathrm{...}& {X}_{nk}^{\text{'}}\end{array}\right)$$

*n* is the sample size. Each row of the matrices
*A* and *B* corresponds to one parameter sample set,
which is a single realization of model parameter values.

Estimates for *S _{i}* and

```
i = 1, 2,
…, params
```

.$${{\displaystyle A}}_{B}^{i}=\left(\begin{array}{cccccc}{X}_{11}& {X}_{12}& \mathrm{...}& {X}_{1i}^{\text{'}}& \mathrm{...}& {X}_{1k}\\ {X}_{21}& {X}_{22}& \mathrm{...}& {X}_{2i}^{\text{'}}& \mathrm{...}& {X}_{2k}\\ \mathrm{...}& \mathrm{...}& \mathrm{...}& \mathrm{...}& \mathrm{...}& \mathrm{...}\\ {X}_{n1}& {X}_{n2}& \mathrm{...}& {X}_{ni}^{\text{'}}& \mathrm{...}& {X}_{nk}\end{array}\right)$$

The formulas to approximate the first- and total-order Sobol indices are as follows.

$${\widehat{S}}_{i}=\frac{\frac{1}{n}{\displaystyle \sum _{j=1}^{n}f{(B)}_{j}\left(f{\left({{\displaystyle A}}_{B}^{i}\right)}_{j}-f{(A)}_{j}\right)}}{V(Y)}$$

$${\widehat{S}}_{Ti}=\frac{\frac{1}{2n}{\displaystyle \sum _{j=1}^{n}{\left(f{(A)}_{j}-f{\left({{\displaystyle A}}_{B}^{i}\right)}_{j}\right)}^{2}}}{V(Y)}$$

`f(`

, *A*)`f(`

, and $$f{\left({{\displaystyle A}}_{B}^{i}\right)}_{j}$$ are the model simulation results using the parameter sample values from
matrices *B*)*A*, *B*, and $${{\displaystyle A}}_{B}^{i}$$.

The matrix *A* corresponds to the `ParameterSamples`

property of the Sobol results object (`resultsObj.ParameterSamples`

). The matrix *B* corresponds to the `SupportSamples`

property (`resultsObj.SimulationInfo.SupportSamples`

).

The $${{\displaystyle A}}_{B}^{i}$$ matrices are stored in the `SimData`

structure of the
`SimulationInfo`

property
(`resultsObj.SimulationInfo.SimData`

). The size of
`SimulationInfo.SimData`

is

, where *NumberSamples*-by-*params* +
2*NumberSamples* is the number of samples and *param* is the number of input parameters. The number of columns is
`2 + `

because the first column of
*params*`SimulationInfo.SimData`

contains the model simulation results using
the sample matrix *A*. The second column contains simulation results using
`SupportSamples`

, which is another sample matrix
*B*. The rest of the columns contain simulation results using $${{\displaystyle A}}_{B}^{1}$$, $${{\displaystyle A}}_{B}^{2}$$, …, $${{\displaystyle A}}_{B}^{i}$$, …, $${{\displaystyle A}}_{B}^{params}$$. See `getSimulationResults`

to retrieve the model simulation results and samples
for a specified *i*th index ($${{\displaystyle A}}_{B}^{i}$$) from the `SimulationInfo.SimData`

array.

[1] Saltelli, Andrea, Paola Annoni, Ivano Azzini, Francesca Campolongo, Marco Ratto, and Stefano Tarantola. “Variance Based Sensitivity Analysis of Model Output. Design and Estimator for the Total Sensitivity Index.” *Computer Physics Communications* 181, no. 2 (February 2010): 259–70. https://doi.org/10.1016/j.cpc.2009.09.018.

You have a modified version of this example. Do you want to open this example with your edits?

You clicked a link that corresponds to this MATLAB command:

Run the command by entering it in the MATLAB Command Window. Web browsers do not support MATLAB commands.

Choose a web site to get translated content where available and see local events and offers. Based on your location, we recommend that you select: .

Select web siteYou can also select a web site from the following list:

Select the China site (in Chinese or English) for best site performance. Other MathWorks country sites are not optimized for visits from your location.

- América Latina (Español)
- Canada (English)
- United States (English)

- Belgium (English)
- Denmark (English)
- Deutschland (Deutsch)
- España (Español)
- Finland (English)
- France (Français)
- Ireland (English)
- Italia (Italiano)
- Luxembourg (English)

- Netherlands (English)
- Norway (English)
- Österreich (Deutsch)
- Portugal (English)
- Sweden (English)
- Switzerland
- United Kingdom (English)