|Home | About | Journals | Submit | Contact Us | Français|
Nonparametric system modeling constitutes a robust method for the analysis of physiological systems as it can be used to identify nonlinear dynamic input–output relationships and facilitate their description. First- and second-order kernels of hippocampal CA3 pyramidal neurons in an in vitro slice preparation were computed using the Volterra–Wiener approach to investigate system changes associated with epileptogenic low-magnesium/high-potassium (low-Mg2+/high-K+) conditions. The principal dynamic modes (PDMs) of neurons were calculated from the first- and second-order kernel estimates in order to characterize changes in neural coding functionality. From our analysis, an increase in nonlinear properties was observed in kernels under low-Mg2+/high-K+. Furthermore, the PDMs revealed that the sampled hippocampal CA3 neurons were primarily of integrating character and that the contribution of a differentiating mode component was enhanced under low-Mg2+/high-K+.
Communication between neurons and neuronal ensembles in the brain is mediated by excitatory and inhibitory interactions. If these interactions are disrupted or become imbalanced, pathological communications can emerge, manifesting as abnormal electrical events or rhythms.
Epilepsy is a neurophysiological disease in which neuronal communications are spontaneously punctuated by abnormal activity consisting of seizures, interictal bursts, or fast ripples . Although the disease has a complex etiology, many in vitro and animal models suggest epileptogenesis is dependent in part on a shift in the dynamic balance of excitation and inhibition. For example, application of pharmacological agents such as bicuculline or gabazine to a hippocampal slice preparation in combination with elevated extracellular potassium produces hyperexcitable conditions that are conducive to the occurrence of seizure-like events (SLEs) . Bicuculline and gabazine are γ-aminobutyric acid receptor antagonists and therefore block the function of inhibitory synapses, while elevated potassium depolarizes the neuronal membrane. Similarly, if a slice preparation is bathed in low-magnesium solution, spontaneous SLEs are generated due to the hyperexcitable conditions derived from reduced Mg2+ blockade of N-methyl-d-aspartate receptors and Ca2+ channels, as well as charge-screening effects . In rodents, administering kainic acid over time leads to the sprouting of extra recurrent excitatory connections which are hypothesized to increase excitability and rhythmic susceptibility, thereby leading to chronic epileptiform activity .
With animal models of epilepsy such as the kainate model, it is possible to observe long-term morphological changes in brain tissue associated with epileptogenesis, such as changes in synaptic connectivity and neuronal cell loss. With in vitro epilepsy models and during the early phases of animal models, however, plasticity is typically not a contributing factor, and subtle changes at the molecular level prime the neuronal network for the emergence of epileptiform activity . Because these changes are not directly observable, they must be quantified by indirect means.
One approach is to address the problem from a dynamical perspective using nonparamet- ric modeling. The dynamics of a nonlinear system can be mapped by administering bandlimited white noise as a stimulus and estimating the Wiener or Volterra kernels of the system , assuming a stationary and causal process. The kernels encode the transfer characteristics of the system, quantitatively defining the form of the system output corresponding to a given input . The first-order kernel describes the linear response, while higher-order kernels describe the nonlinear response. The Volterra–Wiener approach is a black-box approach in that it does not require specific knowledge of the internal system structure or properties (e.g., ion channel conductances), thus the methodology has been successfully applied to model neuronal systems [7–9] as well as other complex physiological systems such as the heart . Once the kernels are identified, they can be eigen-decomposed into a set of component impulse responses, of which the most significant in terms of their relative contribution to the overall system response are the principal dynamic modes (PDMs) [11, 12]. PDMs with monophasic profiles reflect a predominantly integrating character, and n-phasic PDMs (for n>1) in general reflect a differentiating character. Neuronal PDMs are functionally related to how signals are encoded and decoded by respective network neurons, so an alteration in signaling or membrane properties would be expected to produce changes in neuronal coding and hence changes in the shape or polarity of identified modes. Therefore, characterization of PDMs may serve to identify system changes related to pathological communication and yield a surrogate measure of epileptic susceptibility—which may prove particularly useful in circumstances where no obvious focus or neurophysiological defect can be causally linked to the emergence of epileptiform phenomena.
In this article, we applied the Laguerre expansion technique (LET)  to estimate the kernels of pyramidal neurons recorded from rat hippocampal CA3 in slices prepared with standard artificial cerebrospinal fluid (ACSF), simulating normal conditions, followed by estimation of kernels during perfusion with low-magnesium ACSF, simulating epileptogenic conditions. We then compared the PDMs derived from eigen-decomposition of the kernels to identify features that distinguish the pre- and post-treatment specimens.
Brain slices were prepared from 15- to 21-day-old Wistar rats. Animals were anesthetized with halothane and decapitated in accordance with the guidelines of the Animal Care Committee (University Health Network, Toronto). The brain was quickly removed and placed in ice-cold, continuously oxygenated (95% O2 and 5% CO2) artificial cerebrospinal fluid for 2 to 4 min before the brain slices were prepared. ACSF contained (in millimolar units): 125 NaCl, 26 NaHCO3, 2.5 KCl, 1.8 CaCl2, 2 MgCl2, 1.25 NaH2PO4, and 10 glucose. Hippocampal slices were prepared at room temperature as follows: the brain was hemisected along the midsagittal line, and the cerebellum and the forebrain were removed. The dorsal cortex then was cut parallel to the longitudinal axis, and the brain was fixed, ventral side up, to an aluminum block with cyanoacrylate glue. The aluminum block was secured at a 12° angle in a Vibratome (Series 1000, Technical Products International, St. Louis, MO, USA) with the caudal end of the brain facing the blade. Slices (500 μm thick) were incubated for at least 1 h before being transferred to an interface-type chamber for electrophysiological recording. Standard ACSF was used for incubation and for perfusion of slices in the recording chamber. Once the whole-cell patch was established and stabilized in the current clamp mode, the perfusion solution was changed to low-Mg2+/high-K+ ACSF (125 NaCl, 26 NaHCO3, 5 KCl, 1.8 CaCl2, 0.5 MgCl2, 1.25 NaH2PO4, and 10 glucose). ACSF for perfusion of slices was continuously aerated with carbogen (95% O2 and 5% CO2), and the slice in the recording chamber was also aerated with humidified carbogen. The temperature of the recording chamber was regulated to 33°C. Intracellular whole-cell patch recordings at a sampling rate of 10 kHz were made from CA3 pyramidal neurons using borosilicate glass tubing filled with a solution containing: 135 K-Glu, 10 HEPES, 5 KCl, and 2 Mg–ATP (275 ±5 mOsm, pH 7.2 adjusted with KOH). All chemicals used in the experiments were purchased from Sigma-Aldrich Canada Ltd. (Oakville, ON, Canada).
The input–output dynamics of neuronal systems can be coded by Volterra–Wiener functional expansions. An efficient method for estimating the system kernels is the Laguerre expansion technique, first proposed by Wiener  and developed into a practical tool by Marmarelis . The components of the response can be further characterized by finding a set of dynamic modes relating to the kernels through eigen-decomposition of a properly-constructed matrix of kernel values .
The dynamic response of a causal, stationary nonlinear dynamic system is represented mathematically by the discrete-time Volterra series as follows:
where kp is the p th order kernel and x(n) and y(n) are the system input and output, respectively. With LET, the response is expanded on an orthonormal basis of discrete Laguerre functions, transforming Eq. 1 into the following :
where is the convolution of the j th discrete Laguerre function, , with the input,
and represent the Laguerre expansion coefficients of the respective kernels. The expansion coefficients are estimated by least-squares fitting of the output, y(n), to the given stimulus (e.g., Gaussian white noise). For this paper, kernels up to second order were estimated. The number of Laguerre functions used in the expansion, L, and the decay parameter α (from Eq. 4), were selected through an incremental parameter search (sweeping L from 5 to 30 and α from 0.3 to 0.9) by determining which combination of parameters yielded the lowest normalized mean-squared error  in terms of a comparison between the computed response and the recorded biological response.
The neuronal dynamic modes of the system were computed by eigen-decomposition of the following real symmetric square matrix of kernel values:
whereby the response of the second-order model is coded quadratically as a convolution of the kernel values with the time-delayed inputs:
(superscript T represents the transpose) and
Since Q is a real symmetric square matrix, it can be diagonalized such that
whereby D is a diagonal matrix of eigenvalues of Q, and V is the matrix of eigenvectors. To select the principal dynamic modes of the system, the eigenvalues were assessed according to their significance (in terms of relative contribution to the overall response). For each significant eigenvalue, λs, such that , where 0<χ<<1 is the significance level, the corresponding eigenvector is used to derive the component impulse response defining the sth PDM of the system:
where is the discrete Kronecker delta. The sth component mode output is generated through convolution of the sth PDM with the input x(n). Therefore, the quadratic nonlinear model response produced using S selected PDMs is given by
where the asterisk denotes the discrete-time convolution operation.
System identification was carried out on CA3 pyramidal cells that were patched in the stratum pyramidale, and the input to each cell consisted of Gaussian white noise (GWN) current injected in two separate 1.8-s epochs: first, during perfusion of the slice specimen with normal ACSF and, second, following the switch to perfusion with low-Mg2+/high-K+ ACSF. GWN sequences of zero mean and unit variance were initially generated in MATLAB (The MathWorks, Natick, MA, USA) then scaled so that the peak-to-peak magnitude of the injected current did not exceed 1.0 nA in order for action potentials not to be induced under normal (non-epileptogenic) conditions. The low-Mg2+/high-K+ ACSF induced spontaneous spiking activity which was intermittent (Fig. 1). Under epileptogenic conditions, GWN was injected during the quiescent, non-spiking periods, and epochs that were free of spikes used to obtain the kernel estimates. The response of a given neuron to the GWN epochs was recorded and used to estimate the treatment-dependent Volterra kernels up to second order. Stimulation and recording were performed essentially simultaneously through the same patch electrode by partitioning the duty cycle equally between the two tasks at 10 kHz.
The choice of number of the Laguerre functions to be used in the functional expansion affects the shapes of the kernels and hence the profiles of the neuronal modes that are derived from the kernels. The choice of α (Eq. 4) is likewise important as it determines the rate of exponential asymptotic decay of the Laguerre functions and hence must be large enough to cover the duration (memory) of the response. The chosen parameters were selected by minimizing the prediction normalized mean squared error (NMSE) for the sample population of neurons tested (n=5). The parameter values of α=0.80 and L=20 yielded the best NMSE on average, with a mean error of 10.9 ±3.5% for normal conditions and 13.6 ±6.0% for epileptiform conditions. The output of the kernel model for a sample neuron, as computed from the kernel estimates using Eq. 6, is shown against the actual recorded neuronal output in Fig. 2. From the time series, the contrast between the normal and epileptiform activity is distinct; the fluctuations on average appear smaller and more randomly distributed in the normal case, whereas in the epileptogenic case, they are larger and have greater consistency in terms of frequency and amplitude characteristics, appearing more oscillatory.
Estimated first-order kernel traces corresponding to the same neuron for the two different treatment regimes are shown in Fig. 3a. The defining characteristics of the first-order kernel across the sample population of neurons seem to be comparable in each of the treatment regimes: namely, the first-order kernels are characterized by an initial steep rise and sharp peak, followed by a gradual decay tail superimposed with low-amplitude ripples. Nonetheless, several differences are apparent. First, the amplitude and prevalence of kernel tail ripples under epileptogenic conditions are more pronounced. The cross-section of the kernel peak is also perceptibly wider. Typical second-order kernels for normal and epileptiform-conditioned neurons are shown in Fig. 3b, c. The second-order kernels of CA3 neurons under normal conditions did not show any prominent features but displayed mostly low-amplitude complex or noise-associated ripples, whereas the second-order kernels of epileptiform neurons featured prominent ridges along the diagonal rising above the mean level of fluctuations.
The principal dynamic modes were calculated from the kernel estimates, as outlined in Section 2. Because of the quadratic form of the model (Eq. 10), two PDMs were found to be necessary to reconstruct the system output reasonably, especially for μs,0 ~0. The output prediction using a single PDM was generally not acceptable because of the tendency of quadratic rectification to produce a single polarity response when both positive and negative polarities were featured in the original output signal (Fig. 2), therefore yielding a high prediction error (Table 1). By contrast, two PDMs produced adequate output predictions, even when the PDMs were mirror images of one another (which indicates one true PDM for the system). As seen in Table 1, the two-PDM response improved the prediction NMSE by 4.28 ±1.67% for the normal treatment and 3.07 ±2.59% for the epileptogenic treatment over that of the conventional kernel model, as given by Eq. 1. This is presumably because the kernels reflect, to some degree, contributions from noise or artifacts that may be weeded out in the PDM model by selecting only the predominant neuronal modes associated with the true response and omitting those modes whose eigenvalues are insignificant and are likely to be associated with noise or artifacts.
The significance of each PDM was defined by the relative magnitude of its eigenvalue to the sum of all the eigenvalue magnitudes. The criterion as to how to set the level of significance is somewhat arbitrary, though this study has circumvented the issue by not defining a fixed threshold but rather by examining the relative significance of each neuronal mode (see Fig. 4). Overall, the contribution of the first PDM to the system response, as quantified by the magnitude of its eigenvalue, was 55.23 ± 2.63% for the normal treatments (n=5) and 51.29 ±6.35% for the epileptogenic treatments (n=5). Likewise, the significance of the second PDM was 37.81 ±2.94% and 28.22 ±4.37%, respectively. The combined contribution of the first and the second PDMs to the measured system response was 90.0% for normal and 79.5% for the epileptogenic treatments. By comparison, the third PDM contributed only 3.10 ±0.01% to the overall response under epileptogenic conditions but contributed even less at 0.96 ±0.19% for the normal treatment case. Interpreted another way, there was a more than threefold relative increase in the significance of the third PDM between normal and epileptogenic treatments, whereas there was only a modest decrease in the significance of the first two PDMs.
The coding functionality of the various PDMs within the context of system dynamics was not directly evident from visual inspection of the PDMs alone. The input–output mapping of the PDMs was better characterized by convolving them with a rectangular pulse to determine their transformative effect on the input. The first two PDMs of both the normal and the epileptogenic treatments exhibited predominantly integrating characteristics since the convolution produced unidirectional amplitude changes that followed the form of the stimulus (Figs. 5d, e and and6d,6d, e). This was in spite of the somewhat complex ripples or undulations that were encountered in several of the integrating PDMs and which were suggestive of more intricate (possibly differentiating) functionality. When the sample neurons were subjected to epileptogenic treatment, however, convolution of the third PDM with the rectangular pulse yielded component responses that were biphasic or triphasic and which demarcated the rising and/or descending edges of the pulse, respectively (Fig. 6f), indicating differentiating character. During normal conditions, the third PDM was predominantly integrating (Fig. 5f), although slight deflections could be seen at the pulse beginning or end in some neurons. This implies that differentiating character may be masked under normal conditions (at least for the sample population of cells recorded from the CA3 stratum pyramidale) but could be promoted or activated by certain factors or circumstances, including hyperexcitable conditions—though not all factors need to be pathological, as such dynamic functionality may serve other neural coding purposes.
The investigation of changes in nonlinear system dynamic properties associated with epileptogenic conditions as compared to normal or non-epileptic conditions was accomplished by nonparametric modeling, which relies on a black-box approach whereby the properties of the system are measured and characterized by examining the dynamic relationship between the system inputs and outputs. Knowledge of system structure or internal mechanisms generating the observed dynamics is not required , which is distinct from parametric models such as those based upon differential equations. Nonparametric kernel models are suitable for describing coding transformations and are therefore sensitive to changes in system properties linked to the emergence of pathological conditions. Associating kernel features with specific physiological processes or structures, however, is not a transparent task, and hence the tandem use of parametric models that emulate the biological system may provide the means for deducing such associations .
Our present work has revealed that notable changes in the system kernels and PDMs of CA3 pyramidal neurons are induced upon switching to an epileptogenic preparation consisting of low-Mg2+/high-K+ ACSF. These differences include the shape, amplitude (significance), and functional character of system kernels and modes (Figs. 4, ,5,5, and and6).6). The first-order kernels were found to be somewhat broader under low-Mg2+/high-K+. Bardakjian et al.  previously reported that the duration of the first-order kernel in a Hodgkin–Huxley model of two bidirectionally coupled neurons was affected by varying the synaptic strength. Consequently, we can speculate that the kernel broadening may be attributed to changes in synaptic function induced by low-Mg2+/high-K+ conditions. The kernel features and response to GWN current injection of the neurons when subjected to epileptogenic conditions were also more oscillatory (Fig. 3).
The second-order kernels of normal CA3 neurons did not present with any notable features, other than low-amplitude ripples that may be attributed to noise. However, the second-order kernel of epileptiform neurons showed prominent ridges along each side of the diagonal that were significantly larger than the average baseline kernel fluctuations. These ridges are similar to those reported by Naylor , who conducted a study of changes in Wiener kernels in a kindling model of epilepsy. In another study by Bardakjian et al. , the boost in the second-order kernel profile was not the result of epileptogenic conditions but was more a reflection of the switch in GWN stimulation from subthreshold to suprathreshold intensity and therefore the invocation of a nonlinear component response, whereas subthreshold excitation only elicited the linear portion of the response and hence did not yield a significant second-order kernel. In this respect, the observed changes in the second-order kernel under low-Mg2+/high-K+ epileptogenic conditions may be due to a systemic lowering of the relative excitation threshold as a consequence of the induced hyperexcitability in the system. In general, however, the changes associated with pre- and post-epileptogenic treatments in system properties and coding are subtle and are not reflected by any gross morphological changes in tissue properties or synaptic connections, therefore making diagnostic detection of such pathological transitions difficult without applying specialized techniques such as nonlinear system identification.
Neuronal modes are useful because they serve as a basis over which a complex dynamic system response can be expanded and therefore broken into simpler dynamic components that are more easily interpreted. Marmarelis and Orme  have categorized two principal dynamic mode classes in neural systems, namely integration (coding for amplitude/displacement) and differentiation (coding for rates of change). Depending on the manner in which modes are combined, or mixed, and on their relative weighting or significance, the system may be interpreted as being more amplitude sensitive or rate sensitive. Integrating and differentiating functionality in neuronal systems has been identified as being important to network signaling and communication. Vigmond et al.  showed that the electric-field coupling generates waveforms that are similar to the first and second derivatives of the sources, which suggests the involvement of differentiating modalities in cellular communication. Moreover, research has revealed the existence of neurons with receptive fields tuned to stimuli or cues relating to displacement, velocity, acceleration, or a mixture of these attributes. Examples of such neuronal classes can be found in the vestibular system of primates , the whisker-trigeminal system in rodents , the motor cortex of felines , and the crab chordotonal organs , to name a few. These studies demonstrate the significance and universality of integrating and differentiating modalities in neuronal information coding, from the cellular to the system level.
According to our analysis, the prevalent dynamics of sampled hippocampal CA3 pyramidal cells were attributable to a small number of PDMs (two to three), as demonstrated by the successful kernel and PDM model predictions (Fig. 2). The first two dominant PDMs were mainly of integrating character under both normal and epileptogenic conditions (Figs. 5 and and6).6). Since the relative significance of the two dominant PDMs was approximately 90% and 80%, respectively (Fig. 4), the sampled hippocampal CA3 neurons were strongly amplitude-sensitive, but they transitioned to being more rate-sensitive after exposure to hyperexcitable conditions, as indicated by the increase in the significance and differentiating character of the third PDM. As the hippocampus plays important roles in memory formation and spatial navigation , we infer that integrating functionality may be essential for sequential learning and the persistence of episodic memory traces in the highly auto-associative networks of the CA3. This may help to explain why the CA3 pyramidal neurons tested displayed a primarily integrating character. On the other hand, the hippocampus becomes sensitized by high-frequency stimulation, as in the case of long-term potentiation (LTP) , which implies the action of a rate-dependent (differentiating) process. Abnormal high-frequency activity is also associated with the epileptiform condition [25, 26]. In the kindling model of epilepsy, the system threshold for initiating a seizure-like event becomes lower as the number of SLEs is increased , until kindling is no longer required for SLEs to occur. In other words, the system learns to seize on its own after a finite number of induced SLEs. Altogether, this suggests the epileptic brain may become progressively more susceptible to seizures partly due to its heightened rate sensitivity in combination with the mechanisms of LTP.
In conclusion, nonlinear system dynamic properties—and hence the way in which the system handles and processes information through coding transformations—are sensitive to pathological conditions such as those produced in the low-Mg2+/high-K+ model of epilepsy. These changes can be quantified through nonparametric techniques such as estimation of the system kernels and dynamic modes. Neuronal modes code information related to system functionality, and changes to their profile and/or polarity brought about by abnormal conditions may provide relevant markers for identifying or diagnosing dynamic diseases such as epilepsy.