Identify Key Parameters for Estimation (GUI)

Model Description

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.

Explore the Design Space

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.

Evaluate the Model

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.

Statistical Analysis

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.

Select Parameters for Estimation

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')