Abstract:
A DBS system and method for predicting future neurological activity in a subject and administering a corrective electrical stimulation signal to prevent anticipated pathological neuronal activity. The DBS system includes an implantable electrode configured to both record neuronal activity from a target brain area in a subject and administer the corrective electric stimulation signal to the target area. The DBS system also includes a controller configured to determine the characteristics of the corrective electrical stimulation signal based on point process models of healthy and pathological neuronal activity in the target area.

Description:
CROSS-REFERENCE TO RELATED APPLICATIONS 
     This application is based on, claims the benefit of, and incorporates herein by reference in their entirety, PCT International Application PCT/US2009/062072 filed on Oct. 26, 2009 and U.S. Provisional Patent Application Ser. No. 61/108,060 filed on Oct. 24, 2008, and entitled “SYSTEM AND METHOD FOR DYNAMICALLY CONFIGURABLE DEEP BRAIN STIMULATION.” 
    
    
     BACKGROUND OF THE INVENTION 
     Neurodegenerative diseases plague a large number of individuals world wide. For example, an estimated 3-4 million people in the US have Parkinson&#39;s disease (PD), which is a chronic progressive neurodegenerative disease that occurs when dopaminergic neurons in the substantia nigra pars compacta of the midbrain degenerate, causing resting tremor, rigidity, and bradykinesia. Currently, there is no cure or definitive means to stop the progression of many neurodegenerative diseases and PD is just one example. However, medications and surgery can relieve many of the symptoms. Such treatments for PD have been developed based on an improved understanding of basal ganglia (BG) anatomy and physiology. 
     It has been long appreciated that PD follows the degeneration of dopaminergic neurons in the substantia nigra pars compacta. This triggers a cascade of functional changes in the BG that leads to abnormal activity in its output nuclei, the substantia nigra reticulata and globus pallidus internus. Therefore, the goal of traditional treatment is to enhance concentrations of dopamine, or to modify the activity of the output nuclei by creating lesions in target areas, or more recently by using deep brain stimulation (DBS). 
     DBS is a surgical procedure in which a stimulating probe is implanted in a targeted area, typically the subthalamic nucleus (STN), which is connected to an insulated wire that is passed under the skin of the head, neck, and shoulder and terminated at a neurostimulator, which typically sits inferior to the clavicle. At major surgical centers the surgery has become routine. Patient&#39;s motor symptoms generally decrease with treatment and they can regain quality of life and reduce their medications, which have serious side effects. 
     While traditional DBS is a valuable tool for treating neurological disorders, the stimulation signal must be optimized post-operatively.  FIG. 1  provides schematic depiction of the “open loop” feedback system of traditional DBS systems, in which both natural external stimulus x(t) and an administered DBS signal u(t) affect the neuronal activity y(t) of a target region in a subject&#39;s brain. The neuronal activity y(t) in turn affects the behavior b(t) of the subject&#39;s body, which feeds back to affect the target region. The DBS system is not responsive to feedback from either y(t) or b(t). Therefore, identification of a stimulation signal u(t) that elicits an appropriate subject response is achieved by a manual calibration process in which a multitude of different stimulation signals are evaluated by trial-and-error. Manual calibration is costly in terms of medical resources and can limit the number of patients that a neurologist may treat simultaneously. Calibration can take several weeks or months, during which the high expectations of the subject having an implanted DBS system may not be met, thereby leading to subject depression and medication use. Moreover, current DBS waveforms are high-frequency and can leak into neighboring brain areas, causing side effects. Beyond such undesirable side effects, current DBS systems, employing high-frequency stimulation waveforms, consume significant amounts of power and rapidly drain batteries that must be replaced surgically. 
     It would therefore be desirable to have DBS system that would reduce the resources needed for calibration and operation and provide more immediate and effective treatment to the patient. Such a system would improve patient care, reduce medical costs, increase the number of patients that a neurologist may treat simultaneously, and improve the performance characteristics, for example, battery life, of the device. 
     SUMMARY OF THE INVENTION 
     The present invention overcomes the aforementioned drawbacks by providing a system and method for performing deep brain stimulation that is capable of dynamic and continuous self-calibration. The system for performing deep brain stimulation includes an implantable electrode configured to measure a physiological parameter of a subject and deliver an electrical stimulation signal to a target area in the subject. The DBS system also includes a neurostimulator having a memory device having stored thereon a healthy model characterizing healthy neuronal activity in the target area and a pathological model characterizing pathological neuronal activity in the target area. The neurostimulator also includes a controller in communication with the implantable electrode and the memory device and configured to analyze the measured physiological parameter using the healthy model and the pathological model to identify a corrective electrical stimulation signal that, when delivered by the implantable electrode to the target area, reduces pathological neuronal events in the target area. 
     The method for performing deep brain stimulation includes the steps of measuring a physiological parameter from a target area of a subject using an electrode implanted in the subject and analyzing the measured physiological parameter using a healthy model characterizing healthy neuronal activity in the target area and a pathological model characterizing pathological neuronal activity in the target area to identify a corrective electrical stimulation signal that, when delivered to the target area, reduces pathological neural events in the target area. The method further includes the step of administering the corrective electrical stimulation signal to the target area using the implanted electrode. 
     Various other features of the present invention will be made apparent from the following detailed description and the drawings. 
    
    
     
       BRIEF DESCRIPTION OF THE DRAWINGS 
         FIG. 1  is a block diagram depicting a prior-art DBS system; 
         FIG. 2  is a schematic of a DBS apparatus in accordance with the present invention; 
         FIG. 3  is a flowchart setting forth the process used to create models for healthy and diseased subjects in accordance with the present invention; 
         FIG. 4  is a block diagram of a DBS system in accordance with the present invention; and 
         FIG. 5  is graph showing a segment of optimal solutions for p*(t) and q*(t) in accordance with the present invention. 
     
    
    
     DETAILED DESCRIPTION OF THE PREFERRED EMBODIMENT 
     Referring to  FIG. 2 , the present invention provides a DBS system  10  capable of predicting future neural activity and administering electrical stimulation to a target brain area, such as the STN, to prevent anticipated pathological neural activity. The DBS system  10  includes an electrode probe  12  that is capable of both stimulating populations of neurons and measuring single-unit neuronal activity. The probe  12  is typically implanted in the target and connected via an insulated lead  14  to a neurostimulator  16 . It is contemplated that the lead  14  runs under the skin of the head, neck, and shoulder and that the neurostimulator  16  is implanted to sit inferior to the clavicle. The neurostimulator  16  includes a pulse generator  18 , a controller  20 , and a battery pack  22  to power the DBS system  10 . The neurostimulator  16  can also include a memory to store measured neural activity data and models for implementation on the controller  20 . 
     In operation, the DBS system  10  acquires neuronal activity, or spike train, data with the electrode probe  12 . This neuronal activity data is carried via lead  14  to the neurostimulator  16  where it is processed by the controller  20 . The controller  20  analyzes the data and identifies a corrective stimulation signal that will prevent anticipated pathological neural events. The selected stimulation signal is then generated by the pulse generator  18  and delivered via the lead  14  to the electrode probe  12 , which administers the stimulation signal to the target area. Depending on the predicted neural activity, the stimulation signal may inhibit neurons, excite neurons, or do nothing. 
     Referring now to  FIGS. 3 and 4 , to operate as described above, the controller  20  employs control algorithms based on models that characterize the neural activity of healthy and PD subjects. Typically, the models relate neural spike trains for both healthy and diseased subjects to natural extrinsic factors, such as the environment x(t) and behavioral stimuli b(t), and intrinsic dynamics, such as the history of neural spiking activity y(t) and local neural network activity. A method for constructing such models begins at process block  52 , with the definition of functions characterizing neural spike trains in healthy and PD subjects. Specifically, a PD patient model y(t) and a healthy model y H (t) can be defined as follows:
 
 y ( t )= f   PD ( t,X   t   ,Y   t   ,u ( t ))  Eqn. 1;
 
and
 
 y   H ( t )= f   PD ( t,X   t   ,B   y   ,Y   t   H )  Eqn. 2;
         where y(t) is the scalar value of the spike train (generally, either 0 or 1) at a given time t, x(t) is the value of natural extrinsic stimuli, and b(t) is the value of behavioral stimuli. Y t , Y H   t , X t , and B t  are the sequences of historical values over the time interval from 0 to t for y(t), j H (t), x(t) and b(t), respectively. The variables Y t , X t , and B t  may be grouped into H t , which therefore gives the history of spiking activity and its covariates up to time t. At process block  54 , neurophysiological data for the construction of functions f PD  and f H  is acquired by recording neural activity from the same target brain area of both healthy and PD subjects that are executing the same behavioral task. For example, the neuronal activity of subjects having PD can collected in an operating room by functional neurosurgeons during the implantation of a DBS device, while healthy subject data can be collected from primates under similar experimental conditions.       

     At process block  56 , the collected healthy and PD neural activity data are fit to models to provide a basis for estimating of y(t) and y H (t). This can be achieved using a point process paradigm that overcomes difficulties in characterizing neural activity dynamics, particularly those associated with noise and the dependence of neural activity on both intrinsic and extrinsic factors. A point process is a binary stochastic process defined in continuous time, for example, the number of neuronal spikes in a given time interval, and is characterized entirely by a conditional intensity function (CIF). A CIF for a point process model relating the spiking propensity of target brain area neurons to factors associated with environmental conditions, behavioral stimuli, and the neurons&#39; spiking history can be defined as follows: 
     
       
         
           
             
               
                 
                   
                     
                       λ 
                       ⁡ 
                       
                         ( 
                         
                           
                             t 
                             ❘ 
                             
                               H 
                               t 
                             
                           
                           , 
                           
                             u 
                             ⁡ 
                             
                               ( 
                               t 
                               ) 
                             
                           
                         
                         ) 
                       
                     
                     = 
                     
                       
                         lim 
                         
                           Δ 
                           → 
                           0 
                         
                       
                       ⁢ 
                       
                         
                           Pr 
                           ⁡ 
                           
                             ( 
                             
                               
                                 
                                   
                                     N 
                                     ⁡ 
                                     
                                       ( 
                                       
                                         t 
                                         + 
                                         Δ 
                                       
                                       ) 
                                     
                                   
                                   - 
                                   
                                     N 
                                     ⁡ 
                                     
                                       ( 
                                       t 
                                       ) 
                                     
                                   
                                 
                                 = 
                                 
                                   1 
                                   ❘ 
                                   
                                     H 
                                     t 
                                   
                                 
                               
                               , 
                               
                                 u 
                                 ⁡ 
                                 
                                   ( 
                                   t 
                                   ) 
                                 
                               
                             
                             ) 
                           
                         
                         Δ 
                       
                     
                   
                   ; 
                 
               
               
                 
                   Eqn 
                   . 
                   
                       
                   
                   ⁢ 
                   3 
                 
               
             
           
         
       
     
     where N(t) is the number of spikes in a time interval [0,t] for tε(0,T] and t 1  to t n  denotes the time of measured neuronal spikes such that 0&lt;t 1 &lt;t 2 &lt; . . . &lt;t n ≦T. Multivariate point process models may also be employed in accordance with the present invention, for example, λ(t|H t , u(t)) and N(t) may be vector-valued when modeling several neurons. Accordingly, Eqn. 3 defines the probability of a spike in each neuron in any small time interval (t,t+Δ) as follows:
 
 Pr (spike in( t,t +Δ)| H   t )≈λ( t|H   t )Δ  Eqn. 4.
 
     Thus, when Δ is small, Eqn. 4 is approximately equal to the spiking propensity at time t. While the model can be fitted to measured neuronal activity data using any appropriate parametric or nonparametric modeling class, it is contemplated that the present invention employs a generalized linear model (GLM) framework. In a GLM, the log of the CIF is a linear function of model parameters. A GLM is advantageous because it separates the contributions of extrinsic and intrinsic factors to the probability that the neuron will spike at a given time t. A GLM also provides an efficient computational scheme for estimating model parameter and a likelihood framework for conducting statistical inferences based on the estimated model. For example, a GLM framework for fitting collect neuronal activity data to the above CIF employ the following relations: 
     
       
         
           
             
               
                 
                   
                     
                       λ 
                       ⁡ 
                       
                         ( 
                         
                           
                             t 
                             ❘ 
                             
                               H 
                               t 
                             
                           
                           , 
                           
                             u 
                             ⁡ 
                             
                               ( 
                               t 
                               ) 
                             
                           
                           , 
                           θ 
                         
                         ) 
                       
                     
                     = 
                     
                       
                         θ 
                         0 
                       
                       + 
                       
                         
                           ∑ 
                           
                             i 
                             = 
                             1 
                           
                           I 
                         
                         ⁢ 
                         
                             
                         
                         ⁢ 
                         
                           
                             α 
                             i 
                           
                           ⁢ 
                           
                             
                               f 
                               i 
                             
                             ⁡ 
                             
                               ( 
                               
                                 
                                   extrinsic 
                                   ⁢ 
                                   
                                       
                                   
                                   ⁢ 
                                   covariates 
                                 
                                 , 
                                 t 
                               
                               ) 
                             
                           
                         
                       
                       + 
                       
                         
                           ∑ 
                           
                             j 
                             = 
                             1 
                           
                           J 
                         
                         ⁢ 
                         
                             
                         
                         ⁢ 
                         
                           
                             β 
                             j 
                           
                           ⁢ 
                           
                             
                               g 
                               j 
                             
                             ⁡ 
                             
                               ( 
                               
                                 
                                   intrinsic 
                                   ⁢ 
                                   
                                       
                                   
                                   ⁢ 
                                   covariates 
                                 
                                 , 
                                 t 
                               
                               ) 
                             
                           
                         
                       
                     
                   
                   ; 
                 
               
               
                 
                   Eqn 
                   . 
                   
                       
                   
                   ⁢ 
                   5 
                 
               
             
           
         
       
         
         
           
             where the vector θ 0 ={α, β, γ} denotes the unknown parameters, such as spike times, to be estimated from the collected data. In this equation, the logarithm of the CIF is linear for θ 0  and f i  and g j  are basis functions and may be any arbitrary nonlinear function. It should be noted that this GLM is an extension of a multiple linear regression model in which the spike times being estimated need not be Gaussian. 
           
         
       
    
     At process block  60 , the goodness-of-fit of the point process model can optionally be tested. This can be achieved by generating a Kolmogorov-Smirov (KS) plot that compares the empirical cumulative distribution function of time-scaled spike times to the cumulative distribution function of a unit rate exponential. Improved goodness-of-fit is indicated if the KS plot lies on the 45 degree line. Further, a 95 percent confidence bounds can be computed for the degree of agreement using the distribution of the KS statistic. To test the independence of rescaled times, the spike times can be transformed into Gaussian rescaled times with zero means and unit variances. Since lack of correlation is equivalent to independence for Gaussian random variables, the autocorrelation function (ACF) of the Gaussian rescaled times can be plotted and the number of points of the ACF lying outside the 95% confidence intervals can be counted. 
     At process block  62 , the point process models can be related back to the functions y H (t) and y(t). Using time units of msecs and assuming Δ=1, this can be by noting that at any time t, f PD  and f H  are random variables that take on the values 0 or 1 according to the following probabilities:
 
Prob[ f   PD ( t,H   t   ,u ( t ))=1]≈λ PD ( t|H   t   ,u ( t ))  Eqn. 6;
 
and
 
Prob[ f   H ( t,H   t   ,u ( t ))=1)]≈λ H ( t|H   t )  Eqn. 7.
 
     Initial point process models characterizing healthy and PD STN neuronal activity in the absence of DBS stimulation, that is, λ PD (t|H t ) and λ H (t|H t ), can be generated from neuronal activity recordings from PD subjects and healthy primates by assuming u(t)=0. The primates are used as surrogates for healthy humans and the studies are performed under identical conditions, for example, as the PD subjects and primates perform the same task. Studies using such models can quantify prevalent abnormalities in PD activity not present in healthy activity. In particular, the neural activity of PD subjects exhibits 10-30 Hz oscillations, bursting, and persistent directional tuning, all of which may directly related to the well-known PD motor symptoms of resting tremor, bradykinesia, and rigidity. However, for the control algorithm used in the DBS system  10  of  FIG. 2 , it is beneficial to have a predictive model that characterizes the effect of different DBS inputs on spiking activity. 
     The model λ PD (t|H t ,u(t)) can be predicted from λ PD (t|H t ) by making the following assumptions: First, the DBS signal u(t) is a sequence of the values 0, 1, −1, that is, the DBS signal is an aperiodic train of positive and negative pulses, where a value of 1 indicates a pulse with positive height and a value of −1 indicates a pulse with negative height. Specifically, the DBS signal is a time sequence of independent distributed random variable with the following probability distributions:
 
 Pr ( u ( t )=1)= p ( t )  Eqn. 8;
 
 Pr ( u ( t )=0)= q ( t )  Eqn. 9;
 
and
 
 Pr ( u ( t )=−1)=1−( q ( t )+ p ( t ))  Eqn. 10.
 
     Secondly, when the DBS signal u(t) is applied, it adds to the neuronal spike train of each STN neuron it impacts. If the spike train of a single neuron can be denoted as y(t) with a CIF λ PD (t|H t ), then this addition obeys the following rules:
         a. If for a given time t, y(t)=0 and u(t)=1, then y(t)+u(t)=1 (event A)   b. If for a given time t, y(t)=0 and u(t)=0, then y(t)+u(t)=0   c. If for a given time t, y(t)=0 and u(t)=−1, then y(t)+u(t)=0   d. If for a given time t, y(t)=1 and u(t)=1, then y(t)+u(t)=1 (event A)   e. If for a given time t, y(t)=1 and u(t)=0, then y(t)+u(t)=1 (event A)   f. If for a given time t, y(t)=1 and u(t)=−1, then y(t)+u(t)=0       

     Lastly, the new spiking activity z(t)=y(t)+u(t) is a binary sequence of 0&#39;s and 1&#39;s and is point process with the following CIF: 
     
       
         
           
             
               
                 
                   
                     
                       
                         
                           
                             λ 
                             z 
                           
                           ⁡ 
                           
                             ( 
                             
                               t 
                               ❘ 
                               
                                 H 
                                 t 
                               
                             
                             ) 
                           
                         
                         = 
                           
                         ⁢ 
                         
                           
                             P 
                             ⁡ 
                             
                               ( 
                               
                                 
                                   z 
                                   ⁡ 
                                   
                                     ( 
                                     t 
                                     ) 
                                   
                                 
                                 = 
                                 
                                   1 
                                   ❘ 
                                   
                                     H 
                                     t 
                                   
                                 
                               
                               ) 
                             
                           
                           ≡ 
                           
                             
                               λ 
                               PD 
                             
                             ⁡ 
                             
                               ( 
                               
                                 
                                   t 
                                   ❘ 
                                   
                                     H 
                                     t 
                                   
                                 
                                 , 
                                 
                                   u 
                                   ⁡ 
                                   
                                     ( 
                                     t 
                                     ) 
                                   
                                 
                               
                               ) 
                             
                           
                         
                       
                     
                   
                   
                     
                       
                         = 
                           
                         ⁢ 
                         
                           P 
                           ⁡ 
                           
                             ( 
                             
                               
                                 event 
                                 ⁢ 
                                 
                                     
                                 
                                 ⁢ 
                                 a 
                               
                               ❘ 
                               
                                 H 
                                 t 
                               
                             
                             ) 
                           
                         
                       
                     
                   
                   
                     
                       
                         = 
                           
                         ⁢ 
                         
                           
                             
                               p 
                               ⁡ 
                               
                                 ( 
                                 t 
                                 ) 
                               
                             
                             ⁢ 
                             
                               ( 
                               
                                 1 
                                 - 
                                 
                                   
                                     λ 
                                     PD 
                                   
                                   ⁡ 
                                   
                                     ( 
                                     
                                       t 
                                       ❘ 
                                       
                                         H 
                                         t 
                                       
                                     
                                     ) 
                                   
                                 
                               
                               ) 
                             
                           
                           + 
                           
                             
                               p 
                               ⁡ 
                               
                                 ( 
                                 t 
                                 ) 
                               
                             
                             ⁢ 
                             
                               
                                 λ 
                                 PD 
                               
                               ⁡ 
                               
                                 ( 
                                 
                                   t 
                                   ❘ 
                                   
                                     H 
                                     t 
                                   
                                 
                                 ) 
                               
                             
                           
                           + 
                         
                       
                     
                   
                   
                     
                       
                           
                         ⁢ 
                         
                           
                             q 
                             ⁡ 
                             
                               ( 
                               t 
                               ) 
                             
                           
                           ⁢ 
                           
                             
                               λ 
                               PD 
                             
                             ⁡ 
                             
                               ( 
                               
                                 t 
                                 ❘ 
                                 
                                   H 
                                   t 
                                 
                               
                               ) 
                             
                           
                         
                       
                     
                   
                   
                     
                       
                         
                           = 
                             
                           ⁢ 
                           
                             
                               p 
                               ⁡ 
                               
                                 ( 
                                 t 
                                 ) 
                               
                             
                             + 
                             
                               
                                 q 
                                 ⁡ 
                                 
                                   ( 
                                   t 
                                   ) 
                                 
                               
                               ⁢ 
                               
                                 
                                   λ 
                                   PD 
                                 
                                 ⁡ 
                                 
                                   ( 
                                   
                                     t 
                                     ❘ 
                                     
                                       H 
                                       t 
                                     
                                   
                                   ) 
                                 
                               
                             
                           
                         
                         ; 
                       
                     
                   
                 
               
               
                 
                   Eqn 
                   . 
                   
                       
                   
                   ⁢ 
                   11 
                 
               
             
           
         
       
         
         
           
             where H t  is the history of y(t) and any other extrinsic factors up to time t, that is H t ={X t , B t , Y t }. 
           
         
       
    
     Referring now to  FIGS. 2 and 4 , healthy and PD spike train models generated using the above method can serve as the basis for controller algorithms that predict future neuronal activity based on measured factors and determined a corrective electrical stimulation signals to correct for anticipated pathological neuronal activity. As shown schematically in  FIG. 4 , the controller  20  may use such a model as a basis for 1) estimating environmental factors and behavioral stimuli from spike train measurements, 2) determining in real-time the DBS stimulation signal that minimizes a difference between a PD subject&#39;s spike train and the predicted spike train of a healthy subject, and 3) translating this information into continuous voltage signal to be generated by the pulse generator  18  and delivered to the target brain area by the electrode  12 . In contrast to traditional DBS systems, such as that depicted of  FIG. 1 , the present invention can be considered as a “closed loop” system, since controller determines the corrective stimulation signal u(t) based on environmental and behavioral stimuli, x(t) and b(t), respectively, as well as the neural activity of the target region y(t). 
     The environment and behavioral stimuli can be estimated from neuronal data acquired by the electrode  12  by maximizing the following a posteriori probability:
 
 Pr (stimuli|neuronal spiking activity up to time  t )  Eqn. 12;
         though it should be noted that this process may be skipped if the CIF is not a function of unknown environment and behavioral stimuli. The real-time determination of the DBS stimulation signal can be achieved by optimizing a distance function. The high-level control objective of the present invention is to match the neural activity of a PD subject in which the DBS system  10  of  FIG. 2  is implanted to that of a healthy representative. With the point-process paradigm, this translates to matching the CIF of the PD subject with that of a healthy subject. This match between the two CIFs can be achieved by minimizing the following distance relationship:
 
distance{λ PD ( t|H   t   ,u (t)),λ Healthy ( t|H   t )}  Eqn. 13.
       

     Referring to  FIG. 5 , it should be noted that the distance can be expressed using any relevant function. For example, the distance function may be selected to equate the average of each spike train in the least-squares sense by setting the distance as the expected value of the square difference of the PD spike train with DBS input (z(t)=y(t)+u(t)) and a healthy spike train h(t) as follows: 
     
       
         
           
             
               
                 
                   
                     distance 
                     ⁢ 
                     
                       { 
                       
                         
                           
                             λ 
                             PD 
                           
                           ⁡ 
                           
                             ( 
                             
                               
                                 t 
                                 ❘ 
                                 
                                   H 
                                   t 
                                 
                               
                               , 
                               
                                 u 
                                 ⁡ 
                                 
                                   ( 
                                   t 
                                   ) 
                                 
                               
                             
                             ) 
                           
                         
                         , 
                         
                           
                             λ 
                             H 
                           
                           ⁡ 
                           
                             ( 
                             
                               t 
                               ❘ 
                               
                                 H 
                                 t 
                               
                             
                             ) 
                           
                         
                       
                       } 
                     
                   
                   = 
                   
                     E 
                     ⁢ 
                     
                       
                         { 
                         
                           
                             [ 
                             
                               
                                 
                                   z 
                                   ⁡ 
                                   
                                     ( 
                                     t 
                                     ) 
                                   
                                 
                                 - 
                                 
                                   h 
                                   ⁡ 
                                   
                                     ( 
                                     t 
                                     ) 
                                   
                                 
                               
                               ❘ 
                               
                                 H 
                                 t 
                               
                             
                             ] 
                           
                           2 
                         
                         } 
                       
                       . 
                     
                   
                 
               
               
                 
                   Eqn 
                   . 
                   
                       
                   
                   ⁢ 
                   14 
                 
               
             
           
         
       
     
     Thus, the solution to the least-squares problem: 
     
       
         
           
             
               
                 
                   
                     
                       
                         min 
                         
                           u 
                           ⁡ 
                           
                             ( 
                             t 
                             ) 
                           
                         
                       
                       ⁢ 
                       
                         ( 
                         
                           E 
                           ⁢ 
                           
                             { 
                             
                               
                                 [ 
                                 
                                   
                                     
                                       z 
                                       ⁡ 
                                       
                                         ( 
                                         t 
                                         ) 
                                       
                                     
                                     - 
                                     
                                       h 
                                       ⁡ 
                                       
                                         ( 
                                         t 
                                         ) 
                                       
                                     
                                   
                                   ❘ 
                                   
                                     H 
                                     t 
                                   
                                 
                                 ] 
                               
                               2 
                             
                             } 
                           
                         
                         ) 
                       
                     
                     = 
                     
                       
                         min 
                         
                           u 
                           ⁡ 
                           
                             ( 
                             t 
                             ) 
                           
                         
                       
                       ⁢ 
                       
                         ( 
                         
                           E 
                           ⁢ 
                           
                             { 
                             
                               
                                 [ 
                                 
                                   
                                     
                                       y 
                                       ⁡ 
                                       
                                         ( 
                                         t 
                                         ) 
                                       
                                     
                                     + 
                                     
                                       u 
                                       ⁡ 
                                       
                                         ( 
                                         t 
                                         ) 
                                       
                                     
                                     - 
                                     
                                       h 
                                       ⁡ 
                                       
                                         ( 
                                         t 
                                         ) 
                                       
                                     
                                   
                                   ❘ 
                                   
                                     H 
                                     t 
                                   
                                 
                                 ] 
                               
                               2 
                             
                             } 
                           
                         
                         ) 
                       
                     
                   
                   ; 
                 
               
               
                 
                   Eqn 
                   . 
                   
                       
                   
                   ⁢ 
                   15 
                 
               
             
             
               
                 
                   
                       
                   
                   ⁢ 
                   
                     is 
                     ⁢ 
                     
                         
                     
                     ⁢ 
                     given 
                     ⁢ 
                     
                         
                     
                     ⁢ 
                     by 
                     ⁢ 
                     
                       : 
                     
                   
                 
               
               
                 
                     
                 
               
             
             
               
                 
                   
                     
                       
                         
                             
                         
                         ⁢ 
                         
                           
                             E 
                             ⁢ 
                             
                               { 
                               
                                 z 
                                 ⁡ 
                                 
                                   ( 
                                   t 
                                   ) 
                                 
                               
                               } 
                             
                           
                           = 
                           
                             
                               E 
                               ⁢ 
                               
                                 { 
                                 
                                   y 
                                   ⁡ 
                                   
                                     ( 
                                     t 
                                     ) 
                                   
                                 
                                 } 
                               
                             
                             + 
                             
                               E 
                               ⁢ 
                               
                                 { 
                                 
                                   
                                     u 
                                     * 
                                   
                                   ⁡ 
                                   
                                     ( 
                                     t 
                                     ) 
                                   
                                 
                                 } 
                               
                             
                           
                         
                       
                     
                   
                   
                     
                       
                         = 
                         
                           
                             E 
                             ⁢ 
                             
                               { 
                               
                                 h 
                                 ⁡ 
                                 
                                   ( 
                                   t 
                                   ) 
                                 
                               
                               } 
                             
                           
                           ⇒ 
                           
                             
                               λ 
                               PD 
                             
                             ⁡ 
                             
                               ( 
                               
                                 
                                   t 
                                   ❘ 
                                   
                                     H 
                                     t 
                                   
                                 
                                 , 
                                 
                                   
                                     u 
                                     * 
                                   
                                   ⁡ 
                                   
                                     ( 
                                     t 
                                     ) 
                                   
                                 
                               
                               ) 
                             
                           
                         
                       
                     
                   
                   
                     
                       
                         
                           = 
                           
                             
                               
                                 λ 
                                 H 
                               
                               ⁡ 
                               
                                 ( 
                                 
                                   t 
                                   ❘ 
                                   
                                     H 
                                     t 
                                   
                                 
                                 ) 
                               
                             
                             = 
                             
                               
                                 
                                   p 
                                   * 
                                 
                                 ⁡ 
                                 
                                   ( 
                                   t 
                                   ) 
                                 
                               
                               + 
                               
                                 
                                   
                                     q 
                                     * 
                                   
                                   ⁡ 
                                   
                                     ( 
                                     t 
                                     ) 
                                   
                                 
                                 ⁢ 
                                 
                                   
                                     λ 
                                     PD 
                                   
                                   ⁡ 
                                   
                                     ( 
                                     
                                       t 
                                       ❘ 
                                       
                                         H 
                                         t 
                                       
                                     
                                     ) 
                                   
                                 
                               
                             
                           
                         
                         ; 
                       
                     
                   
                 
               
               
                 
                   Eqn 
                   . 
                   
                       
                   
                   ⁢ 
                   16 
                 
               
             
           
         
       
         
         
           
             where p*(t), q*(t), and u*(t) denote optimal values for p(t), q(t), and u(t). It should thus be noted that there are infinitely many optimal solutions for p*(t), q*(t) that define the stochastic properties of u*(t). Specifically, the set of all possible optimal solutions lie on the line segment  80  of  FIG. 5 . This line segment moves up and down with a slope of −1 as the values of λ H (t|H t ) and λ PD (t|H t ) change over time, with the constraint that λ H (t|H t ) and λ PD (t|H t ) are always between 0 and 1.lf, for example, λ H (t|H t ) is close to zero and λ PD (t|H t ) is close to one, then both p*(t) and q*(t) are close to zero, which implies that u*(t)=−1 with a probability close to 1. Thus, the healthy spike train y H (t) and PD spike train y(t) will respectively be 0 and 1 with high probability and the DBS signal should inhibit y(t) so that z(t)=y(t)+u(t)=0. 
           
         
       
    
     The distance function can alternatively be dependent on the actual behavior of the PD subject in which the DBS system  10  of  FIG. 2  is implanted. Instead of matching the neural activity patterns of the subject to those of a healthy reprehensive, the distance between the motor behavior, for example, movement time, movement velocity, or tremor frequency resulting from the motor behavior, of a healthy and the PD subject can be minimized. For example, this can be achieved by minimizing the following relationship:
 
distance{λ PD ( t|H   t   ,u ( t )), λ H ( t|H   t )}= E{[V ( t ) PD   −V ( t ) H ] 2 }  Eqn. 17;
         where V(t), for example, is the movement velocity as a function of time. Since V(t) is a function of neural activity, a model of the body, as depicted in  FIG. 4 , should be estimated either from measured data or biophysical laws. Further, though single input-output models for PD subject STN activity are described above, other input-output point process models with CIFs λ H (t|H t ,u(t)) can be selected and subjected to similar optimization procedures.       

     Referring again to  FIG. 2 , regardless of the distance function used, the information provided by the controller  20  following minimization can be translated by the pulse generator  18  into a continuous voltage capable of eliciting an appropriate neural activity from the PD subject. As described above, the corrective stimulation signal can take on one of three possible pulse height values, {−1, 0, 1}, at any given time bin. Therefore, a value of u(t)=−1 indicates that the pulse generator  18  should produce an electrical signal that inhibits a target neuron, while a value of u(t)=1 indicates that the pulse generator  18  should produce an electrical signal that excites the target neuron. Likewise, a value of u(t)=0 indicates that the pulse generator  18  should do nothing. 
     Because the present invention administers a stimulation signal to the subject based on observed and predicted physiological parameters of a PD subject rather than a stimulation pattern set by a physician, it can be considered as a “self-calibrating” device. Accordingly, the present can test subject response to a broad set of stimulation signal waveforms and, based on their observed effect, develop waveforms that provide improved patient response and device performance. For example, the present invention could automatically develop a stimulation signal for a particular subject using stimulation signal frequencies lower that those used by traditional DBS devices. This is advantageous, since a lower frequency stimulation signals offers reduced power consumption, thereby prolonging device battery life, and reduced patient side effects by substantially limiting the leakage of stimulation signal to brain areas surrounding the target area. Repeatedly adjusting the stimulation signal to this degree would be prohibitively time-consuming if done using tradition, that is, manual, calibration techniques. 
     The present invention has been described in terms of the various aspects and features, and it should be appreciated that many equivalents, alternatives, variations, and modifications, aside from those expressly stated, are possible and within the scope of the invention. Therefore, the invention should not be limited to a particular described embodiment.