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;
//Note: 131 Ism spikes,
#include <stdlib.h>
#include <time.h>
#include <iostream>
#include <fstream>
#include <math.h>
#include <mpi.h>
#include "control.h"
#include "integrator.h"

using namespace std;

int main( int argc, char** argv )
{

int myRank, poolSize;
MPI_Init( &argc, &argv );
MPI_Comm_size( MPI_COMM_WORLD, &poolSize );
MPI_Comm_rank( MPI_COMM_WORLD, &myRank );

srand( time( NULL ) );

if( myRank == 0 ) { //if I am the master
   MPI_Status statusMaster;
   int i, j, destination, tagSent, tagReceived, counter;
   int ctrlGenome[POPSIZE][NCTRL]={0};
   int geneSent[NCTRL]={0};
   double resultReceived[1]={0};
   double score[POPSIZE]={0};
   double temp=0;

   //input and output files
   ofstream ctrlFile( "process.dat", ios::out );
   ofstream scoreFile( "score.dat", ios::out );
   if( !scoreFile || !ctrlFile ) {
      cerr << "Output files can't be created"; exit(1);
      }

   //genome evolution 
   popInitializer( ctrlGenome ); //genome initialization
   scoreFile << "#Generation & minimum score & mean score" << endl;

   counter=1;
   while( counter<=MAXGEN ) {
      //send out one genome to each slave
      for( destination=1; destination<=POPSIZE; destination++ ) {
         tagSent=destination-1;
         for( i=0; i<NCTRL; i++ ) geneSent[i] = ctrlGenome[tagSent][i];
         MPI_Send( geneSent, NCTRL, MPI_INT, destination, tagSent,
                   MPI_COMM_WORLD );
         }

      //get the result from each slave 
      for( i=1; i<=POPSIZE; i++ ) {
         MPI_Recv( resultReceived, 1, MPI_DOUBLE, MPI_ANY_SOURCE, 
		   MPI_ANY_TAG, MPI_COMM_WORLD, &statusMaster );
         tagReceived = statusMaster.MPI_TAG;
         score[tagReceived]=resultReceived[0];//NOTE: score minimization
         }

      //sort the scores in ascending order
      piksr2( POPSIZE, score, ctrlGenome ); 

      //output the scores
      ctrlFile << "Generation no. " << counter << endl;
      for( i=0; i<POPSIZE; i++ ) {
         ctrlFile << "Control: ";
         for( j=0; j<NCTRL; j++ ) ctrlFile << ctrlGenome[i][j] << " ";
         ctrlFile << endl << "Score = " << score[i] << endl;
         }
      scoreFile << counter << " " << score[0] << " ";//output the best score
      for( i=0; i<POPSIZE; i++ ) temp+=score[i]/POPSIZE;
      scoreFile << temp << endl; //output the mean score
      temp = 0;

      //crossover and mutation
      crossOver( ctrlGenome );
      muTation( ctrlGenome );

      counter++; //proceed to the next generation
      }
   ctrlFile.close();
   scoreFile.close();
   }

else { //if I am not the master
   MPI_Status statusSlave;
   int i, j, tagReceived;
   int geneReceived[NCTRL]={0}; 
   double resultSent[1]={0};

   for( i=1; i<=MAXGEN; i++ ) {
      MPI_Recv( geneReceived, NCTRL, MPI_INT, 0, MPI_ANY_TAG,
                MPI_COMM_WORLD, &statusSlave ); //receive a genome
      tagReceived = statusSlave.MPI_TAG;

      //genome evaluation
      if( integrator( geneReceived, resultSent ) == 0 ) 
         MPI_Send( resultSent, 1, MPI_DOUBLE, 0, tagReceived, MPI_COMM_WORLD ); 
      else {
         cerr << "Integration failed with control" << endl;
	 for( j=0; j<NCTRL; j++ ) cerr << geneReceived[j] << " ";
	 cerr << endl << endl;
         resultSent[0]=0;
         MPI_Send( resultSent, 1, MPI_DOUBLE, 0, tagReceived, MPI_COMM_WORLD ); 
	 }
      }
   }

MPI_Finalize( );

return 0;
}


//random control current generator
void popInitializer( int ctrl[][NCTRL] )
{
int i, j;
                                                                                
for( i=0; i<POPSIZE; i++ ) {
   for( j=0; j<NCTRL-2; j++ ) //set distribution function
      ctrl[i][j] = randGen( CTRL_L, CTRL_H, CTRL_S );
   ctrl[i][NCTRL-2]=randGen( DUR_L, DUR_H, DUR_S ); //set pulse duration
   ctrl[i][NCTRL-1]=randGen( AMP_L, AMP_H, AMP_S ); //set pulse amplitude
   }
}
                                                                                
                                                                                
//random number generator
int randGen( int low, int high, int step )
{
int nstep, rstep;
                                                                                
nstep = int( ( high - low ) / step ) + 1;
                                                                                
rstep = int( double( nstep ) * rand( ) / ( RAND_MAX + 1.0 ) );
                                                                                
return ( low + rstep * step );
}

//sorting
void piksr2( int n, double arr[ ], int brr[ ][ NCTRL ] )
{
int i, j, k, b[ NCTRL ] = { 0 };
double a;
                                                                                
for( j = 1; j < n; j++ ) {
   a = arr[ j ];
   for( k = 0; k < NCTRL; k++ )
      b[ k ] = brr[ j ][ k ];
   i = j - 1;
   while( i >= 0 && arr[ i ] > a ) {
      arr[ i + 1 ] = arr[ i ];
      for( k = 0; k < NCTRL; k++ )
         brr[ i + 1 ][ k ] = brr[ i ][ k ];
      i--;
      }
   arr[ i + 1 ] = a;
   for( k = 0; k < NCTRL; k++ )
      brr[ i + 1 ][ k ] = b[ k ];
   }
}


//crossover
void crossOver( int arr[ ][ NCTRL ] )
{
int a, b, c, d, i, j, temp;
int brr[POPSIZE][NCTRL];
                                                                                
for( i=0; i<POPSIZE; i++ )
   for( j=0; j<NCTRL; j++ )
      brr[i][j]=0;
                                                                                
//crossover
for( i=NSAVED; i<POPSIZE-NMUT; i++ ) {
   a = int( double( NCTRL ) * rand( ) / ( RAND_MAX + 1.0 ) );
   b = int( double( NCTRL ) * rand( ) / ( RAND_MAX + 1.0 ) );
   c = int( double( NUSED ) * rand( ) / ( RAND_MAX + 1.0 ) );
   d = int( double( NUSED ) * rand( ) / ( RAND_MAX + 1.0 ) );
                                                                                
   while( a==b ) a = int( double( NCTRL ) * rand( ) / ( RAND_MAX + 1.0 ) );
   if( a > b ) {
      temp = a;
      a = b;
      b = temp;
      }
   while( c==d ) c = int( double( NUSED ) * rand( ) / (RAND_MAX+1.0) );
                                                                                
   for( j=0; j<NCTRL; j++ ) brr[i][j] = arr[c][j];
   for( j=a; j<b; j++ ) brr[i][j] = arr[d][j];
   }
                                                                                
//put crossover result back to arr[][]
for( i=NSAVED; i<POPSIZE-NMUT; i++ )
   for( j=0; j<NCTRL; j++ )
      arr[i][j]=brr[i][j];
}


//mutation
void muTation( int arr[ ][ NCTRL ] )
{
int a=0, b=0, i, j;
                                                                                
//two point mutation
for( i=POPSIZE-NMUT; i<POPSIZE; i++ ) {
   a = int( double( NUSED ) * rand( ) / ( RAND_MAX + 1.0 ) );
   for( j=0; j<NCTRL; j++ ) arr[i][j]=arr[a][j];
   for( j=0; j<PMUT; j++ ) {
      b = int( double( NCTRL ) * rand( ) / ( RAND_MAX + 1.0 ) );
      if( b==(NCTRL-2) ) arr[i][b]=randGen( DUR_L, DUR_H, DUR_S ); //width
      else if( b==(NCTRL-1) ) arr[i][b]=randGen( AMP_L, AMP_H, AMP_S ); //amp
      else arr[i][b] = randGen( CTRL_L, CTRL_H, CTRL_S ); //PDF
      }
   }
}

Loading data, please wait...