Nonlinear dynamical model based control of in vitro hippocampal output

Introduction

A neural prosthesis is a prosthetic device that interfaces with the nervous system to improve or restore impaired neural function (Berger et al., 1994; Schwartz, 2004; Patil and Turner, 2008). The neuroprosthetic technology has been advancing rapidly (Bernotas et al., 1986; Creasey et al., 2004; Mayberg et al., 2005; Hochberg et al., 2006; Allison et al., 2007; Stacey and Litt, 2008). Neural prostheses can be categorized according to the directions of the signal communication between the device and the nervous system (Turner et al., 2005; Song et al., 2007). The first category of neural prostheses attempts to decode neural signals and then to activate an external object. An example would be the neuroprobes decoding motor cortex signals to control a robotic arm (Donoghue, 2002; Nicolelis, 2003; Taylor et al., 2003). The second kind of neural prostheses encodes external sensory stimuli and intends to activate the nervous system. Examples are cochlear implants and artificial retinas (Middlebrooks et al., 2005; Weiland et al., 2005). The third kind of neural prostheses, which forms a bi-directional closed-loop system with the nervous system, receives incoming neural signals from one nervous region and sends its output to activate another nervous system region (Berger et al., 2001, 2011). For the neural prosthesis that involves stimulation to the nervous system, the output system responses could be influenced by the stimulation parameters such as location, intensity, and frequency. Because the signal transformation in the nervous system is nonlinear, it is also important to consider the nonlinearity between stimulation patterns and the output responses. Without considering this nonlinear relationship, large deviations between the device-evoked responses and the desired responses are expected. In practice, such deviations can be mitigated by tuning the stimulation parameters (Lauer et al., 2000; O'Suilleabhain et al., 2003; McIntyre et al., 2004; Tellez-Zenteno et al., 2006; Rupp and Gerner, 2007; Albert et al., 2009; McLachlan et al., 2010). This optimization procedure is typically performed manually and empirically, e.g., assuming a static and linear relation between the stimulation pattern and the desired responses, and then searching for the optimal ratio between the stimulation intensity and the outcome responses via a trial-and-error procedure. To formally solve this important problem, one needs to develop a rigorous stimulation paradigm that takes the (nonlinear dynamical) relationship between stimulation signals and system responses into account (Liu and Oweiss, 2010; Liu et al., 2011).

We are in the process of developing a neural prosthesis to restore the long-term memory formation function of the hippocampus that is lost in Alzheimer's disease, stroke, epilepsy, or other neurological disorders. Our concept of such a prosthetic device is a biomimetic model of the input–output nonlinear dynamics of the hippocampus—a model that captures how hippocampal circuitry re-encodes, or transforms, incoming spatio-temporal patterns of neural activity (i.e., short-term memories) into outgoing spatio-temporal patterns of neural activity (i.e., long-term memories) (Squire, 1992; Berger et al., 2001, 2005; Burgess et al., 2002). We have shown in rodents, both in vitro (Chan et al., 2004; Hsiao et al., 2006) and in vivo (Song et al., 2007, 2009; Berger et al., 2011, 2012), that a nonlinear hippocampal model is capable of predicting accurately the output signals based on the ongoing input signals in the hippocampus. In this study, we extend this concept by developing a rigorous stimulation paradigm with control theory, and then implementing it rat hippocampal slices.

The intrinsic circuitry of the hippocampus consists of three major subregions: dentate gyrus (DG), CA3, and CA1 as shown in Figure 1A. This trisynaptic circuit can be maintained in a transverse slice preparation (Andersen et al., 1969, 2000; Amaral and Witter, 1989). The signal transformations in all three regions are highly nonlinear and dynamical (Berger et al., 1988; Sclabassi et al., 1988; Bartesaghi et al., 2006). From an engineering perspective, the hippocampal circuit can be viewed as a cascade of input–output transfer functions between the DG, CA3, and CA1 subregions (Figure 1B). In the context of extracellular recording as in this study, the evoked field potentials in each subsystem are measured as input–output signals. For example, the CA3 response (field excitatory postsynaptic potentials amplitude, fEPSP) can be used as the input signal to CA1, and the CA1 response can be considered as the final system output. A schematic diagram of such a hippocampal prosthesis is shown in Figure 1C, where CA3 is damaged, thus the signal transmission from DG to CA1 cannot be completed. In the replacement scenario, the prosthesis model processes the DG signals and generates optimal stimulations to elicited desired output response in the CA1 region.

FIGURE 1

Figure 1. (A) A rat hippocampal slice and its major intrinsic pathways. The input signals from perforant path fibers excite dentate granule cells. Dentate output, in turn, excites CA3 pyramidal cells through mossy fibers. Output from CA3 is transmitted to CA1 pyramidal cells through Schaffer collaterals. This so-called “trisynaptic pathway” is the principal network involved in hippocampal neuronal information processing. (B) A block diagram showing the trisynaptic pathway in a hippocampal slice. (C) A schematic diagram of a hippocampal prosthesis model functionally replacing the original pathway, where CA3 is damaged, so the signal transmission cannot be completed. This bi-directional prosthetic device receives incoming neural signals from one hippocampal region (DG) and sends its output to stimulate another hippocampal region (CA1). (D) The proposed modeling-control paradigm to optimize the stimulation patterns. In this framework, the desired CA1 output is first predicted with the DG signal by the trajectory model, and then converted to the desired stimulation patterns through the inverse model. The desired stimulation patterns then drives the output system (CA1) to the desired output responses.

The successful implementation of such a device depends on three sequential components. First, the device must capture incoming neuronal signals reliably from the input region. Second, it must mimic the damaged system precisely through a computational model. Finally, the device should reproduce the desired responses in the output region through electrical stimulation. Thus, through bi-directional communication with the brain, the prosthetic device could essentially bypass the damaged region and substitute the lost function.

This paper describes the procedure of deriving optimal stimulation patterns using an inverse control concept (Houk, 1988; Widrow and Walach, 1996; Camacho and Bordons, 2003; Normann, 2007). The “trajectory model” is a model that predicts the desired output response based on the input patterns. This model can be developed using available knowledge or built directly from experimental input–output data. The stimulation-response properties of the output system is described as the “plant model.” The “inverse plant model” describes a system whose transfer function is the inverse transformation of the plant model (Widrow et al., 1993; Widrow and Plett, 1997; Karniel et al., 2001). This inverse transformation can be determined once the input–output transformation of the plant model is fully explored. Once these three models (i.e., trajectory, plant, and inverse plant models) are built, the signals flow like what is shown in Figure 1D. Signals recorded from the DG (input) system pass through the trajectory model to predict the desired output. The inverse plant model is then used to derive the desired stimulation amplitudes from desired output. Finally, the CA1 (output) region generates the controlled output responses. Results show that the strategy described in this paper is able to control CA1 output activities (shown in Figure 1D as “Controlled CA1 responses”) to replicate the CA1 activities recorded from the hippocampal slice with intact trisynaptic pathway (shown in Figure 1B as “Trisynaptic CA1 responses”).

Materials and Methods

The proposed modeling-control paradigm was verified using an in vitro rat hippocampal slice preparation. Section “Experimental Procedures” provides an explanation of the methodology used to prepare the hippocampal slices, and the description of our electrophysiology experimental setup. Section “Modeling-Control Paradigm Implementation and Data Collection” describes the estimations and validations of the trajectory model, the plant model and the inverse plant model, and the associated data collection and analysis procedures. The overall experimental protocol is described in section “Modeling-Control Framework Experiment Protocol.”

Experimental Procedures

Acute hippocampal slice preparation

Hippocampal slices from 8 to 10-week-old male Sprague-Dawley rats (250–300 gm) were prepared. The animals were first anesthetized with halothane (Halocarbon Laboratory, USA) and then decapitated. Their skulls were rapidly removed and the brain was carefully extracted. Hippocampi were separated from the cortices in an iced sucrose buffer solution (Sucrose 206 mM; KCl 2.8 mM; NaH₂PO₄ 1.25 mM; NaHCO₃ 26 mM; Glucose 10 mM; MgSO₄ 2 mM; Ascorbic Acid 2 mM). Hippocampal slices 400 micrometers thick were sliced transversely from the ventral hippocampi using a vibratome (Leica VT1000S, Germany). The slices were incubated for at least 1 h in 2 mM MgSO₄ artificial cerebral spinal fluid (aCSF) at room temperature, to equilibrate. During each electrophysiological recording session, one slice at a time was transferred to the planar multielectrode array. The array attached with a circular plastic chamber and perfused with normal aCSF (NaCl 128 mM; KCl 2.5 mM; NaH₂PO₄ 1.25 mM; NaHCO₃ 26 mM; Glucose 10 mM; MgSO₄ 1 mM; Ascorbic Acid 2 mM; CaCl₂ 2 mM) maintained at room temperature (24~26°C). In the recording chamber, each slice was held down by a metallic ring with nylon mesh attached to it. The positioning of the slice was accomplished by manipulating the ring with a small brush. All the solutions were bubbled with 95% O₂ and 5% CO₂ mixed gas. The protocol described above was approved by the Department of Animal Resources and Institutional Animal Care and Use Committee at the University of Southern California.

Electrophysiological recording setup and procedures

Electrophysiology data were collected through an extracellular recording technique using an MEA60 system (Multi Channel Systems, Germany), as seen in Figure 2. This system consisted of pre-amplifiers (1200× gain), a data acquisition device (MC_Card), and an 8-channel stimulus generator (STG1008), all operated using software provided by Multi Channel Systems (MC_Rack V3.2.0 and MC_Stimulus V2.0.6). A conformal 60-channel planar multielectrode array was made specifically for this study. The geometry of this conformal array was designed to match the cytoarchitecture of the hippocampus slices (Figure 2C) and was platinum based. Details in fabrication and the arrangement of the array can be found in Gholmieh et al. (2006) and Taketani and Baurdy (2006). Collected data were sampled at a frequency of 10 kHz per channel and were recorded using MC_Rack. The MEA60 system was assembled over an inverted microscope (Leica DM-IRB, Germany). In each experiment, the position of the slice on the MEA was captured by a digital image capture system (Diagnostic Instruments, Spot RT Digital Camera, USA) with SPOT (V4.6.4.3) software and Adobe Photoshop (Adobe V7.0, USA).

FIGURE 2

Figure 2. A photo of the electrophysiological recording system. (A) The MEA60 system and (B) the conformal planar MEA (Gholmieh et al., 2006; Taketani and Baurdy, 2006). (C) A photomicrograph of a hippocampal slice on the conformal MEA. The set alignment of this array is according to rat hippocampal cytoarchitecture covering major subregions of DG, CA3, and CA1. The waveforms represent the trisynaptic response of the hippocampal slice recorded in each region. The white lines indicate the amplitude measurement of DG population spike amplitude and CA1 fEPSP amplitude (see section “FARIT-Induced Trisynaptic Data Collection and Analysis”).

Stimulation and data collection procedures

In this study, biphasic currents with a 100 μs duration in each phase were applied to all stimulation patterns. Different stimulation trains were programmed in MC_Stimulus and used to study the nonlinear properties of different regions. There was a 5–7 min waiting period between each stimulus train. The evoked neural responses were simultaneously recorded from different regions. The channels were first selected based on the placement of the recording electrodes on the cytoarchitecture of the slices (i.e., DG channels must be in the DG region, CA3 channels must be in the CA3 region, etc). Among those channels, the channels with the largest response (dendritic population spike or EPSP) amplitudes are further selected and analyzed. The main purpose of this procedure is to find the most representative responses for each region, channels with small responses or inappropriate placements are not analyzed because their recordings may reflect a non-cell-body placement or a mixture of activities from multiple regions.

Modeling-Control Paradigm Implementation and Data Collection

DG-CA1 trajectory model implementation

FARIT-induced trisynaptic data collection and analysis. An external bipolar electrode of twisted Nichrome wires was used to elicit the trisynaptic response. Paired-pulse or quadruplet-pulse electrical stimulation was applied to the perforant pathway of each slice using the external electrode to generate electrophysiological responses throughout the trisynaptic pathway (evoked field potentials in DG, CA3, and CA1, as seen in Figure 2C). When the full trisynaptic response was observed, we stimulated the slice with a series of fixed-amplitude, random inter-impulse-interval trains (FARITs). Four 300-pulse Poisson distributed FARITs of a fixed current intensity (biphasic, 150–300 μA) were delivered to the perforant path (1200 impulses; range of intervals: 2 ms to 5 s; mean frequency: 2 Hz). Response amplitudes from selected channels in DG and CA1 regions were analyzed to build the DG-CA1 trajectory model. The neuron response measurement in DG was population spike amplitude, the amplitude was calculated by averaging the distance between the negative peak and the first positive peak (measure “a” in Figure 2C) and the distance between the negative peak and the second positive peak (measure “b” in Figure 2C) (Houk, 1988). To measure responses in the CA1 regions, the field potential amplitude was defined as the negative peak of the waveform (measure “c” in Figure 2C).

Trajectory model configuration. A single-input, single-output discrete model was derived from Volterra series as expressed below (Marmarelis and Orme, 1993; Marmarelis, 2004):

$\begin{array}{cc}\begin{array}{l}y(n)=k_0+{\displaystyle \sum _{m=0}^M{k}_1}(m)x(n-m)\\ +{\displaystyle \sum _{m_1=0}^M{\displaystyle \sum _{m_2=0}^M{k}_2}}(m_1,m_2)x(n-m_1)x(n-m_2)+\dots \end{array}& (1)\end{array}$

The zeroth order kernel k₀ is the value of output y when the input is absent. First order kernels k₁ describe the relationship between each single input x(n − m) and output y. Second order kernels k₂ describe the relationship between the output y and each unique pair of input x(n − m₁), x(n − m₂). The term n represents time of occurrence of the present impulse in the input–output sequence and m represents the interval of the impulses prior to the present impulse within the kernel memory window M, m = 0 denotes the present input. The input to the system can be expressed as a series of variable-amplitude, random-interval delta functions:

$\begin{array}{cc}x\left(t_i\right)={\displaystyle \sum _{i=1}^I{A}_i}\delta \left(t-t_i\right)& (2)\end{array}$

where i is the index number of impulses and I is the total number of impulses. The time of occurrence of the i^th impulse is t_i. In the DG-CA1 trajectory model experiment, DG population spike amplitude were used as input (A_i) and CA1 fEPSP amplitude were used as output y(n). Because the input amplitude is varied, in order to isolate influence from present input, we considered the zero-lag terms in the original Volterra series (1) independently, as follows:

$\begin{array}{l}y(n)=k_0+k_1(0)x(n)+{\displaystyle \sum _{m=1}^M{k}_1}(m)x(n-m)\\ +k_2(0,0)x(n)x(n)+{\displaystyle \sum _{m_1=1}^M{k}_2}(m_1,0)x(n-m_1)x(n)\\ +{\displaystyle \sum _{m_2=1}^M{k}_2}(0,m_2)x(n)x(n-m_2)\\ +{\displaystyle \sum _{m_1=1}^M{\displaystyle \sum _{m_2=1}^M{k}_2}}(m_1,m_2)x(n-m_1)x(n-m_2)+\dots \end{array}$

and can be then rearranged as:

$\begin{array}{cc}\begin{array}{l}y(n)=k_0+k_1(0)x(n)+k_2(0,0)x{\left(n\right)}^2+{\displaystyle \sum _{m=1}^M{k}_1}(m)x(n-m)\\ +{\displaystyle \sum _{m_1=1}^M{\displaystyle \sum _{m_2=1}^M{k}_2}}(m_1,m_2)x(n-m_1)x(n-m_2)\\ +2{\displaystyle \sum _{m=1}^M{k}_2}(m)x(n)x(n-m)+\dots \end{array}& (3)\end{array}$

The first three terms on right represent the static input–output curve. The last three terms describe the nonlinear dynamical effect of the inputs on the output. In order to reduce the number of open parameters, an estimation of the kernels is facilitated by expanding them on the orthonormal Laguerre basis functions L (Marmarelis, 1993):

$L_l(m)={\alpha }^{(m-l)/2}{(1-\alpha )}^{1/2}{\displaystyle \sum _{k=0}^l{(-1)}^k}\left(\begin{array}{c}m\\ k\end{array}\right)\left(\begin{array}{c}l\\ k\end{array}\right){\alpha }^{l-k}{\left(1-\alpha \right)}^k$

where α is the Laguerre parameter (0 < α < 1) and affects the time extent of the basis functions. The convolution of Laguerre basis functions L and inputs x can be represented as

$v_l(t_i)={\displaystyle \sum _{t_i-\mu <t_j\le t_i}{A}_j{L}_l(t_i-t_j)}$

where A_j is the input spike amplitude in (2), t_i is the time of occurrence of the current impulse in the input–output sequence and t_j is the time of occurrence of the j^th impulse prior to the present impulse within the kernel memory window μ. The adapted Laguerre expansion of Volterra kernels with L basis functions can be rewritten as:

$\begin{array}{cc}\begin{array}{l}y(t_i)=c_0+A_i{c}_1(0)+A_i^2{c}_2(0,0)+{\displaystyle \sum _{l=1}^L{c}_1}(l)v_l(t_i)\\ +{\displaystyle \sum _{l_1=1}^L{\displaystyle \sum _{l_2=1}^L{c}_2}}(l_1,l_2)v_{l_1}(t_i)v_{l_2}(t_i)\\ +2A_i{\displaystyle \sum _{l=1}^L{c}_2}(l)v_l(t_i)+\dots \end{array}& (4)\end{array}$

where c₀, c₁, c₂,… are the kernel expansion coefficients. Since the number of basis functions can be made much smaller than the memory length, the number of open parameters is greatly reduced by this expansion technique. The kernel expansion coefficients (c₀, c₁, c₂,…) can be estimated via the least-squares method, and can be used to reconstruct the Volterra kernels (k_i) using Laguerre basis functions (L)

$k_0=c_0,k_1={\displaystyle \sum _{l=1}^L{c}_l}{L}_l,k_2={\displaystyle \sum _{l_i=1}^L{\displaystyle \sum _{l_i=1}^L{c}_{l_1,l_2}}}{L}_{l_1,l_2}$

CA1 plant model implementation

RARIT-induced monosynaptic data collection and analysis. After collecting the FARIT data for building DG-CA1 trajectory model, in the same slice, paired-pulse stimulation was applied to the stratum radiatum from a pair of stimulation electrodes in the conformal array in order to elicit the monosynaptic CA1 response. The pair of stimulation electrodes was selected according to their location and their ability to evoke typical paired-pulse facilitation. In this set of experiments, the amplitudes of the FARITs were modified to formulate a random-amplitude, random-interval trains (RARITs, Gaussian distributed, mean amplitude: 150 μA, which is the mean CA1 evoked postsynaptic potential amplitude observed in the FARIT-induced trisynaptic dataset). Once the pair of stimulation electrodes were determined, four 300-pulse RARITs were delivered to the slice. A channel from the CA1 region was selected and fEPSP amplitudes were analyzed for suitability in constructing the CA1 plant model.

CA1 plant model configuration. The same Laguerre-Volterra (LV) modeling approach described in section “Trajectory Model Configuration” was applied to build a CA1 plant model. In this set of experiments, the amplitudes of the RARITs (from previous section) were used as measures of the input signal A_i in (2), and the fEPSP amplitudes of CA1 were used as measures of output. Once the input–output transformation of the CA1 system was fully explored, the inverse model can be further derived.

Inverse CA1 plant model implementation

Inverse CA1 plant model configuration. The inverse model was built to transform the output (i.e., desired output of a CA1 region) to the input (i.e., desired input stimulation to a CA1 region). To develop the inverse model based on the LV model, the original Equation (4) was rearranged to:

$\begin{array}{cc}\begin{array}{l}[c_2(0,0)]A_i^2+\left[c_1(0)+2{\displaystyle \sum _{l=1}^L{c}_2(l)v_l(t_i)}\right]A_i\hfill \\ \begin{array}{l} +[c_0+{\displaystyle \sum _{l=1}^L{c}_1(l)v_l(t_i)}\\ +{\displaystyle \sum _{l_1=1}^L{\displaystyle \sum _{l_2=1}^L{c}_2(l_1,l_2)v_{l_1}(t_i)\times v_{l_2}(t_i)-y(t_i)}}]=0\end{array}\hfill \\ \hfill \end{array}& (5)\end{array}$

In (5), the desired output y and the coefficients c₀, c₁, c₂,… were obtained during process of model estimation. All the convolution terms could also be determined using the coefficients and previous stimulation amplitudes A_{i − 1}, A_{i − 2},…. Once all the terms are determined, (5) became a quadratic equation with unknown desired input stimulation (A). It can be simplified as:

$\begin{array}{cc}a{A}^2+bA+c=0& (6)\end{array}$

where

$\begin{array}{l}a=c_2(0,0),b=c_1(0)+2{\displaystyle \sum _{l=1}^L{c}_2}(l)v_l(t_i),\\ c=c_0+{\displaystyle \sum _{l=1}^L{c}_1}(l)v_l(t_i)+{\displaystyle \sum _{l_1=1}^L{\displaystyle \sum _{l_2=1}^L{c}_2}}(l_1,l_2)v_{l_1}(t_i)v_{l_2}(t_i)-y(t_i)\end{array}$

such that the transformation of the inverse model (output to input) became an operation of solving A in (6). In this study, the roots of the quadratic equation can all be calculated from

$\begin{array}{cc}A=\frac{-b+\sqrt{b^2-4ac}}{2a}& (7)\end{array}$

The validation of this inverse model implementation is shown in the result (section “Inverse CA1 Plant Model Implementation and Validation Results”). All the calculated stimulation amplitude were used to recompose to the new stimulation trains, called desired-amplitude RITs (DARITs) as described below.

DARIT-induced monosynaptic data collection and analysis. In this set of experiments, the amplitude of the RARITs (as used in section “RARIT-Induced Monosynaptic Data Collection and Analysis”) were reformed using the optimal stimulation amplitudes calculated from the inverse CA1 plant model (from previous section “Inverse CA1 Plant Model Configuration”), named DARITs. Four 300-pulse DARITs were delivered to the slice through the same pair of stimulation electrodes as RARIT experiments. A channel from the CA1 region was selected and the fEPSP amplitudes were analyzed.

Model validation

In this study, data were evaluated using the Variance Accounted For (VAF) and the Normalized Mean Square Error (NMSE) as described below:

$\begin{array}{l}VAF=(1-\text{var}(Y_i-X_i)/\text{var}(Y_i))\\ NMSE={\displaystyle \sum _i{(Y_i-X_i)}^2}/{\displaystyle \sum _i{Y}^2}\end{array}$

where X is the predicted amplitude of the model, Y is the amplitude analyzed from the recorded data, and var is the variance of the data. Specific data sets were chosen for comparison and are presented in the result section. To evaluate the prediction power of the models, we have used a cross-validation method, i.e., independent datasets are used for model estimation and model prediction. All model goodness-of-fit reported in this paper are obtained using this method.

Modeling-Control Framework Experiment Protocol

The experiment protocol to verify our modeling-control paradigm for an in vitro hippocampal prosthesis involves following steps:

1. Stimulating the perforant path (DG input) with FARITs, and analyzing the trisynaptic responses in DG and CA1; applying DG population spike patterns as input and CA1 fEPSP patterns as output, building a DG-CA1 trajectory model using an LV kernel modeling approach (Figure 3A).

2. Stimulating the Schaffer collaterals (CA1 input) with RARITs, and analyzing the monosynaptic responses in CA1; applying RARITs patterns as input and CA1 fEPSP patterns as output, building a CA1 plant model using the LV kernel modeling approach, and then formulating an inverse CA1 plant model (Figure 3B).

3. Applying DG patterns as the input to DG-CA1 trajectory model to predict desired CA1 output patterns; applying the predicted desired CA1 output as the input to inverse CA1 plant model, to derive the optimal stimulation patterns; stimulating at the Schaffer collaterals with the derived optimal stimulation patterns, and then analyzing the responses in CA1 (Figure 3C).

4. Comparing the CA1 fEPSP amplitudes from Step 1 to those in Step 3.

FIGURE 3

Figure 3. A schematic diagram of the modeling-control experimental protocol. (A) Stimulating the perforant path with FARITs, and analyzing the trisynaptic responses in DG and CA1, to build a DG-CA1 trajectory model. (B) Stimulating the Schaffer collaterals with RARITs, and analyzing the monosynaptic responses in CA1, to build a CA1 plant model. The inverse CA1 plant model can then be formulated. (C) Applying DG patterns (from A) as the input to DG-CA1 trajectory model (built in A) to predict desired CA1 output patterns; applying the predicted desired CA1 output as the input to the inverse CA1 plant model (built in B) to derive the optimal stimulation patterns; stimulating at the Schaffer collaterals with the derived optimal stimulation patterns, and then analyzing the responses in CA1. CA1 fEPSP amplitudes from A and C can then be compared.

Results

The diagrams and performance of the DG-CA1 trajectory and CA1 plant model prediction, and the inverse CA1 plant model implementation are presented in this section. The presented protocol was conducted in six experiments. In each experiment, two sets of data were collected from a hippocampal slice. The first dataset was composed of the FARIT-induced trisynaptic data and was used to build the DG-CA1 trajectory model. The second dataset was composed of the RARIT-induced monosynaptic data and was used to build the CA1 plant model. The resulting two built models and their predictions are presented. This section also include the implementation of the inverse CA1 plant model, and lastly, the validation of the modeling-control paradigm used for regulating CA1 nonlinear dynamics.

DG-CA1 Trajectory Model and the Prediction Results

The FARIT-induced hippocampal trisynaptic data were analyzed for use in building the DG-CA1 trajectory model. The amplitudes of evoked DG population spikes were used as measures of the input to the system, and the amplitudes of evoked CA1 fEPSPs were used as measures of output of the system. An LV kernel model was applied to study the nonlinearity of this system. Examples of the first and the second order LV kernels are shown in Figures 4A,B, respectively.

FIGURE 4

Figure 4. (A) The first order and (B) the second order LV kernel of the DG-CA1 trajectory model. The singular points and lines showing on the edge of each figure (indicated by arrows) represent the effect of present input. (C) A segment of comparison between a FARIT-induced trisynaptic CA1 response amplitude and the amplitude predicted by the DG-CA1 trajectory model. (D) The Q–Q plot of the data distribution between actual trisynaptic CA1 responses recorded from the slice and the outputs predicted by the DG-CA1 trajectory model.

It should be noted that in Figure 4A, the singular point represents the K₁(0) term in (3). The different polarity in Figure 4A manifests the importance of isolating the zero lag terms. In Figure 4B, the singular point represents the K₂(0,0) term in (3), and the two lines indicated by arrows represent the ∑^M_{m1 = 1}∑^M_{m2 = 1}k₂(m₁, m₂) term, while one of the input pairs is the present input (m = 0).

Model estimation was completed using population spike amplitudes and the intervals of the input–output sequences. From all datasets, the slice response amplitudes were analyzed and compared with the predicted amplitudes. The mean VAF was 65.97 ± 17.30%. In Figure 4C, a segment of FARIT-induced trisynaptic CA1 fEPSP amplitudes is compared to its counterpart predicted from the DG-CA1 trajectory model. The result is further confirmed by an overall comparison between the actual responses and model predicted outputs, as shown in Figure 4D. The quantile–quantile (Q–Q) plot demonstrates that the actual trisynaptic CA1 responses recorded from the slice are accurately predicted by the DG-CA1 trajectory model.

CA1 Plant Model and the Prediction Results

The RARIT-induced CA1 monosynaptic data were analyzed for use in building the CA1 plant model. The random amplitudes of the RARITs were used as measures of input into the system, and the amplitudes of evoked CA1 fEPSPs were used as measures of the output of the system. The LV kernel model was applied to study the nonlinearity of the CA1 system. Examples of the first and the second order LV kernels are shown in Figures 5A,B, respectively.

FIGURE 5

Figure 5. (A) The first order and (B) the second order LV kernel of the CA1 plant model. The singular points and lines showing on the edge of each figure (indicated by arrows) represent the effect of present input. (C) A segment of comparison between RARIT-induced monosynaptic CA1 response amplitudes and the amplitudes predicted by the CA1 plant model. (D) The Q–Q plot of the data distribution between actual monosynaptic CA1 responses recorded from the slice and the outputs predicted by the CA1 plant model.

Model estimation was completed using stimulation intensities and intervals of the input–output sequences. The VAF between slice response and model prediction was 85.56 ±13.91%, averaged from six datasets. This shows that the CA1 plant model can accurately predict CA1 amplitudes based on stimulation amplitudes. A segment of the RARIT-induced monosynaptic CA1 fEPSP amplitudes is compared to its counterpart predicted by the CA1 plant model, as shown in Figure 5C. Figure 5D displays the Q–Q plot of the overall monosynaptic CA1 responses and CA1 plant model predicted results.

Inverse CA1 Plant Model Implementation and Validation Results

The implementation of the inverse CA1 plant model is accomplished using RARIT-induced monosynaptic data. The purpose for formulating such an inverse model is to transform output predictions into input stimulations. The output predictions were acquired from the CA1 plant model, and applied as the y in (5). The coefficients c₀, c₁, c₂,… were obtained during process of model estimation. Three terms involved the convolution of Laguerre basis functions and input amplitude were unknown, which include

${\displaystyle \sum _{l=1}^L\text{}}{c}_2(l)v_l(t_i),{\displaystyle \sum _{l=1}^L\text{}}{c}_1(l)v_l(t_i),\text{ and }{\displaystyle \sum _{l_1=1}^L{\displaystyle \sum _{l_2=1}^L{c}_2}}(l_1,l_2)v_{l_1}(t_i)v_{l_2}(t_i).$

Based on our experimental design, no stimulation existed before the stimulation train was sent, so these unknown terms were equal to zero. Thus, the first stimulation amplitude can be calculated by (7). After the first stimulation amplitude was calculated, it was then applied to convolve with the Laguerre basis function and formulate the unknown terms for calculating the next stimulation amplitude. The operations for solving the root were run through all data points in order to process the transformation from output into input. As a result, the inverse model allows us to convert the desired output response amplitudes to input stimulation amplitudes in a dynamic, recursive manner.

The validation of this inverse model was completed by comparing the calculated stimulation amplitudes with the RARIT amplitudes. The scatter plot in Figure 6 shows that the calculated stimulation amplitudes and the RARIT stimulation amplitudes are identical, showing that: (1) the real roots could all be calculated; and (2) this inverse model implementation can successfully derive optimal stimulations based on desired response amplitudes.

FIGURE 6

Figure 6. The comparison between stimulation amplitude in RARIT that sent into CA1 and the amplitude calculated from inverse CA1 plant model.

Modeling-Control Results

Following the protocol in this modeling-control framework experiment, CA1 desired output is first predicted through DG-CA1 trajectory model, and then applied into the inverse CA1 plant model to derive the optimal stimulation amplitudes. These amplitudes were used to formulate DARITs and were then sent into the slice, and the monosynaptic CA1 responses were recorded. The proposed modeling-control framework was intended to evoke CA1 to produce activities similar to the original CA1 activities. Thus, DARIT-induced monosynaptic CA1 amplitudes were compared with FARIT-induced trisynaptic CA1 response amplitudes. Two examples are shown in Figure 7.

FIGURE 7

Figure 7. A comparison of CA1 fEPSP amplitudes in response to FARIT-induced trisynaptic response (CA1-trysynaptic in blue diamonds) and DARIT-induced monosynaptic response (CA1-Controlled in red squares). The NMSE from 6 datasets was 15.41 ±8.35%.

Each panel illustrates results from one experiment: amplitudes of fEPSPs recorded from the CA1 region are shown as a function of 50 impulses chosen from among 2400 impulses of the stimulation trains (1200 administered with FARIT stimulation; 1200 administered with DARIT stimulation). In order to collapse the x axis to comprise more data points, time intervals between impulses are not represented in the figures; only “Input Event” number (sequence of sample impulses) is shown. Data for the FARIT-induced trisynaptic CA1 responses (CA1-trisynaptic) are shown in blue diamonds; data for the DARIT-induced monosynaptic CA1 responses (CA1-Controlled) are shown in red squares. As seen in Figure 7, the variation in CA1 fEPSP amplitudes was also captured in our model controlled paradigm. The accuracy was evaluated using NMSE of the amplitude, and the average NMSE was 15.41 ±8.35%. A Q–Q plot compared through the entire data sets is shown in Figure 8.

FIGURE 8

Figure 8. A Q–Q plot demonstration of the data comparison between the trisynaptic CA1 responses and the paradigm controlled CA1 responses.

Discussions

One of the essential objectives of a neural prosthetic device is to recreate the desired neural responses. While it is important to develop a reliable hardware model to represent the computational functions of a system, the control between device stimuli and actual responses is equally important. For example, in the application of deep brain stimulation (DBS), many efforts have been made in calibrating the stimulation parameters to achieve the desired effect (Mayberg et al., 2005; Okun et al., 2005). DBS devices depend on a trial-and-error process for finding the optimal stimulation pattern. Patients must repeatedly perform an exercise for a neurologist to adjust the stimulation parameters such as voltage, amplitude, pulse width, frequency, and electrode position (Moro et al., 2002; Volkmann et al., 2002; O'Suilleabhain et al., 2003). Another example is the application of functional electrical stimulation (FES) (Riener, 1999; Duffell et al., 2008; Donovan-Hall et al., 2011). The basic principle of FES is the generation of action potentials in the uninjured lower motor neurons by external electrical stimulation. This device faces problems such as muscle fatigue, spasticity, and limited force in the stimulated muscle. Using control strategies is one way to avoid internal disturbances and improve the time-consuming trial-and-error adjustment (Matjacic et al., 2003; Braz et al., 2009). Current neural prostheses face the same problem—the stimulation signals need to be adjusted manually or empirically based on the output response. In this article, we describe a rigorous approach to generate the stimulation patterns using a modeling-control framework. In the hippocampal slice preparation, with the purpose of restoring the CA1 output responses observed in the intact trisynaptic (DG to CA3 to CA1) circuitry, a nonlinear trajectory model was built to predict the CA1 desired output based on DG input. The predicted CA1 output was then converted to optimal stimulation through an inverse plant model of CA1 (i.e., an inverse transformation of CA1 input–output properties). Thus, the stimulation was essentially derived based on the desired output response, and was used to reactivate the CA1 response. An experimental validation of this modeling-control paradigm using hippocampal slices is provided. One of our preliminary studies was to stimulate CA1 region with non-optimal stimulation parameters, which means the nonlinearity of CA1 input–output relationship was not take into concern. The average NMSE from four experiments was 35.23 ±18.21%, which is much higher than the modeling-control paradigm results.

In current experimental paradigm, the only open parameters is the stimulation intensity, since the desired output responses are single EPSPs, the stimulations are the standard biphasic pulses, and there is no frequency (multiple impulses will elicit undesired multiple EPSPs). We are aware of the fact that in other applications (e.g., DBS), the phase and frequency are equally important parameters and also could be optimized. Our current modeling-control paradigm can be extended and used as a platform for the optimization of those parameters in the future. We understand that the stimulation site is also critical for this kind of devices, sometimes the misplacement of the electrode lead could cause poor efficacy or adverse effects (Richardson et al., 2009; van den Munckhof et al., 2010). Current clinical DBS surgeries were assisted with preoperative images analysis (MRI or CT), and intra-operatively guided with computerized stereotactic techniques. The lead could also be switched with limited-adjustability to compensate the inappropriate placement issue. By applying the stimulation electrode array in current experimental setup, the optimal stimulation site and its influence to the model can be further evaluated.

Our demonstrations also show that implementing a bidirectional neural prosthesis implicitly replaces the damaged system. In this approach, we do not need to explicitly estimate the transformational property of the CA3 region in the trisynaptic circuit. As long as we have the trajectory responses of the output system and once we identify the nonlinear input–output relationship of the output system, we are able to drive it to the desired output through its inverse model. One potential issue here would be “How to know the desired trajectory responses in the intact system?” In our opinion, there are several solutions/mechanisms that can mitigate this problem. First, we may develop a “generic model” from data recorded in normal animals. Previous studies have shown that there are significant amount of common features in the functional input–output properties across different animals, despite the animal-to-animal variations. In the case described in this study, all trajectory models are qualitatively similar in terms of the kernel polarity, kernel duration and kernel shape. Using model derived from other animals is imperfect, but at least provides a good approximation. Second, in behaving animal applications, the “imperfect” outputs generated by the “imperfect model” will be read out by the downstream brain regions. Neural plasticity, which is ubiquitous in the central nervous system, may play a role in adapting the system to the imperfect outputs or model. Third, more sophisticated computational methods such as reinforcement learning can potentially be used to develop self-adaptive or co-adaptive models.

The paradigm introduced in this paper did not include the error feedback. This was based on the assumption that the error observed in the output responses is instantaneous and does not influence the future output. In order to extend the paradigm to a closed-loop feedback system (Bernotas et al., 1986; Houk, 1988; Veltink et al., 1992; Abbas and Riener, 2001; Liu et al., 2011), the output error, that may caused by the interface between electrodes and nervous systems, the variation of the system, or the internal disturbance, need to be considered. The trajectory model developed in this study can be used as a reference model, which provide the desired output responses to compare with the actual responses recorded from output system, to calculate the error signal (Figure 9). In this case, the feedback error will be used as an external input signal and sent to adjust the properties of the inverse model. The influence of the previous errors on the current output may be taken in to account in a dynamic manner. In such a scheme, the optimal stimulation signals are calculated by the inverse plant model based on both the input and error signals.

FIGURE 9

Figure 9. A schematic diagram of a neuroprosthesis implemented with a feedforward and feedback controller.

Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

This work was supported in part by the NSF (BMES ERC and BITS Program), DARPA (HAND Project), ONR (Adaptive Neural System Program), NIBIB, and the Brain Restoration Foundation.

References

Abbas, J. J., and Riener, R. (2001). Using mathematical models and advanced control systems techniques to enhance neuroprosthesis function. Neuromodulation 4, 187–195.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Albert, G. C., Cook, C. M., Prato, F. S., and Thomas, A. W. (2009). Deep brain stimulation, vagal nerve stimulation and transcranial stimulation: an overview of stimulation parameters and neurotransmitter release. Neurosci. Biobehav. Rev. 33, 1042–1060.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Allison, B. Z., Wolpaw, E. W., and Wolpaw, J. R. (2007). Brain-computer interface systems: progress and prospects. Expert Rev. Med. Dev. 4, 463–474.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Amaral, D. G., and Witter, M. P. (1989). The three-dimensional organization of the hippocampal formation: a review of anatomical data. Neuroscience 31, 571–591.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Andersen, P., Soleng, A. F., and Raastad, M. (2000). The hippocampal lamella hypothesis revisited. Brain Res. 886, 165–171.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Andersen, P., Bliss, T. V., Lomo, T., Olsen, L. I., and Skrede, K. K. (1969). Lamellar organization of hippocampal excitatory pathways. Acta Physiol. Scand. 76, 4A–5A.

Pubmed Abstract | Pubmed Full Text

Bartesaghi, R., Migliore, M., and Gessi, T. (2006). Input-output relations in the entorhinal cortex-dentate-hippocampal system: evidence for a non-linear transfer of signals. Neuroscience 142, 247–265.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Berger, T. W., Ahuja, A., Courellis, S. H., Deadwyler, S. A., Erinjippurath, G., Gerhardt, G. A., et al. (2005). Restoring lost cognitive function. IEEE Eng. Med. Biol. Mag. 24, 30–44.

Pubmed Abstract | Pubmed Full Text

Berger, T. W., Baudry, M., Brinton, R. D., Liaw, J.-S., Marmarelis, V. Z., Yoondong Park, A., et al. (2001). Brain-implantable biomimetic electronics as the next era in neural prosthetics. Proc. IEEE 89, 993–1012.

Berger, T. W., Chauvet, G., and Sclabassi, R. J. (1994). A biologically based model of functional properties of the Hippocampus. Neural Netw. 7, 1031–1064.

Berger, T. W., Eriksson, J. L., Ciarolla, D. A., and Sclabassi, R. J. (1988). Nonlinear systems analysis of the hippocampal perforant path-dentate projection. III. Comparison of random train and paired impulse stimulation. J. Neurophysiol. 60, 1095–1109.

Pubmed Abstract | Pubmed Full Text

Berger, T. W., Hampson, R. E., Song, D., Goonawardena, A., Marmarelis, V. Z., and Deadwyler, S. (2011). A cortical neural prosthesis for restoring and enhancing memory. J. Neural Eng. 8, 046017.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Berger, T. W., Song, D., Chan, R. H. M., Marmarelis, V. Z., LaCoss, J., Wills, J., et al. (2012). A hippocampal cognitive prosthesis: multi-input, multi-output nonlinear modeling and VLSI implementation. IEEE Trans. Neural Syst. Rehabil. Eng. 20, 198–211.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Bernotas, L. A., Crago, P. E., and Chizeck, H. J. (1986). A discrete-time model of electrically stimulated muscle. IEEE Trans. Biomed. Eng. 33, 829–838.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Braz, G. P., Russold, M., and Davis, G. M. (2009). Functional electrical stimulation control of standing and stepping after spinal cord injury: a review of technical characteristics. Neuromodulation 12, 180–190.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Burgess, N., Maguire, E., and O'Keefe, J. (2002). The human hippocampus and spatial and episodic memory. Neuron 35, 625–641.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Camacho, E. F., and Bordons, C. (2003). Model Predictive Control. London, New York: Springer, xxii, 405.

Chan, C.-H., Srinivasan, V., Hsiao, M.-C., Khanna, S., Wills, J. D., Gholmieh, G., et al. (2004). “Real time hardware neural spike amplitude extraction,” in IEEE International Workshop on Biomedical Circuits and Systems 2004 (Singapore), 437–440.

Creasey, G. H., Ho, C. H., Triolo, R. J., Gater, D. R., DiMarco, A. F., Bogie, K. M., et al. (2004). Clinical applications of electrical stimulation after spinal cord injury. J. Spinal Cord Med. 27, 365–375.

Pubmed Abstract | Pubmed Full Text

Donoghue, J. P. (2002). Connecting cortex to machines: recent advances in brain interfaces. Nat. Neurosci. 5(Suppl), 1085–1088.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Donovan-Hall, M. K., Burridge, J., Dibb, B., Ellis-Hill, C., and Rushton, D. (2011). The views of people with spinal cord injury about the use of functional electrical stimulation. Artif. Organs 35, 204–211.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Duffell, L. D., Donaldson, N. D. N., Perkins, T. A., Rushton, D. N., Hunt, K. J., Kakebeeke, T. H., et al. (2008). Long-term intensive electrically stimulated cycling by spinal cord-injured people: effect on muscle properties and their relation to power output. Muscle Nerve 38, 1304–1311.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Gholmieh, G., Courellis, S., Dimoka, A., Wills, J. D., LaCoss, J., Granacki, J. J., et al. (2004). An algorithm for real-time extraction of population EPSP and population spike amplitudes from hippocampal field potential recordings. J. Neurosci. Methods 136, 111–121.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Gholmieh, G., Soussou, W., Han, M., Ahuja, A., Hsiao, M.-C., Song, D., et al. (2006). Custom-designed high-density conformal planar multielectrode arrays for brain slice electrophysiology. J. Neurosci. Methods 152, 116–129.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Hochberg, L. R., Serruya, M. D., Friehs, G. M., Mukand, J. A., Saleh, M., Caplan, A. H., et al. (2006). Neuronal ensemble control of prosthetic devices by a human with tetraplegia. Nature 442, 164–171.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Houk, J. C. (1988). Control strategies in physiological systems. FASEB J. 2, 97–107.

Pubmed Abstract | Pubmed Full Text

Hsiao, M.-C., Chan, C., Srinivasan, V., Ahuja, A., Erinjippurath, G., and Zanos, T. P. (2006) VLSI implementation of a nonlinear neuronal model: a ‘neural prosthesis’ to restore hippocampal trisynaptic dynamics. Conf. Proc. IEEE Eng. Med. Biol. Soc. 4396–4399.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Karniel, A., Meir, R., and Inbar, G. F. (2001). Best estimated inverse versus inverse of the best estimator. Neural Netw. 14, 1153–1159.

Pubmed Abstract | Pubmed Full Text

Lauer, R. T., Peckham, P. H., Kilgore, K. L., and Heetderks, W. J. (2000). Applications of cortical signals to neuroprosthetic control: a critical review. IEEE Trans. Rehabil. Eng. 8, 205–208.

Pubmed Abstract | Pubmed Full Text

Liu, J., Khalil, H. K., and Oweiss, K. G. (2011). Neural feedback for instantaneous spatiotemporal modulation of afferent pathways in bi-directional brain-machine interfaces. IEEE Trans. Neural Syst. Rehabil. Eng. 19, 521–533.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Liu, J., and Oweiss, K. (2010). “Feedback control of the spatiotemporal firing patterns of neural microcircuits,” in 49th IEEE Conference on Decision and Control (Atlanta, GA), 4679–4684.

Marmarelis, V. Z. (1993). Identification of nonlinear biological systems using Laguerre expansions of kernels. Ann. Biomed. Eng. 21, 573–589.

Pubmed Abstract | Pubmed Full Text

Marmarelis, V. Z. (2004). Nonlinear Dynamic Modeling of Physiological Systems. Hoboken, NJ: Wiley-Interscience, xvi, 541.

Marmarelis, V. Z., and Orme, M. E. (1993). Modeling of neural systems by use of neuronal modes. IEEE Trans. Biomed. Eng. 40, 1149–1158.

Pubmed Abstract | Pubmed Full Text

Matjacic, Z., Hunt, K., Gollee, H., and Sinkjaer, T. (2003). Control of posture with FES systems. Med. Eng. Phys. 25, 51–62.

Pubmed Abstract | Pubmed Full Text

Mayberg, H. S., Lozano, A. M., Voon, V., McNeely, H. E., Seminowicz, D., Hamani, C., et al. (2005). Deep brain stimulation for treatment-resistant depression. Neuron 45, 651–660.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

McIntyre, C. C., Savasta, M., Walter, B. L., and Vitek, J. L. (2004). How does deep brain stimulation work? Present understanding and future questions. J. Clin. Neurophysiol. 21, 40–50.

Pubmed Abstract | Pubmed Full Text

McLachlan, R. S., Pigott, S., Tellez-Zenteno, J. F., Wiebe, S., and Parrent, A. (2010). Bilateral hippocampal stimulation for intractable temporal lobe epilepsy: impact on seizures and memory. Epilepsia 51, 304–307.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Middlebrooks, J. C., Bierer, J. A., and Snyder, R. L. (2005). Cochlear implants: the view from the brain. Curr. Opin. Neurobiol. 15, 488–493.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Moro, E., Esselink, R. J. A., and Xie, J. (2002). The impact on Parkinson's disease of electrical parameter settings in STN stimulation. Neurology 59, 706–713.

Pubmed Abstract | Pubmed Full Text

Nicolelis, M. A. L. (2003). Brain-machine interfaces to restore motor function and probe neural circuits. Nat. Rev. Neurosci. 4, 417–422.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Normann, R. A. (2007). Technology insight: future neuroprosthetic therapies for disorders of the nervous system. Nat. Clin. Pract. Neurol. 3, 444–452.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Okun, M. S., Tagliati, M., Pourfar, M., Fernandez, H. H., Rodriguez, R. L., Alterman, R. L., et al. (2005). Management of referred deep brain stimulation failures: a retrospective analysis from 2 movement disorders centers. Arch. Neurol. 62, 1250–1255.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

O'Suilleabhain, P. E., Frawley, W., Giller, C., and Dewey, R. B. (2003). Tremor response to polarity, voltage, pulsewidth and frequency of thalamic stimulation. Neurology 60, 786–790.

Pubmed Abstract | Pubmed Full Text

Patil, P. G., and Turner, D. A. (2008). The development of brain-machine interface neuroprosthetic devices. Neurotherapeutics 5, 137–146.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Richardson, R. M., Ostrem, J. L., and Starr, P. (2009). Surgical repositioning of misplaced subthalamic electrodes in Parkinson's disease: location of effective and ineffective leads. Stereotact. Funct. Neurosurg. 87, 297–303.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Riener, R. (1999). Model-based development of neuroprosthesis for paraplegic patients. Philos. Trans. R. Soc. Lond. B Biol. Sci. 354, 877–894.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Rupp, R., and Gerner, H. J. (2007). Neuroprosthetics of the upper extremity–clinical application in spinal cord injury and challenges for the future. Acta Neurochir. Suppl. 97(Pt 1), 419–426.

Pubmed Abstract | Pubmed Full Text

Schwartz, A. B. (2004). Cortical neural prosthetics. Annu. Rev. Neurosci. 27, 487–507.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Sclabassi, R. J., Eriksson, J. L., Port, R. L., Robinson, G. B., and Berger, T. W. (1988). Nonlinear systems analysis of the hippocampal perforant path-dentate projection. I. Theoretical and interpretational considerations. J. Neurophysiol. 60, 1066–1076.

Pubmed Abstract | Pubmed Full Text

Song, D., Chan, R. H. M., Marmarelis, V. Z., Hampson, R. E., Deadwyler, S. A., and Berger, T. W. (2007). Nonlinear dynamic modeling of spike train transformations for hippocampal-cortical prostheses. IEEE Trans. Biomed. Eng. 54(6 Pt 1), 1053–1066.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Song, D., Chan, R. H. M., Marmarelis, V. Z., Hampson, R. E., Deadwyler, S. A., and Berger, T. W. (2009). Nonlinear modeling of neural population dynamics for hippocampal prostheses. Neural Netw. 22, 1340–1351.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Squire, L. R. (1992). Memory and the hippocampus: a synthesis from findings with rats, monkeys, and humans. Psychol. Rev. 99, 195–231.

Pubmed Abstract | Pubmed Full Text

Stacey, W. C., and Litt, B. (2008). Technology insight: neuroengineering and epilepsy-designing devices for seizure control. Nat. Clin. Pract. Neurol. 4, 190–201.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Taketani, M., and Baurdy, M. (2006). Advances in Network Electrophysiology: Using Multi-Electrode Arrays, 1st Edn. Michigan: Springer, 127–152.

Taylor, D. M., Tillery, S. I. H., and Schwartz, A. B. (2003). Information conveyed through brain-control: cursor versus robot. IEEE Trans. Neural Syst. Rehabil. Eng. 11, 195–199.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Tellez-Zenteno, J. F., McLachlan, R. S., Parrent, A., Kubu, C. S., and Wiebe, S. (2006). Hippocampal electrical stimulation in mesial temporal lobe epilepsy. Neurology 66, 1490–1494.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Turner, D. A., Dimitrov, D. F., and Ncolelis, M. A. L. (2005). “Neuroprosthetics and clinical realization of brain-machine interfaces, Chapter 7,” in Modern neurosurgery: Clinical Translation of Neuroscience Advances, 1st Edn. ed D. A. Turner (Danvers, MA: CRC Press), 117–138.

van den Munckhof, P., Contarino, M. F., Bour, L. J., Speelman, J. D., de Bie, R. M., and Schuurman, P. R. (2010). Postoperative curving and upward displacement of deep brain stimulation electrodes caused by brain shift. Neurosurgery 67, 49–53; discussion 53–54.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Veltink, P. H., Chizeck, H. J., Crago, P. E., and El-Bialy, A. (1992). Nonlinear joint angle control for artificially stimulated muscle. IEEE Trans. Biomed. Eng. 39, 368–380.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Volkmann, J., Herzog, J., Kopper, F., and Deuschl, G. (2002). Introduction to the programming of deep brain stimulators. Mov. Disord. 17(Suppl 3), S181–S187.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Weiland, J. D., Liu, W., and Humayun, M. S. (2005). Retinal prosthesis. Annu. Rev. Biomed. Eng. 7, 361–401.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Widrow, B., and Bilello (1993). “Adaptive inverse control,” in Proceedings of the 1993 IEEE International Symposium on Intelligent Control 1993 (Chicago, IL), 1–6.

Widrow, B., and Plett, G. L. (1997). “Nonlinear adaptive inverse control,” in Proceedings of the 36th IEEE Conference on Decision and Control 1997, Vol. 2 (San Diego, CA), 1032–1037.

Widrow, B., and Walach, E. (1996). Adaptive Inverse Control. Upper Saddle River, NJ: Prentice Hall PTR, xvii, 508.