## Find Important Parameters for Receptor Occupancy with Global Sensitivity Analysis Using SimBiology Model Analyzer

This example shows how to identify important parameters in a target-mediated drug disposition
(TMDD) model using the **Global Sensitivity Analysis** (GSA) program in
**SimBiology Model Analyzer**. In this example, you compute Sobol indices and elementary
effects, and also perform multiparametric GSA (MPGSA) to find model parameters that the target or
receptor occupancy (*TO*) is sensitive to. For more information on supported
sensitivity analysis features, see Sensitivity Analysis in SimBiology.

The GSA program requires Statistics and Machine Learning Toolbox™.

### Load TMDD Model

At the MATLAB^{®} command line, enter:

`simBiologyModelAnalyzer("tmdd_with_TO.sbproj")`

The app opens, and the **TMDD** model is loaded in the
**Models** folder of the **Browser** panel.

### Perform Parameter Scan to Assess Model Behavior

Before performing the global sensitivity analysis on the model, run a parameter scan to get
a general idea of how the model response (*TO*) behaves with respect to
parameters of interest. For this example, scan the following model parameters:
*km*, *kdeg*, and *kon*.

Select

**Program**>`Run Scan`

.Go to the

**Generate Samples**step of the program. Disable auto plot generation by clicking the plot icon.In the

**Parameter Set**table, double-click the empty cell in the**Component Name**column. Enter*km*. Change**Type**,**Spacing**,**Min**,**Max**, and**# Of Steps**options as shown next.Set up

*kdeg*and*kon*similarly to km, except use 0.5 as the**Max**value.Change

**Parameter Combination**to`cartesian`

. For details about the combination methods, see Combine Simulation Scenarios in SimBiology.Click

**Simulation Settings**on the**Home**tab of the app.In the

**Simulation Time**section, change**StopTime**to 2.In the

**States To Log**section, click**Clear All**. Then select**TO**only.

Click

**Close**.Click

**Run**on the**Home**tab.The program then opens the

**Plot1**tab. Click**Time**in the**Property Editor**pane on the right to see the time courses of the model response (*TO*) with respect to different parameter combinations of*km*,*kdeg*, and*kon*.Categorize the plot to make it easier to identify the sources of sensitivity. In the

**Property Editor**pane, under**Visual Channels**, change the style of*km*,*kdeg*, and*kon*as in the following figure. You need to also change the style of**Scenarios**from`Color`

to empty.Click

**Plot Settings**and click**Link Y-Axis**to use the same Y-axis limit for all plots and to make the plots more comparable. The plot indicates that*TO*is sensitive to changes in*kon*and*km*values. It seems to be less sensitive to variations in*kdeg*. Optionally, you can also scan one parameter at a time by disabling and enabling the parameters in the**Parameter Set**table to see the effect of each parameter.Follow the next steps to perform GSA on these parameters to get more insights into relative contributions of individual parameters that contribute most to the overall model behavior.

### Compute Sobol Indices

The section shows you how to perform GSA by computing the first-order and total-order Sobol indices of model
parameters using the **Global Sensitivity Analysis** program.

Select

**Program**>`Global Sensitivity Analysis`

. By default, the program uses**Sobol indices**as a GSA analysis method.Define the model parameters of interest in the

**Sensitivity Inputs**table. Double-click the empty cell under**Component Name**column and enter*km*.By default, the program uses the uniform distribution to sample parameter values. To see a list of supported distributions, click

**Uniform**. For this example, use the uniform distribution and set the**Lower**and**Upper**bounds to`1e-3`

and 0.1, respectively.Similarly, set up

*kdeg*and*kon*as sensitivity inputs. Use 0.5 as the**Upper**value.In

**Sampling Options**, use the default values.In the

**Sensitivity Outputs**table, enter*TO*.Click the Run button of the

**Sobol Indices**step.The status bar at the bottom of the program provides the run progress and the number of runs that have been finished.

The program opens the

**Datasheet1**and**Plot2**tabs. By default,**Plot2**shows the**GSA Time**plot of the first-order and total-order Sobol indices for each input parameter.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. Based on the first-order and total-order plots, both

*km*and*kon*seem to be the sensitive parameters to*TO*.*km*seems to become more sensitive at later time points while*kon*is more sensitive before time = 1. The fraction of unexplained variance plot is a flat line, which means that there is almost no unexplained variance and most variances are accounted for in the first-order and total-order plots.Click the

**Datasheet1**tab to view the tables of GSA program results. The first table shows the time dependent responses of Sobol indices. In the second table, you can find out if any model simulation failed during the computation by checking the**Number of Successful Samples**row. In this example, there were no failed simulations as all 1000 samples were successfully simulated.

### Compute Elementary Effects

The following steps show you how to reuse the same program setup in the previous step to
compute the means and standard deviations of the elementary effects of
input parameters with respect to the model response *TO*.

Given that computing GSA results is often computationally expensive, save the previous Sobol indices results before switching over to elementary effects. In the

**Browser**panel, expand**Program2**. Right-click**LastRun**and select**Save Data**. Enter the name of the data as`sobol_results`

. Click**OK**.Go back to the

**Program2**tab. In the**Sobol Indices**step, under the**Analysis**section, select`Sobol indices`

>`Elementary effects`

.The

**Sensitivity Inputs**table is automatically set up with the same set of parameters and their corresponding upper and lower bounds as the previous step.*TO*is the only sensitivity output.Use the default values for

**Grid Settings**. Click the info icon to display additional information for each option. For details about how the algorithm uses these settings, see Elementary Effects for Global Sensitivity Analysis.Use the default value (

`AbsoluteEffects = true`

) for the**Output Settings**as well. The program uses the absolute elementary effects by default because the elementary effects can average out when you calculate the mean otherwise.Click the Run button of the

**Elementary Effects**step.The program then opens the

**Datasheet2**and**Plot3**tabs. By default,**Plot3**shows the**GSA Time**plot of the mean and standard deviation of elementary effects for each input parameter.The mean of elementary effects explains whether variations in input parameter values have any effect on the model response

*TO*. The mean plots indicate*km*becomes more sensitive at later time points and*kon*is more sensitive before time = 0.5. This trend is similar to the Sobol GSA results as well.*kdeg*also shows some sensitivity but might be insignificant because mean values are much smaller than those of*km*and*kon*. The standard deviation of effects explains whether the sensitivity change is dependent on the location in the parameter domain. For instance, the*km*standard deviation plot indicates larger deviations for larger parameter values in the later stage (time > 0.5).Click the

**Datasheet2**tab to view the tables of GSA program results. The first table shows the time dependent responses of elementary effects. In the second table, you can find out if any model simulation failed during the computation by checking the**Number of Successful Samples**row. In this example, there were no failed simulations as all 1000 samples were successfully simulated.

### Perform Multiparametric GSA (MPGSA)

The following steps show you how to reuse the same program set up in the previous step to
perform multiparametric global
sensitivity analysis to answer the question of whether the input parameters have any
effects on the exposure (area under the curve of the *TO* profile) threshold
for the target occupancy.

Save the previous elementary effects GSA results before switching over to MPGSA. In the

**Browser**panel, expand**Program2**. Right-click**LastRun**and select**Save Data**. Enter the name of the data as`ee_results`

. Click**OK**.Click the

**Program2**tab. In the**Elementary Effects**step, under the**Analysis**section, select`MPGSA (multi-parametric global sensitivity analysis)`

.The stop time is set to 2 and the

**Sensitivity Inputs**table is automatically set up with the same number of samples (1000) and the same set of parameters, with their corresponding upper and lower bounds as the previous GSA analysis.Use default settings for

**Sampling Options**and**Correlation Matrix**.In the

**Classifiers**table, double-click the empty cell in the**Expression**column and enter:`trapz(time,TO) <= 0.5`

. This expression defines an exposure (area under the curve of the*TO*profile) threshold for the target occupancy.Use the default value (

`0.05`

) for**Significance Level**.Click the Run button of the

**MPGSA**step. You might see a warning about dimensional analysis not being able to perform. For the purposes of this example, you can ignore the warning.**Tip**The GSA program enables you to perform

**MPGSA**as a standalone analysis or as an added step following the**Sobol Indices**or**Elementary Effects**analysis. A benefit of adding**MPGSA**as a subsequent step is that it can reuse the simulation results from the previous step if possible and saves computation time. To add MPGSA as a subsequent step, click the green plus icon at the top of the GSA program and click`MPGSA`

.The MPGSA subsequent step reuses the simulation results whenever the model response (such as

**TO**) defined in the classifier was used as the sensitivity output in the previous step.

The program then opens

**Datasheet3**with GSA results and**Plot4**showing the bar plot of the MPGSA results for each input parameter. To see the empirical cumulative distribution functions (eCDFs) curves of the results, click**eCDF**in**Plot Settings**.For each parameter, two curves show eCDFs of the accepted and rejected samples. Except for

*km*, the other two parameters do not seem to show a significance difference between the two curves. These plots qualitatively show that*km*and*kon*affect the classification of samples while*kdeg*does not. Quantitatively, the maximum absolute distance between two eCDFs curves is called the Kolmogorov-Smirnov (K-S) distance. The*km*plot shows a large K-S distance.To compute the K-S distance between the two eCDFs, SimBiology uses a two-sided test based on the null hypothesis that the two distributions of accepted and rejected samples are equal. See

`kstest2`

(Statistics and Machine Learning Toolbox) for details. If the K-S distance is large, then the two distributions are different, meaning that the classification of the samples is sensitive to variations in the input parameter. On the other hand, if the K-S distance is small, then variations in the input parameter do not affect the classification of samples.To assess the significance of the K-S statistic rejecting the null-hypothesis, you can examine the p-values by looking at the bar plot. Click

**Bar**in the**Plot Settings**.The bar plot shows two bars for each parameter: one for the K-S distance (K-S statistic) and another for the corresponding p-value. You reject the null hypothesis if the p-value is less than the significance level. A cross (x) is shown for any p-value that is almost 0.

The p-values for

*km*and*kon*are less than the significance level (0.05), supporting the alternative hypothesis that the accepted and rejected samples come from different distributions. The p-value of*kdeg*is larger than the significance level supporting the null hypothesis that the accepted and rejected samples come from the same distribution. To conclude, the classification of the samples is sensitive to*km*and*kon*, not*kdeg*. These results are in agreement with the previous GSA results which identify*km*and*kon*as sensitive parameters but not*kdeg*.To see the exact p-value corresponding to each parameter, click

**Datasheet3**.

## See Also

SimBiology Model Analyzer | `sbiosobol`

| `sbioelementaryeffects`

| `sbiompgsa`

| `ecdf`

(Statistics and Machine Learning Toolbox) | `kstest2`

(Statistics and Machine Learning Toolbox)