Main Content

This example shows how to use sensitivity analysis to narrow down the number of parameters that you need to estimate when fitting a model. This example uses a model of the vestibulo-ocular reflex, which generates compensatory eye movements.

The vestibulo-ocular reflex (VOR) enables the eyes to move at the same speed and in the opposite direction as the head, so that vision is not blurred when the head moves during normal activity. For example, if the head turns to the right, the eyes turn to the left at the same speed. This happens even in the dark. In fact, the VOR is most easily characterized by measurements in the dark, to ensure that eye movements are predominantly driven by the VOR.

Head rotation is sensed by organs in the inner ears, known as semicircular canals. These detect head motion and transmit signals about head motion to the brain, which sends motor commands to the eye muscles, so that eye movements compensate for head motion. We would like to use eye movement data to estimate model parameters for these various stages. The model we will use is shown below. There are four parameters in the model: `Delay`

, `Gain`

, `Tc`

, and `Tp`

.

```
open_system('sdoVOR')
```

The file `sdoVOR_Data.mat`

contains uniformly sampled data of stimulation and eye movements. If the VOR were perfectly compensatory, then a plot of eye movement data, when flipped vertically, would overlay exactly on top of a plot of head motion data. Such a system would be described by a gain of 1 and a phase of 180 degrees. However, real eye movements are close, but not perfectly compensatory.

load sdoVOR_Data.mat; % Column vectors: Time HeadData EyeData

We will use Sensitivity Analysis UI to see how well the model output fits the data, and explore which model parameters have the most influence on the goodness of fit. To open Sensitivity Analysis UI, in the **Apps** tab, click **Sensitivity Analyzer** under **Control Systems** to launch the **Sensitivity Analyzer**.

To associate the data with the model, click **New Requirement** and select a Signal Matching requirement. This specifies an objective function consisting of the sum of squared error between the data and model output. In the Signal Matching dialog, specify the output as `[Time EyeData]`

, and specify the input as `[Time HeadData]`

.

To view the eye movement data, navigate to the data browser on the left side of the UI, right-click the `SignalMatching`

requirement, and select Plot & Simulate. The bottom plot shows the stimulation, consisting of a series of pulses. The top plot shows eye movement data, which resembles but does not exactly match the stimulation. It also shows that the model simulated output does not match the eye movement data, because model parameters need to be estimated.

The model attempts to capture the phenomena which cause the difference between head movements and eye movements. Here we will explore the design space formed by the model parameters. To specify the parameters to explore in Sensitivity Analysis UI, click **Select Parameters** and create a new parameter set. Select all model parameters: `Delay`

, `Gain`

, `Tc`

and `Tp`

.

Explore the design space by generating parameter values. Click **Generate Values** and select random values. For repeatability of the example, reset the random number generator.

```
rng('default')
```

Since there are 4 parameters, we will generate 40 samples.

The `Delay`

parameter models the fact that there is some delay in communicating the signals from the inner ear to the brain and eyes. This delay is due to the time needed for chemical neurotransmitters to traverse the synaptic clefts between nerve cells. Based on the number of synapses in the vestibulo-ocular reflex, this delay is expected to be around 5 ms. We will model it with a uniform distribution with a lower bound of 2 ms and an upper bound of 9 ms.

The `Gain`

parameter models the fact that in the dark, the eyes do not move quite as much as the head. We will model it with a uniform distribution with a lower bound of 0.6 and an upper bound of 1.

The `Tc`

parameter models the dynamics associated with the semicircular canals, as well as some additional neural processing. The canals are high-pass filters, because after a subject has been put into rotational motion, the neurally active membranes in the canals slowly relax back to resting position, so the canals stop sensing motion. Thus, after the stimulation undergoes transition edges, eye movements tend to depart from the stimulation over time. Based on mechanical characteristics of the canals, combined with additional neural processing which prolongs this time constant to improve the accuracy of the VOR, we will model `Tc`

with a normal (i.e., bell curve) distribution with a mean of 15 seconds and a standard deviation of 3 seconds.

Finally, the `Tp`

parameter models the dynamics of the oculomotor plant, i.e. the eye and the muscles and tissues attached to it. The plant can be modeled by two poles, however it is believed that the pole with the larger time constant is cancelled by precompensation in the brain, to enable the eye to make quick movements. Thus in the plot, when the stimulation undergoes transition edges, the eye movements follow with only a little delay. We will model `Tp`

with a uniform distribution with a lower bound of 0.005 seconds and an upper bound of 0.05 seconds.

When the sample values are generated, they appear in a table in the Sensitivity Analysis UI. To plot them, select `ParamSet`

in the data browser, click the **Plots** tab, and make a scatter plot. The sampling above used default options, and these are reflected in the scatter plot. For parameters modeled by a uniform distribution, the histograms appear approximately uniform. However, parameter `Tc`

was modeled by a normal distribution, and its histogram has a bell curve profile. If Statistics and Machine Learning Toolbox™ is available, many other distributions may be used, and sampling can be done using Sobol or Halton low-discrepancy sequences. The off-diagonal plots show scatter plots between pairs of different variables. Since we did not specify cross-correlations between parameters, the scatter plots appear uncorrelated. However, if parameters were believed to be correlated, this can be specified using Correlation Matrix tab in the dialog for generating random parameter values.

Now that we have generated values for the parameter set and specified a requirement (`SignalMatching`

), we can evaluate the model. In the **Sensitivity Analysis** tab, click **Evaluate Model**.

The model is run once for each set of parameter values, and the results scatter plot is updated as new computations become available. Evaluation could also be sped up using parallel computing. After evaluation is complete, all results are also displayed in a table.

From the scatter plot of evaluation results, the `SignalMatching`

requirement seems to vary systematically as a function of `Gain`

and `Tc`

, but not `Delay`

or `Tp`

. Something similar can be seen in a contour plot. Select the `EvalResults`

variable in the data browser, click the **Plots** tab, and make a contour plot. The requirement does not vary systematically from left to right as a function of `Delay`

, but it does vertically as a function of `Gain`

.

We can use statistical analysis to quantify how much each parameter influences the requirement. Click the **Statistics** tab and select both correlation and standardized regression; and both linear and ranked analysis types. If Statistics and Machine Learning Toolbox is available, partial correlation and Kendall correlation can also be selected. Click **Compute Statistics** to carry out the calculations and show a tornado plot. The tornado plot displays results from top to bottom in order of which parameter most influences the requirement. The statistical values range from -1 to 1, where the magnitude indicates how much the parameter influences the requirement, and the sign indicates whether an increase in the parameter value corresponds to an increase or decrease in the requirement value. By most measures, this `SignalMatching`

requirement is more sensitive to `Gain`

and `Tc`

, and less sensitive to `Delay`

and `Tp`

.

For parameter estimation, we need to specify starting values for the parameters. Click the evaluation results table and click the `SignalMatching`

column header to sort results. Select the row of parameter values that minimizes the `SignalMatching`

requirement. Right-click on the row and extract these parameter values. A new variable, `ParamValues`

, is shown in the data browser.

To transition from sensitivity analysis to parameter estimation, navigate to the **Sensitivity Analysis** tab, click **Optimize**, and open a parameter estimation session. In the dialog that appears, specify that you want to use the parameter values in `ParamValues`

, and the `SignalMatching`

requirement.

Since we found above that parameters `Gain`

and `Tc`

have the most influence on the value of `SignalMatching`

, we would like to estimate only these two parameters, since the time for estimation increases with the number of parameters being estimated. In Parameter Estimation UI, click **Select Parameters** and select only `Gain`

and `Tc`

for estimation.

Since the experiment definition has been imported from `SignalMatching`

and the parameter values have been imported from `ParamValues`

, we have everything needed for estimation. Click **Estimate** to carry out parameter estimation for `Gain`

and `Tc`

. Because we are only estimating the two most influential parameters, estimation converges quickly and the model output closely matches the data. As was the case with model evaluations in sensitivity analysis, parallel computing could be used to speed up estimation.

In summary, Sensitivity Analysis UI was used to explore the parameter design space and determine that two parameters, `Gain`

and `Tc`

, were substantially more influential than the others. A start point for estimation was also determined. This start point and the requirement of obtaining a good fit to experimental data were imported into Parameter Estimation UI. Estimation completed quickly because only two parameters needed to be estimated, and the model output fit the data with very little residual error.

Close the model.

```
bdclose('sdoVOR')
```