Optimal deep brain stimulation of the subthalamic nucleus-a computational study (Feng et al. 2007)

 Download zip file 
Help downloading and running models
Accession:93449
Here, we use a biophysically-based model of spiking cells in the basal ganglia (Terman et al., Journal of Neuroscience, 22, 2963-2976, 2002; Rubin and Terman, Journal of Computational Neuroscience, 16, 211-235, 2004) to provide computational evidence that alternative temporal patterns of DBS inputs might be equally effective as the standard high-frequency waveforms, but require lower amplitudes. Within this model, DBS performance is assessed in two ways. First, we determine the extent to which DBS causes Gpi (globus pallidus pars interna) synaptic outputs, which are burstlike and synchronized in the unstimulated Parkinsonian state, to cease their pathological modulation of simulated thalamocortical cells. Second, we evaluate how DBS affects the GPi cells' auto- and cross-correlograms.
References:
1 . Terman D, Rubin JE, Yew AC, Wilson CJ (2002) Activity patterns in a model for the subthalamopallidal network of the basal ganglia. J Neurosci 22:2963-76 [PubMed]
2 . Rubin JE, Terman D (2004) High frequency stimulation of the subthalamic nucleus eliminates pathological thalamic rhythmicity in a computational model. J Comput Neurosci 16:211-35 [PubMed]
3 . Feng XJ, Shea-Brown E, Greenwald B, Kosut R, Rabitz H (2007) Optimal deep brain stimulation of the subthalamic nucleus--a computational study. J Comput Neurosci 23:265-82 [PubMed]
Model Information (Click on a link to find other models with that property)
Model Type: Realistic Network;
Brain Region(s)/Organism: Basal ganglia;
Cell Type(s): Globus pallidus neuron;
Channel(s): I T low threshold; I Sodium; I Potassium;
Gap Junctions:
Receptor(s): Glutamate; Gaba;
Gene(s):
Transmitter(s): Gaba; Glutamate;
Simulation Environment: C or C++ program;
Model Concept(s): Parkinson's; Deep brain stimulation;
Implementer(s): Feng, Xiao-Jiang [xfeng at mahler.princeton.edu];
Search NeuronDB for information about:  Glutamate; Gaba; I T low threshold; I Sodium; I Potassium; Gaba; Glutamate;
//general parameters 
#define PI 3.1415926
#define NCELL 8 //number STN, GPe, or GPi cells
#define STNGPE NCELL*2 //sum of STN and GPe cells
#define NEQN 6 //number of ODEs per cell
#define NTCCELL 2 //number of TC cells
#define NTCEQN 3 //number of ODEs per TC cell
#define P_cm 1.0

//STN parameters
#define P_phin 0.75
#define P_phih 0.75
#define P_phica 0.75 //0.75 in Terman's code, not available from his paper.
#define P_eps 3.75e-5 //5.0e-5 in Terman's code
#define P_kca 22.5
#define P_phir 0.2
#define P_alpha 1.0 
#define P_thetag 30.0
#define P_beta 0.05 
#define P_gGS 0.9 
#define P_vGS -100.0 
#define P_gl 2.25
#define P_vl -60.0
#define P_gk 45.0
#define P_vk -80.0
#define P_gna 37.5
#define P_vna 55.0
#define P_gt 0.5
#define P_vca 140.0
#define P_gca 0.5
#define P_gahp 9.0
#define P_k1 15.0
#define P_thetam -30.0
#define P_sigmam 15.0
#define P_thetaa -63.0
#define P_sigmaa 7.8
#define P_thetas -39.0
#define P_sigmas 8.0
#define P_thetan -32.0
#define P_sigman 8.0
#define P_thetah -39.0
#define P_sigmah -3.1
#define P_thetar -67.0
#define P_sigmar -2.0
#define P_thetab 0.4
#define P_sigmab -0.1
#define P_thetaH -39.0
#define P_sigmaH 8.0
#define P_taun0 1.0
#define P_taun1 100.0
#define P_thn -80.0
#define P_sigmant -26.0
#define P_tauh0 1.0
#define P_tauh1 500.0
#define P_thh -57.0
#define P_sigmaht -3.0
#define P_taur0 40.0
#define P_taur1 17.5
#define P_thr 68.0
#define P_sigmart -2.2
#define ISTN 25.0 

//GPe parameters
#define P_phing 0.1
#define P_phihg 0.05
#define P_phirg 1.0
#define P_alphag 1.0 
#define P_thetagg 20.0
#define P_betag 0.1 
#define P_epsg 1.0e-4 
#define P_kcag 15.0 
#define P_gSG 0.6 //0.3 in Terman 2004
#define P_vSG 0.0 
#define P_gGG 0.0005 //normal = 2.0, PD = 0.05 or 0.0005
#define P_vGG -80.0 
#define P_glg 0.1
#define P_vlg -55.0
#define P_gkg 30.0
#define P_vkg -80.0
#define P_gnag 120.0
#define P_vnag 55.0
#define P_gtg 0.5
#define P_vcag 120.0
#define P_gcag 0.15
#define P_gahpg 30.0
#define P_k1g 30.0
#define P_thetamg -37.0
#define P_sigmamg 10.0
#define P_thetaag -57.0
#define P_sigmaag 2.0
#define P_thetasg -35.0
#define P_sigmasg 2.0
#define P_thetang -50.0
#define P_sigmang 14.0
#define P_thetahg -58.0
#define P_sigmahg -12.0
#define P_thetarg -70.0
#define P_sigmarg -2.0
#define P_thetaHg -57.0
#define P_sigmaHg 2.0
#define P_taun0g 0.05
#define P_taun1g 0.27
#define P_thng -40.0
#define P_sigmantg -12.0
#define P_tauh0g 0.05
#define P_tauh1g 0.27
#define P_thhg -40.0
#define P_sigmahtg -12.0
#define P_taurg 30.0
#define IGPE -13.3 //normal = -2.0, PD = -13.0

//GPi parameters
#define P_phini 0.1
#define P_phihi 0.05
#define P_phiri 1.0
#define P_alphai 2.0
#define P_thetagi 20.0
#define P_betai 0.08
#define P_epsi 1.0e-4
#define P_kcai 15.0
#define P_gSI 0.2 //0.3 in Terman 2004
#define P_vSI 0.0 
#define P_gGI 1.7 //1.0 in Terman 2004
#define P_vGI -100.0 
#define P_gli 0.1
#define P_vli -55.0
#define P_gki 30.0
#define P_vki -80.0
#define P_gnai 120.0
#define P_vnai 55.0
#define P_gti 0.5
#define P_vcai 120.0
#define P_gcai 0.15
#define P_gahpi 30.0
#define P_k1i 30.0
#define P_thetami -37.0
#define P_sigmami 10.0
#define P_thetaai -57.0
#define P_sigmaai 2.0
#define P_thetasi -35.0
#define P_sigmasi 2.0
#define P_thetani -50.0
#define P_sigmani 14.0
#define P_thetahi -58.0
#define P_sigmahi -12.0
#define P_thetari -70.0
#define P_sigmari -2.0
#define P_thetaHi -57.0
#define P_sigmaHi 2.0
#define P_taun0i 0.05
#define P_taun1i 0.27
#define P_thni -40.0
#define P_sigmanti -12.0
#define P_tauh0i 0.05
#define P_tauh1i 0.27
#define P_thhi -40.0
#define P_sigmahti -12.0
#define P_tauri 30.0
#define IGPI 3.0 

//TC parameters (obtained from Terman 2004)
#define P_gIT 0.08 //0.06 in Terman 2004
#define P_vIT -85.0
#define P_glt 0.05
#define P_vlt -70.0
#define P_gkt 5.0
#define P_vkt -90.0
#define P_gnat 3.0
#define P_vnat 50.0
#define P_gtt 5.0
#define P_vcat 0.0
#define P_thetamt -37.0
#define P_sigmamt 7.0
#define P_thetaat -60.0
#define P_sigmaat 6.2
#define P_thetaht -41.0
#define P_sigmahtt -4.0
#define P_thetart -84.0
#define P_sigmartt -4.0
#define ISM 10.0 //5.0 in Terman 2004, 7.5-10 may be good
#define PSM 50.0 //25.0 in Terman 2004, 50 in his figures
#define DSM 5.0
#define SHIFT 2.5 //phase shift when counting spike number
#define THLD_TC -30.0 // threshold for TC voltage spikes

//integrator parameters
#define NEQ (NCELL*NEQN*3+NTCCELL*NTCEQN) //number of ODEs in total
#define LRW ( 22 + NEQ * ( NEQ + 9 ) )
#define LIW ( 20 + NEQ )
#define ITOL 2
#define RTOL 1.0e-10 
#define ATOL 1.0e-12 
#define ATOLV 1.0e-10 //absolute tolerance for v
#define ITASK 1
#define ISTATE 1
#define IOPT 0
#define JT 2
#define TOUT 0.1 //first output time
#define TSTEP 0.1 //output step size
#define NSTEP 65536
#define NDELAY 25000 //No. of steps before the control is applied

//current input parameters
#define ILOW 1.0 //lowest current pulse distance
#define ISTEP 5.0 //step size for increasing pulse distance
#define DEV_DBS 0.1 //white noise of I_{DBS}
#define DEV_NET 0.01 //heterogeneity in the network itself (0 to 1)

//parameters for counting spike time
#define NSPK 2000 //maximum no. of Gpi spikes
#define VTHLD 0.5 //Gpi synaptic output threshold for counting spikes
                                                                                
//parameters for calculating auto- and cross-correlations
#define NSTART 10000 //when to start calculating correlations (in steps)
#define NEND 55000 //when to end calculating correlations (in steps)
#define DURATION 8192 //duration for correlation calculation (in steps)
#define BIN_RESLN 100 //binning resolution of Gpi spikes (in steps)
#define NBIN DURATION/BIN_RESLN+1 //maximum no. of Gpi bins
#define COMBI NCELL-1 //no. of comparisons for crosscorrelation
                                                                                
//parameters for convolution
#define DURATION2 DURATION*2
#define MSTEP 21 //width of the response function (in bins)
#define WIDTH 1 //Gaussian width (in bins)

//parameters for cost function calculation
#define W1 200.0 //1.0 //weight for autocorrelation term
#define W2 100.0 //weight for crosscorrelation term
#define W3 0.0002 //weight for Gpi spike number in first bin term
#define W4 0 //0.1 //current cost weight

//function declaration

int integrator( int [], double [] ); //my integrator

extern "C" void lsoda_( void( int *, double *, double *, double * ), int *,
double *, double *, double *, int *, double *, double *, int *, int *, int *,
double *, int *, int *, int *, void( int *, double *, double *, int *, int *,
double *, int * ), int * ); //ODE integrator

void fex( int *, double *, double *, double * ); //ODE container

void jdum( int *, double *, double *, int *, int *, double *, int * ); //dummy

//STN currents
double Il( const double v );
double Ik( const double v, const double n );
double Ina( const double v, const double h );
double It( const double v, const double r );
double Ica( const double v );
double Iahp( const double v, const double ca );

//STN functions
double minf( const double v );
double ainf( const double v );
double sinf( const double v );
double ninf( const double v );
double hinf( const double v );
double rinf( const double v );
double binf( const double r );
double Hinf( const double v );
double taun( const double v );
double tauh( const double v );
double taur( const double v );

//GPe currents
double Ilg( const double v );
double Ikg( const double v, const double n );
double Inag( const double v, const double h );
double Itg( const double v, const double r );
double Icag( const double v );
double Iahpg( const double v, const double ca );

//GPe functions
double minfg( const double v );
double ainfg( const double v );
double sinfg( const double v );
double ninfg( const double v );
double hinfg( const double v );
double rinfg( const double v );
double Hinfg( const double v );
double taung( const double v );
double tauhg( const double v );

//GPi currents
double Ili( const double v );
double Iki( const double v, const double n );
double Inai( const double v, const double h );
double Iti( const double v, const double r );
double Icai( const double v );
double Iahpi( const double v, const double ca );

//GPi functions
double minfi( const double v );
double ainfi( const double v );
double sinfi( const double v );
double ninfi( const double v );
double hinfi( const double v );
double rinfi( const double v );
double Hinfi( const double v );
double tauni( const double v );
double tauhi( const double v );

//TC currents
double Ilt( const double v );
double Ikt( const double v, const double h );
double Inat( const double v, const double h );
double Itt( const double v, const double r );
double taurt( const double v );
double tauht( const double v );

//TC functions
double minft( const double v );
double ainft( const double v );
double hinft( const double v );
double rinft( const double v );
double Heaviside( const double var );

//other functions
//inhibition of GPe to STN
double IsynGS(const double, int, const double [ ]); 
//neighboring inhibition of GPe to GPe
double IsynGG(const double, int, const double [ ] ); 
//synaptic current from GPi to TC
double IsynIT1( const double [], const double );
double IsynIT2( const double [], const double );
//DBS current to STN
double Idbs( const unsigned long int, const int ); 
//cost value calculation
//Ism current to TC
double Ism( const double );
//Ism current to TC with delay, used in peak counting
double Ism2( const double );
//heterogeneity / random number generator
double hetero( const int, const double );
//random square pulse DBS input generated from specified PDF
void randInput( const double [], double [] );
//random number generator for specified probability distribution function
double randGen2( const double [], const int, const double, const double );
//auto- and cross-correlation calculation
double costCal( const double [][NSTEP] );
//count the number of spikes
int spkCount ( const double [], double [] );
//calculate autocorrelation
void autoCal( const double [], int [] );
//convolution calculation
void convCal( const int [], float [] );
//calculate standard deviation
double stdCal ( const float[] );
                                                                                
void four1( float [], unsigned long int, int ); //fft
void realft( float [], unsigned long, int );
void twofft( float [], float [], float [], float [], unsigned long );
void convlv( float [], unsigned long, float [], unsigned long, int, float [] );

Loading data, please wait...