Computational and Mathematical Methods in MedicineVolume 2012 (2012), Article ID 580795, 5 pagesdoi:10.1155/2012/580795

Research Article

Multiscale Autoregressive Identification of Neuroelectrophysiological Systems

Timothy P. Gilmour,¹ Thyagarajan Subramanian,² Constantino Lagoa,¹ and W. Kenneth Jenkins¹

¹Electrical Engineering Department, Pennsylvania State University, University Park, PA 16802, USA²Neurology Department, Penn State Hershey Medical Center, 500 University Drive, Hershey, PA 17033, USA

Received 8 September 2011; Revised 1 December 2011; Accepted 2 December 2011

Academic Editor: Henggui Zhang

Copyright © 2012 Timothy P. Gilmour et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract.

Electrical signals between connected neural nuclei are difficult to model because of the complexity and high number of paths within the brain. Simple parametric models are therefore often used. A multiscale version of the autoregressive with exogenous input (MS-ARX) model has recently been developed which allows selection of the optimal amount of filtering and decimation depending on the signal-to-noise ratio and degree of predictability. In this paper, we apply the MS-ARX model to cortical electroencephalograms and subthalamic local field potentials simultaneously recorded from anesthetized rodent brains. We demonstrate that the MS-ARX model produces better predictions than traditional ARX modeling. We also adapt the MS-ARX results to show differences in internuclei predictability between normal rats and rats with 6OHDA-induced parkinsonism, indicating that this method may have broad applicability to other neuroelectrophysiological studies.

1. Introduction

VARIOUS types of methods have been used to assess the degree of similarity or shared information between two signals. The methods used depend on the type of the presumptive system which processes the one “input” signal into the other “output” signal. Two basic classifications of systems are whether they are memoryless or not, and whether they are linear or not.

Common linear memoryless methods include the cross-correlation in the time domain or coherence in the frequency domain. Common linear models with memory include autoregressive models with exogenous input (ARX), autoregressive moving average (ARMA), Box-Jenkins, Output-Error, and linear state-space models.

Higher-order nonlinear methods with memory are also sometimes used, such as polyspectral models, nonlinear ARX models, neural networks, Hammerstein-Wiener models, and Volterra models, but these are more difficult to train.

Other statistical evaluations of the similarity focus on the transfer of information between two signals, rather than explicit modeling and prediction. Examples of these analyses include Granger causality analysis, time-delay mutual information, and transfer entropy. The information theoretic analyses can measure nonlinear as well as linear effects.

Complex systems such as the brain are difficult to analyze because of the huge number of individual neuronal/synaptic paths between nuclei, the nonlinear nature of neuronal connections, and the operation at multiple time scales.

One approach is to use simple low-order linear models to approximate the transfer function relationship, such as autoregressive with exogenous (ARX) models. The advantage of using linear ARX models is that there is no need to estimate nonlinearity parameters, and less training data is required. However, the performance of such models depends crucially on the model order, scale, and prefiltering. To mitigate this dependency, a multiscale version of the autoregressive with exogenous input (MS-ARX) model has recently been developed by Nounou and colleagues [1]. The MS-ARX model allows automatic selection of the optimal scale for the ARX prediction.

In this paper, we adapt and apply the MS-ARX model to evaluate the degree of information transfer between cortical electroencephalogram (EEG) and subthalamic nucleus (STN) local field potential (LFP) signals. In a rat model of Parkinson’s disease, the multiscale ARX approach showed significant differences in connectivity compared to normal.

2. Methods

2.1. Autoregressive System Identification

The ARX model is a common method to represent output signals from an unknown system by using a linear combination of past output signal values and past input values. We will be following the notation of Nounou and colleagues in our model description. The equation for the ARX model is where is the output, is the input, and are the estimated system coefficients, and and are the maximum orders of the autoregressive and input filters, respectively. Equation (1) may be written in matrix form as where The weight parameters and may be solved using least squares: The maximum filter lengths and may be estimated by minimizing some criterion such as the Akaike information criterion (AIC): where is the number of model parameters, and is the likelihood function quantifying the model goodness of fit.

2.1.1. Wavelet Decomposition

Signals may be decomposed into a multiscale time-frequency representation by projecting the signal onto an orthonormal set of basic functions. These functions correspond to a particular scale and translation of a prototype scaling function and wavelet function , given by For the Haar wavelet used in this paper,

2.1.2. Multiscale ARX Modeling

We applied the multiscale ARX approach presented by Nounou and colleagues. Briefly, the input (EEG) and output (LFP) data were first split in half into a training and a validation set. Second, both sets were decomposed using Haar wavelets into multiple scaled approximations in addition to the original undecimated scale. Third, at each scale an ARX model was trained using the model structure selected by an AIC minimization. Fourth, the computed ARX model from each scale was converted to the original sampling rate using the following theorem proved by Nounou et al.: an ARX transfer function at scale is equivalent to the transfer function at scale 0, the original undecimated scale. Fifth, the optimal scale ARX model was selected as the one which provided the smallest mean-square error (MSE) on the validation set (Figure 1). The MSE is defined as

Figure 1: Block diagram of MS-ARX system-identification configuration. In general, the system (dotted box) is unknown and so the true output and measurement noise are unknown and only is measurable. The multirate equivalence theorem stated by Nounou et al. allows precomputation of the scaled wavelet ARX prediction blocks (dashed boxes), reducing computational complexity.

2.1.3. Neural Data Collection

The motor cortex has been shown to project into the subthalamic nucleus (STN), thus implying a system of unknown electrical parameters with the EEG as the input and the STN local field potential (LFP) as the output [2]. Furthermore, studies in preclinical models of Parkinson’s disease have shown increased correlation between neighboring neurons [3] and increased coherence in the 15–30 Hz band between basal ganglia nuclei [4, 5].

To examine the connection strength between these disparate brain areas using MS-ARX, we simultaneously recorded voltage data from the motor cortex EEG and STN LFP of anesthetized normal rats and hemiparkinsonian (HP) rats. Parkinsonism was induced by the vendor (Charles River) by 6-hydroxydopamine injection (12 μg in 4 μL, injected into coordinates AP −1.5 mm, ML +1.8, DV −7.5 from dura at 0.67 μL/minute, cf. [6, 7]) and was verified by apomorphine-induced rotation testing and histological verification as detailed elsewhere [7, 8]. Briefly, rats were injected with apomorphine HCl (0.2 mg/kg, subcutaneously) at 3 and 5 weeks after 6-OHDA exposure and the number of contralateral turns over 35 minutes were counted with an automated rotameter. Rats averaging more than 7 turns per minute were included in the HP group. After recording and euthanasia, the brain was frozen, sectioned coronally, stained for tyrosine hydroxylase to confirm 95% unilateral lesioning of the substantia nigra pars compacta, and stained with cresyl violet to confirm electrode localization. All procedures were approved by the Pennsylvania State University Institutional Animal Care and Use Committee.

For EEG recording, animals were deeply anesthetized with urethane, with nominal initial dose 1.3 g/kg (i.p.) and additional doses given as needed to maintain surgical anesthesia. Stainless steel screws were implanted above bregma and above motor cortex (AP  mm, ML +2.5), and the EEG signal recorded as the potential difference between these screws was subsequently amplified and filtered between 0.1 Hz and 500 Hz (3500, A-M Systems) and digitized at an initial rate of 1000 samples per second. An occipital screw was used as the reference electrode. The LFP signal was taken from the tip of the tungsten microelectrode (1-2 MΩ, FHC Inc), filtered between 5 Hz–500 Hz, amplified, and digitized. Offline, both signals were filtered between 0.1 Hz and 100 Hz and downsampled to 200 samples per second.

Recordings were taken from 9 normal rats (37 distinct STN sites) and 8 HP rats (26 distinct STN sites). Electrode tracts were histologically confirmed.

We used a maximum of 260 taps in our AIC structure selection step. This maximum was empirically selected based on the observation that the peak frequencies in our data were usually between 0.8–1.3 Hz or higher, thus allowing at least one cycle period within the ARX filter length at scale zero. The maximum wavelet decomposition level was 4. Only recordings with robust slow-wave activity were included in the analysis [9].

3. Results

Table 1 shows the MSEs at different scales. The undecimated scale was the optimal scale for approximately half of the recordings. The other recordings saw better ARX prediction performance at higher wavelet scales (more heavily filtered wavelet approximations).

Table 1: Mean MSE at different scales across all recordings. Numbers in parentheses are the percent of neuronal recordings which selected that particular scale as optimum.

Scale
MSE (mean ± SEM)

Normal
HP

(51%)
(56%)

(11%)
(24%)

(3%)
(4%)

(8%)
(8%)

(27%)
(8%)

The mean MSE at the optimum scale was significantly lower in the HP group compared to the normal group (Figure 2, , rank-sum test). Also, the ratio (computed for each individual recording) of the mean of the absolute value of the best-scale AR coefficients () to the mean of the absolute value of the best-scale exogenous input coefficients () was significantly lower in the HP group (Table 2, , rank-sum test). Figures 3 and 4 show samples of the wavelet decimation and MS-ARX predictions, illustrating the differences between scales.

Table 2: Parameter summary for normal and hemiparkinsonian recordings. Each row shows the mean (± SEM) absolute value of the parameter at the optimum scale for each recording (*denotes rank-sum test between Normal and HP groups).

Mean of parameter
Normal
HP

AR coeffs.

Exogenous input (X) coeffs.

Ratio of AR to X coeffs.

Number of nonzero AR coeffs.

Number of nonzero X coeffs.

Best scale

Figure 2: Mean MSE across all recordings at their optimal scale (*denotes rank-sum test between Normal and HP groups).

(a)

(b)

(c)

(d)

Figure 3: Sample EEG and LFP waveforms showing the original scale and three successive levels of scaled wavelet decimations.

(a)

(b)

Figure 4: (a) Overlaid sample EEG input, LFP output, and the MS-ARX optimal scale prediction (scale 3). (b) Overlaid LFP output and predictions from all scales.

4. Discussion

The MS-ARX technique showed improved prediction accuracy compared to the traditional ARX approach (which uses scale 0 only). This makes sense because the MS-ARX approach uses cross-validation to automatically select the best tradeoff between smoothing and preservation of signal details.

The results seen of decreased MSE and increased proportion of exogenous input weights also indicate that linear prediction of the STN LFP based on the cortical EEG is more accurate in the HP condition and is based more on cortical input. This indicates a greater amount of similarity between the population-based cortical and STN electrical activity in the HP case. Future studies should investigate whether this increased predictability is correlated to parkinsonian pathophysiology or symptomatic severity. Brown and colleagues have shown that there is increased coherence between the cortical EEG and STN LFP in the 15–30 Hz “beta” band in PD [4, 5, 10, 11], and also that STN LFP beta activity is correlated to clinical symptoms [12–14], suggesting that the linear predictability in those scales may also be correlated to clinical symptoms. However, unlike the coherence, the MS-ARX prediction measure we describe can also measure nonperiodic content similarity between the cortical EEG and STN LFP. Thus it may provide an additional useful tool to investigate nonperiodic corticosubthalamic interactions. The recent optogenetic study by Deisseroth and colleagues showed that high-frequency stimulation of the motor cortex achieved similar symptomatic amelioration as STN stimulation, presumably affecting the STN through the hyperdirect pathway [2, 15]. However, this effect was only seen in high-frequency stimulation, suggesting that the temporal scale of synchronization may be important. The multiscale predictability analysis presented here can aid in analyzing these scale-differential effects.

In conclusion, the MS-ARX method is well-adapted to the high-level analysis of neural signals from different brain nuclei at multiple scales.

Acknowledgments

This work was supported in part by the NIH NINDS RO1NS42402, HRSA DIBTH0632, PA Tobacco Settlement Funds Biomedical Research Grant, PSUHMC Movement Disorders Brain Repair Fund, NCCAMR21 AT001607 to T. Subramanian, and GRSA to T. Gilmour via PA Tobacco Settlement Funds (the Pennsylvania Department of Health specifically disclaims responsibility for any analyses, interpretations, or conclusions).

References

M. N. Nounou and H. N. Nounou, “Improving the prediction and parsimony of ARX models using multiscale estimation,” Applied Soft Computing Journal, vol. 7, no. 3, pp. 711–721, 2007.

A. Nambu, H. Tokuno, and M. Takada, “Functional significance of the cortico-subthalamo-pallidal 'hyperdirect' pathway,” Neuroscience Research, vol. 43, no. 2, pp. 111–117, 2002.

G. Heimer, I. Bar-Gad, J. A. Goldberg, and H. Bergman, “Dopamine replacement therapy reverses abnormal synchronization of pallidal neurons in the 1-methyl-4-phenyll-1,2,3,6-tetrahydropyridine primate model of parkinsonism,” Journal of Neuroscience, vol. 22, no. 18, pp. 7850–7855, 2002.

A. Sharott, P. J. Magill, D. Harnack, A. Kupsch, W. Meissner, and P. Brown, “Dopamine depletion increases the power and coherence of β-oscillations in the cerebral cortex and subthalamic nucleus of the awake rat,” European Journal of Neuroscience, vol. 21, no. 5, pp. 1413–1422, 2005.

N. Mallet, A. Pogosyan, A. Sharott et al., “Disrupted dopamine transmission and the emergence of exaggerated beta oscillations in subthalamic nucleus and cerebral cortex,” Journal of Neuroscience, vol. 28, no. 18, pp. 4795–4806, 2008.

U. Ungerstedt and G. W. Arbuthnott, “Quantitative recording of rotational behavior in rats after 6-hydroxy-dopamine lesions of the nigrostriatal dopamine system,” Brain Research, vol. 24, no. 3, pp. 485–493, 1970.

C. A. Lieu, A. R. Kunselman, B. V. Manyam, K. Venkiteswaran, and T. Subramanian, “A water extract of Mucuna pruriens provides long-term amelioration of parkinsonism with reduced risk for dyskinesias,” Parkinsonism and Related Disorders, vol. 16, no. 7, pp. 458–465, 2010.

T. P. Gilmour, B. Piallat, C. A. Lieu et al., “The effect of striatal dopaminergic grafts on the neuronal activity in the substantia nigra pars reticulata and subthalamic nucleus in hemiparkinsonian rats,” Brain, vol. 134, no. 11, pp. 3276–3289, 2011.

E. A. Clement, A. Richard, M. Thwaites, J. Ailon, S. Peters, and C. T. Dickson, “Cyclic and sleep-like spontaneous alternations of brain state under urethane anaesthesia,” PLoS One, vol. 3, no. 4, Article ID e2004, 2008.

M. Cassidy, P. Mazzone, A. Oliviero et al., “Movement-related changes in synchronization in the human basal ganglia,” Brain, vol. 125, no. 6, pp. 1235–1246, 2002.

D. Williams, M. Tijssen, G. Van Bruggen et al., “Dopamine-dependent changes in the functional connectivity between basal ganglia and cerebral cortex in humans,” Brain, vol. 125, no. 7, pp. 1558–1569, 2002.

A. A. Kühn, A. Kupsch, G. H. Schneider, and P. Brown, “Reduction in subthalamic 8–35 Hz oscillatory activity correlates with clinical improvement in Parkinson's disease,” European Journal of Neuroscience, vol. 23, no. 7, pp. 1956–1960, 2006.

A. A. Kühn, A. Tsui, T. Aziz et al., “Pathological synchronisation in the subthalamic nucleus of patients with Parkinson's disease relates to both bradykinesia and rigidity,” Experimental Neurology, vol. 215, no. 2, pp. 380–387, 2009.

A. A. Kühn, D. Williams, A. Kupsch et al., “Event-related beta desynchronization in human subthalamic nucleus correlates with motor performance,” Brain, vol. 127, no. 4, pp. 735–746, 2004.

V. Gradinaru, M. Mogri, K. R. Thompson, J. M. Henderson, and K. Deisseroth, “Optical deconstruction of parkinsonian neural circuitry,” Science, vol. 324, no. 5925, pp. 354–359, 2009.