Excitability of DA neurons and their regulation by synaptic input (Morozova et al. 2016a, 2016b)

 Download zip file 
Help downloading and running models
Accession:206380
This code contains conductance-based models of Dopaminergic (DA) and GABAergic neurons, used in Morozova et al 2016 PLOS Computational Biology paper in order to study the type of excitability of the DA neurons and how it is influenced by the intrinsic and synaptic currents. We identified the type of excitability by calculating bifurcation diagrams and F-I curves using XPP file. This model was also used in Morozova et al 2016 J. Neurophysiology paper in order to study the effect of synchronization in GABAergic inputs on the firing dynamics of the DA neuron.
Reference:
1 . Morozova EO, Myroshnychenko M, Zakharov D, di Volo M, Gutkin B, Lapish CC, Kuznetsov A (2016) Contribution of synchronized GABAergic neurons to dopaminergic neuron firing and bursting. J Neurophysiol 116:1900-1923 [PubMed]
Model Information (Click on a link to find other models with that property)
Model Type: Neuron or other electrically excitable cell;
Brain Region(s)/Organism:
Cell Type(s): Substantia nigra pars compacta DA cell; Ventral tegmental area dopamine neuron; Ventral tegmental area GABA neuron ;
Channel(s): I K,Ca; I Calcium; I Na,t; I Potassium;
Gap Junctions:
Receptor(s): GabaA; NMDA; AMPA;
Gene(s):
Transmitter(s): Gaba; Glutamate; Dopamine;
Simulation Environment: MATLAB; XPP;
Model Concept(s): Action Potentials; Bifurcation; Bursting; Synaptic Convergence;
Implementer(s): Morozova, Ekaterina O [emorozov at indiana.edu]; Kuznetsov, Alexey ;
Search NeuronDB for information about:  Substantia nigra pars compacta DA cell; GabaA; AMPA; NMDA; I Na,t; I K,Ca; I Calcium; I Potassium; Dopamine; Gaba; Glutamate;
#include <stdio.h>
#include <stdlib.h>
#include "mex.h"    
 #include "matrix.h"   
 #include <string.h> 
#include <math.h>
#include <time.h>
 #include <iostream>

//-------------parameters --------------------------
double gL=0.18, gbarK=1, gbarCa=2.5, gbarKCa=7.8, gSNa=0.13; // maximal conductances for pacemaking
double gbarDR=0, gbarNa=0, gNa; // maximal conductances for spike producing current
// set gbarDR=2 and gbarNa=50 to add spikeproducing currents
double gDR, gCa, gKCa, gK; //conductances
double gampa, gnmda, gbarnmda, ggaba; //synaptic conductances
double nmdasig, ampasig, allgaba; //synaptic activation functions
double k=160, buf=0.00023, zF=0.019298, tc=0.52, r=0.2; //for Ca equation
double Mg=0.5, me=0.062; //for gnmda 
double EL=-35, EK=-90, ECa=50, ENMDA=0, EAMPA=0, EGABA=-90, ENa=55; // reversal potentials
double VHK=-10, VSK=7; //K current
double th=0.05; //for Na curent
double Vcah=-52, Sca=3; // for Ca current
double vhna=-50, slna=5; // for subthreshold Na current
double aact=1., tades=6.1, adesrel=40, adeact=1.6, nact=7., ndeact=170.; //for AMPA 
double nmdathresh=0.2, nmdasl=0.03, ampathresh=0.2, ampasl=0.03; //for NMDA
double dt=0.02, TT, tgl; //step, total time, global time
int N, Neq=9; // number of steps, number of equations
double inp; // input
double caleak=0.1; // fraction of Ca current in a leak current
double Iapp; //applied current
double vhh=-95, slh=8, th0=625, vtauh=-112, Eh=-20, gh, gbarh; // for Ih current
double *y, *f;
//---------------------------------------------------------
void system (void)
{    
// Na current
double alpham=-0.32*(y[0]+39)/(exp(-(y[0]+39)/4)-1);
double betam=0.28*(y[0]+4)/(exp((y[0]+4)/5)-1);
double minf=alpham/(alpham+betam);
double alphah=0.2*th*exp(-((y[0]+47.)/18.));
double betah=25.*th/(1.+(exp(-(y[0]+24.)/5.)));
gNa=gbarNa*(pow(minf,3))*y[5];
// Delayed rectifier K current
double alphan=-0.0032*(y[0]+5.)/(exp(-(y[0]+5.)/10.) - 1.);
double betan=0.05*exp(-((y[0]+10.)/16.));
double nsqr=y[4]*y[4];
gDR=gbarDR*(nsqr*nsqr);
// for K current
gK=gbarK/(1. + exp(-(y[0]-VHK)/VSK));
// for Ca current
double alphac=((fabs(y[0]-Vcah))>0.00001)? (-0.0032*(y[0]-Vcah)/(exp(-(y[0]-Vcah)/Sca) - 1.)) : (-0.0032*0.00001/(exp(-0.00001/Sca)-1.));
double betac=0.05*exp(-(y[0]-Vcah+5)/40.);
double csinf=alphac/(alphac+betac);
double csqr=csinf*csinf;
gCa=gbarCa*(csqr*csqr);
// for Ca-dependent K current
double ksq=k*k;
double casq=y[3]*y[3];
gKCa=gbarKCa*(casq*casq)/((casq*casq) + (ksq*ksq));
// for NMDAR current
nmdasig=1/(1+exp(-(y[1]-nmdathresh)/nmdasl));
gnmda=gbarnmda/(1+0.28*Mg*exp(-me*y[0]));
// for AMPAR current
ampasig=1/(1+exp(-(y[2]*y[6]-ampathresh)/ampasl));
// for subthreshold Na current
double na=1/(1+exp(-(y[0]-vhna)/slna));
// for H-current
double hinf=1/(1+exp((y[0]-vhh)/slh));
double tauh=th0*exp(0.075*(y[0]-vtauh))/(1+exp(0.083*(y[0]-vtauh)));
gh=gbarh*y[7];

f[0]= gnmda*(ENMDA-y[0])+ gampa*(EAMPA-y[0])+ ggaba*(EGABA-y[0]) + gCa*(ECa-y[0])
+ (gKCa+gK+gDR)*(EK-y[0])+ gL*(EL-y[0])+ gNa*(ENa-y[0]) + gSNa*na*(ENa-y[0])
+ gh*(Eh-y[0])+Iapp;// Voltage
f[1]=inp*(1-y[1])/nact-(1-inp)*y[1]/ndeact; // NMDA
f[2]=inp*(1-y[2])/aact-(1-inp)*y[2]/adeact; // AMPA activtion
f[3]= 2.*buf*((gCa+caleak*gL)*(ECa - y[0])/zF - y[3]/tc)/r; //Ca
f[4]= alphan*(1.-y[4])-betan*y[4]; // DR
f[5]= alphah*(1.-y[5])-betah*y[5]; // Na
f[6]=(1-inp)*(1-y[6])/adesrel-inp*y[6]/tades; //AMPA desensitization
f[7]= (hinf-y[7])/tauh; //Ih
}
void euler(double ggaba_1, double gbarnmda_1, double gampa_1)
{
 int i=0;
 ggaba=ggaba_1;
 gbarnmda=gbarnmda_1;
 gampa=gampa_1;
 system();
  while (i<Neq){y[i]+=dt*f[i]; i++;}
} 

//nlhs - Number of expected output mxArrays
//plhs - Array of pointers to the expected output mxArrays
//nrhs - Number of input mxArrays
//prhs - Array of pointers to the input mxArrays. Do not modify any prhs values in your MEX-file. Changing the data in these read-only mxArrays can produce undesired side effects.

void mexFunction( int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[] )
{
mxArray *tempy = mxCreateDoubleMatrix(1,Neq,mxREAL);
y=mxGetPr(tempy);
mxArray *tempf = mxCreateDoubleMatrix(1,Neq,mxREAL);
f=mxGetPr(tempf);
// inputs
double *ggaba1, *gbarnmda1,*gampa1, *Iapp1;      
TT   = *(mxGetPr(prhs[0]));
ggaba1 = (mxGetPr(prhs[1]));
gbarnmda1 = (mxGetPr(prhs[2]));
gampa1 = (mxGetPr(prhs[3]));
gbarh = *(mxGetPr(prhs[4]));
Iapp1 = (mxGetPr(prhs[5]));

int kk;
tgl=0;
dt=0.02;
N=(int)(TT/dt);
// outputs    
double *Vm;  
double *Ca;
plhs[0] = mxCreateDoubleMatrix(1,N,mxREAL);
plhs[1] = mxCreateDoubleMatrix(1,N,mxREAL);
Vm = mxGetPr(plhs[0]);
Ca = mxGetPr(plhs[1]);

y[0]=-60; y[1]=0.; y[2]=0; y[3]=40; y[4]=0; y[5]=0; y[6]=1.; y[7]=0.06; // initial conditions
	
kk=0;
for(int i=0; i<N; i++)
{
tgl+=dt;
inp=0;
ggaba=ggaba1[i];
Iapp=Iapp1[i];
gbarnmda=gbarnmda1[i];
gampa=gampa1[i];
euler(ggaba,gbarnmda,gampa);
Vm[kk]=y[0];
Ca[kk]=y[3];

if(kk>N) break;
kk++;
   }
   return;
} // end mexFunction()

Loading data, please wait...