Cholinergic and nicotinic regulation of DA neuron firing (Morozova et al 2020)

 Download zip file 
Help downloading and running models
The model describes the modulation of firing properties of DA neurons by acetylcholine (ACh) and nicotine in 5 cases: knock-out of ß2-containing nAChRs, ß2-containing nAChRs only on DA neurons, the nAChRs only on GABA neurons, the nAChRs on both DA and GABA neurons and “wild” type (the AChRs on DA, GABA and Glu neurons). The distinct responses to ACh and nicotine could be explained by distinct temporal patterns of these inputs: pulsatile vs. continuous.
1 . Morozova E, Faure P, Gutkin B, Lapish C, Kuznetsov A (2020) Distinct temporal structure of nicotinic ACh receptor activation determines responses of VTA neurons to endogenous ACh and nicotine eNeuro, accepted
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): Ventral tegmental area dopamine neuron; Ventral tegmental area GABA neuron ;
Channel(s): I K,Ca; I Calcium; I K;
Gap Junctions:
Receptor(s): Cholinergic Receptors;
Transmitter(s): Acetylcholine; Gaba; Dopamine; Glutamate;
Simulation Environment: MATLAB;
Model Concept(s): Bursting;
Implementer(s): Morozova, Ekaterina O [emorozov at];
Search NeuronDB for information about:  Cholinergic Receptors; I K; I K,Ca; I Calcium; Acetylcholine; Dopamine; Gaba; Glutamate;
#include <stdio.h>
#include <stdlib.h>
#include "mex.h" // Standard include library for the MEX-files.
 #include "matrix.h" // Standard include library for accessing MATLAB data struct.
 #include <string.h> 
#include <math.h>
#include <time.h>
 #include <iostream>

//---------- parameters for 1-compartmental model------------
double Cgaba=10;
double gL, gbarK=1/2, gbarCa=2.5/2, gbarKCa=7.8/2, gbarDR=0*10, gbarERG=0*4.8, gNa=0,gSNa=0.13; //conductance amplitudes
double gDR, gERG, 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 ergtaumax=300, ergtaumin=62, Vnh=-50.4, Sn=1.8, sh=13; //for ERG
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; //for gK
double th=0.05; //for gNa
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; // old for NMDA
double nmdathresh=9, nmdasl=1.3, ampathresh=9, ampasl=1.3; //for NMDA
double TT, dt=0.02, tgl; // time, step, global time
int N, NG=20, Neq=15+NG*4; // N of steps,  N of GABA neurons, N of equations,
double amg, bmg, minfgg, ang, bng, ahg, bhg, gnmdagg, gspikeg, vgnz; // for GABA population
double glg=0.1*10, gna=22*10, gdrg=7*10, tng=1/2.4, thg=5, tbn=0.7, as=12, bs=0.1; // for GABA population
double gmf, dg=0; //strength of Gap junction
double gampag=0, gnmdabarg=0; //GABA synaptic conductances
double inp, inpachda, inpachgaba, inpnic; // input
double vhna=-50, slna=5; // for subthresh Na current
double Vcah=-52, Sca=3; // for Ca current
double caleak=0.1; // fraction of Ca current in a leak current
double vhh=-95, slh=8, th0=625, vtauh=-112, Eh=-20, gh, gbarh; // Ih (Okamoto)
//--------------for Ach a2b2 receptor from paper---------------
double EC50_actB2=30, Hill_actB2=1.05, w_actNicB2_da=3, w_actNicB2_gaba=3;
double EC50_desB2=0.061; // half-maximum concentration of desensitization by Nic (microM)
double Hill_desB2=0.5; // Hill coefficient of desensitization
double tmin_desB2=500; // minimal des. time constant (500 msec)
double tmax_desB2=600*10^3; // maximal des. time constant (10 min)
double Ktau_desB2=0.11; // half maximum conc of desensitization time constant (microM)
double Htau_desB2=3; // Hill coefficient of des. time constant
double gachgaba, gachda, gbarachgaba, gbarachda;
double EAch=0; // reversal potential
double achtau, achinfda, achinfgaba;
double achtaudes, achinfdes;

double *y;
double *f;
double *rngL=new double[NG];
double *glg1=new double[NG];

void system (void)
double alphan=-0.0032*(y[0]+5.)/(exp(-(y[0]+5.)/10.) - 1.);
double betan=0.05*exp(-((y[0]+10.)/16.));
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 alphah=0.2*th*exp(-((y[0]+47.)/18.));
double betah=25.*th/(1.+(exp(-(y[0]+24.)/5.)));
double ergtau=ergtaumin + ergtaumax*(1/(1 + exp((y[0]-Vnh)/Sn))-1/(1+exp((y[0]-Vnh+sh)/Sn)));
double erginf=1/(1+exp(-(y[0]-Vnh-3)/Sn));
double nsqr=y[5]*y[5];
double minf=alpham/(alpham+betam);
double ergsqr=y[4]*y[4];
gK=gbarK/(1. + exp(-(y[0]-VHK)/VSK));
//double nmdasig=1/(1+exp(-(y[1]-nmdathresh)/nmdasl));
//double ampasig=1/(1+exp(-(y[2]*y[6]-ampathresh)/ampasl));
double nmdasig=1/(1+exp(-(inp-nmdathresh)/nmdasl));
double ampasig=1/(1+exp(-(inp-ampathresh)/ampasl));
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;
double ksq=k*k;
double casq=y[3]*y[3];
gKCa=gbarKCa*(casq*casq)/((casq*casq) + (ksq*ksq));
double na=1/(1+exp(-(y[0]-vhna)/slna)); // subthresh Na
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)));
//double tauh=550+1860*exp(-(pow((y[0]-(-83.6))/19.4,2)));
gachda = gbarachda*y[8]*(1-y[10]);
gachgaba = gbarachgaba*y[9]*(1-y[10]);
achtau=1; // 5 ms // was 1 for some reasone ( everything worked), changed back to 5
achinfdes = 1/(1+pow(EC50_desB2/inpnic,Hill_desB2));
achtaudes = tmin_desB2 + tmax_desB2 * 1/(1+pow(inpnic/Ktau_desB2,Htau_desB2));
achinfda = 1/(1+pow(EC50_actB2/(inpachda+w_actNicB2_da*inpnic),Hill_actB2));
achinfgaba = 1/(1+pow(EC50_actB2/(inpachgaba+w_actNicB2_gaba*inpnic),Hill_actB2));

for (int ig=0; ig<NG; ig++){

f[0]= gnmda*y[1]*(ENMDA-y[0])+ gampa*y[2]*y[6]*(EAMPA-y[0])+ ggaba*allgaba*(EGABA-y[0]) + gCa*(ECa-y[0])
+ (gKCa+gK+gERG/*+gDR*/)*(EK-y[0])+ gL*(EL-y[0])+ gNa*(ENa-y[0]) + gSNa*na*(ENa-y[0]) + gachda*(EAch-y[0]) + gh*(Eh-y[0]);// + ghNa*(ENa-y[0]) + ghK*(EK-y[0]);
f[1]=nmdasig*(1-y[1])/nact-(1-nmdasig)*y[1]/ndeact; // NMDA
f[2]=ampasig*(1-y[2])/aact-(1-ampasig)*y[2]/adeact; // AMPA
f[3]= 2.*buf*((gCa+caleak*gL)*(ECa - y[0])/zF - y[3]/tc)/r; //Ca
f[4]= (erginf-y[4])/ergtau; //ERG
f[5]= alphan*(1.-y[5])-betan*y[5]; // DR
f[6]=(1-ampasig)*(1-y[6])/adesrel-ampasig*y[6]/tades; //ampa desensitization
f[7]= alphah*(1.-y[7])-betah*y[7]; // Na
f[8]= (achinfda-y[8])/(achtau); //Ach
f[9]= (achinfgaba-y[9])/achtau; //Ach
f[10]= (achinfdes-y[10])/achtaudes; //Ach
f[11]= (hinf-y[11])/tauh; //Ih

// population of GABA neurons
for (int ig=0; ig<NG; ig++){gmf+=y[12+ig*4];} //for gap junctions

for (int ig=0; ig<NG; ig++)
double amg=0.1*(y[12+ig*4]+30.0)/(1.0-exp(-(y[12+ig*4]+30.0)/10.0));
double bmg=4.0*exp(-(y[12+ig*4]+55.0)/18.0);
double minfgg=amg/(amg+bmg);
double ahg=0.07*exp(-(y[12+ig*4]+53.0)/20.0);
double bhg=1.0/(1.0+exp(-(y[12+ig*4]+23.0)/10.0));
double ang=0.01*(y[12+ig*4]+29.0)/(1.0-exp(-(y[12+ig*4]+29.0)/10.0));
double bng=tbn*0.125*exp(-(y[12+ig*4]+39.0)/80.0);
double gnmdagg=gnmdabarg/(1+0.28*Mg*exp(-me*y[12+ig*4]));
double gspikeg=1/(1+exp(-y[12+ig*4]/2));

f[12+ig*4]=(gnmdagg*nmdasig*(ENMDA-y[12+ig*4]) + gampag*ampasig*(EAMPA-y[12+ig*4])
    -gdrg*pow(y[13+ig*4],4)*(y[12+ig*4]+90) + dg*(gmf-y[12+ig*4]) + gachgaba*(EAch-y[12+ig*4]))*(1/Cgaba);
f[15+ig*4]=as*gspikeg*(1-y[15+ig*4])-bs*(1-gspikeg)*y[15+ig*4]; // GABA receptor activation;

void euler()
 int i=0;
  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[] )
	srand(1); //edit!   
//if(y) delete[] y;        
// y=new double[Neq];
mxArray *tempy = mxCreateDoubleMatrix(1,Neq,mxREAL);
//if(f) delete[] f;        
// f=new double[Neq];
mxArray *tempf = mxCreateDoubleMatrix(1,Neq,mxREAL);

double *input;
double *ggaba1, *gbarnmda1, *gampa1;

TT   = *(mxGetPr(prhs[0]));
input = mxGetPr(prhs[1]);
ggaba1 = (mxGetPr(prhs[2]));
gbarnmda1 = (mxGetPr(prhs[3]));
gampa1 = (mxGetPr(prhs[4]));

double *inputachda;
double *inputachgaba;
double *inputnic;
inputachda = mxGetPr(prhs[5]);
inputachgaba = mxGetPr(prhs[6]);
inputnic = mxGetPr(prhs[7]);
gbarachgaba = *(mxGetPr(prhs[8]));
gbarachda = *(mxGetPr(prhs[9]));
gL = *(mxGetPr(prhs[10]));
gbarh = *(mxGetPr(prhs[11]));

 int kk;
double *Vm;   
double *mas_Vgaba=new double[NG*N];
double *allgaba1;
double *Achcurrent;
double *nmdasig1;

plhs[0] = mxCreateDoubleMatrix(1,N,mxREAL);
plhs[1] = mxCreateDoubleMatrix(1,N,mxREAL);
plhs[2] = mxCreateDoubleMatrix(1,N,mxREAL);
plhs[3] = mxCreateDoubleMatrix(1,N,mxREAL);
plhs[4] = mxCreateDoubleMatrix(1,NG*N,mxREAL);

    Vm = mxGetPr(plhs[0]);
   mas_Vgaba = mxGetPr(plhs[4]);
   allgaba1 = mxGetPr(plhs[1]);
   Achcurrent = mxGetPr(plhs[2]);
   nmdasig1 = mxGetPr(plhs[3]);

	 for(int i=0; i<NG;i++)
 { for(int j=0;j<N;j++) mas_Vgaba[i*N+j] = 0; }

    y[0]=-60.; y[1]=0.; y[2]=0; y[3]=50; y[4]=0.1; y[5]=0; y[6]=1.;y[7]=0; y[8]=0; y[9]=0; y[10]=0; y[11]=0.06;
 for (int ig=0; ig<NG; ig++){rngL[ig] = (double)rand()/RAND_MAX; y[12+ig*4]=-40+100*(rngL[ig]-0.5); glg1[ig]=(0.1+(0.1*(rngL[ig]-0.5)))*10; y[13+ig*4]=0; y[14+ig*4]=0; y[15+ig*4]=0.;}
    for(int i=0; i<N; i++)
      // inp=0;
       inp = input[i];
       Achcurrent[i]=gachda; //*(EAch-y[8]);
    for(int ig=0;ig<NG;ig++)

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

Loading data, please wait...