Ionic mechanisms of bursting in CA3 pyramidal neurons (Xu and Clancy 2008)

 Download zip file 
Help downloading and running models
Accession:114337
"... We present a single-compartment model of a CA3 hippocampal pyramidal neuron based on recent experimental data. We then use the model to determine the roles of primary depolarizing currents in burst generation. The single compartment model incorporates accurate representations of sodium (Na+) channels (NaV1.1) and T-type calcium (Ca2+) channel subtypes (CaV3.1, CaV3.2, and CaV3.3). Our simulations predict the importance of Na+ and T-type Ca2+ channels in hippocampal pyramidal cell bursting and reveal the distinct contribution of each subtype to burst morphology. We also performed fastslow analysis in a reduced comparable model, which shows that our model burst is generated as a result of the interaction of two slow variables, the T-type Ca2+ channel activation gate and the Ca2+-dependent potassium (K+) channel activation gate. The model reproduces a range of experimentally observed phenomena including afterdepolarizing potentials, spike widening at the end of the burst, and rebound. Finally, we use the model to simulate the effects of two epilepsy-linked mutations: R1648H in NaV1.1 and C456S in CaV3.2, both of which result in increased cellular excitability."
Reference:
1 . Xu J, Clancy CE (2008) Ionic mechanisms of endogenous bursting in CA3 hippocampal pyramidal neurons: a model study. PLoS One 3:e2056 [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): Hippocampus CA3 pyramidal GLU cell;
Channel(s): I Na,t; I L high threshold; I N; I T low threshold; I A; I K; I M; I K,Ca; I Calcium;
Gap Junctions:
Receptor(s): GabaB;
Gene(s): Nav1.1 SCN1A; Cav3.2 CACNA1H; Cav3.1 CACNA1G; Cav3.3 CACNA1I;
Transmitter(s):
Simulation Environment: C or C++ program;
Model Concept(s): Bursting; Ion Channel Kinetics;
Implementer(s):
Search NeuronDB for information about:  Hippocampus CA3 pyramidal GLU cell; GabaB; I Na,t; I L high threshold; I N; I T low threshold; I A; I K; I M; I K,Ca; I Calcium;
/*  Cell model: CA3 hippocampal neuron
 *
 *  Created by jun xu @ Clancy Lab of Cornell University Medical College on 3/27/05
 *  
 *  Geometry: single-compartment model modified on 04/19/07 */ 

#include <iostream>
#include <cmath>
#include <cstdlib>
#include <cstdio>
#include <iomanip>
#include "nr.h"
using namespace std;
#define nseg 1 

	/* Creation of Data File */
	    FILE *out_one;
	    void data_file ();	
	    int printdata;	// printing counter
	    int printval;  // intervals to output results	

// define constant parameters
        double F = 96520; // Farady constant (coulomb/mol)
        double R = 8.3134; // gas constant (J/K.mol)
        double PI = 3.14;
        double celsius = 30;
        double T = celsius+273.14;  // absolute temperature (K)

// define membrane properties
        double Rm = 60000; // 
        double cm = 1.0*1.0e-3; //units
        double Ra = 200*1.0e-4; //units
        double vrest = -65; // (mv)
        double vinit = vrest;
        double vk = -91;
        double vna = 50;
        

// define initial ionic concentration
        double cai_init = 5.0e-5; // (mM)
        double cao = 2; // (mM)

//	cout<<"ok"<<endl;

	/* Voltage */
	    double v[nseg];       // Membrane voltage 

	/* Time and spatial Step */
	    double dt;      // Time step 
        double t;       // Time 

// define ionic conductances (S/cm^2)
        double gna11 = 2.0;
        double gcatbar =0.45;    //0.45
        double gcat31 = 0.45;    //0.45
        double gcat32 = 0.45;    //0.45
        double gcat33 = 0.45;    //0.45
        double gcatbar_C4565 = 1.6*0.45;
        double gcanbar = 0.0025;
        double gcalbar = 0.0025;
        double gkdr =0.08;  // 0.08
        double gka =0.0001;
        double gkm =0.00002;
        double gkc =0.00005;  
        double gkahp =0.0000018;  // 0.0000018 testing
        double gpas = 1/Rm;

// define ionic concentration  & currents
        double cai[nseg];
        double ina[nseg], ina11[nseg];
	    double icat[nseg],icat_alpha1G[nseg],icat_alpha1H[nseg],icat_alpha1I[nseg];
        double iC4565[nseg];
        double ican[nseg],ical[nseg],ica[nseg];
        double ikdr[nseg],ika[nseg],ikm[nseg],ikc[nseg],ikahp[nseg],ipas[nseg];
        double ion[nseg],current[nseg];
        double Ist;

// gating variables
// Na1.1 gates
       double Na_MC3[nseg], Na_MC[nseg], Na_MO[nseg], Na_MI[nseg];
       double Na_MIs[nseg], Na_MC2[nseg], IC3[nseg], IC2[nseg]; 
       double Na_MIs2[nseg], M_C3[nseg], M_C2[nseg], M_C[nseg], M_O[nseg]; 
       double M_OP[nseg], Na_OP[nseg], Na_O[nseg], OPs[nseg], M_IF[nseg];

double ina11_wt = 0.0;
double ina11_mt = 0.0;

// Ca N and L type channel gates
       double mn[nseg],hn[nseg],ml[nseg];

// CaT gates
       double nth[nseg],lth[nseg],ntg[nseg],ltg[nseg],nti[nseg],lti[nseg];

// CaT mutant gates
       double ntC4565[nseg],ltC4565[nseg];
     
// K gates
       double nkdr[nseg],lkdr[nseg],nka[nseg],lka[nseg],mkm[nseg],okc[nseg],wkahp[nseg];       

// gloable variables
       int seg; // segment number from 1: number of segments
// geometry parameters
       double area[nseg], diam[nseg], dx[nseg], g[nseg];

void WT_SCN1A_function ();
void cat_human_alpha1G();
void cat_human_alpha1H();
void cat_human_alpha1I();
void C4565();
void can();
void cal();
void kdr();
void ka();
void km();
void kc();
void kahp_new();
void pas();
void ca_dyn();
void geo();
void cable();

int main()
{
	/* Opening of Datafiles */
	out_one = fopen("out_one", "w");
	printdata = 50;
	printval = 10;

	t = 0.0;        
    double tstop = 200;
    dt = 0.01;

    for(int i = 0; i<= nseg-1; i++)
    {
      v[i] = vrest;
      cai[i] = cai_init;
     }

// define neuron geometry: CA3 neuron
	
	geo();

	/* Beginning of Time Loop */
	for ( t = 0.0; t <= tstop; t+=dt)
	{
          if(t>=2 && t <= 7)
          {
          Ist = -0.2*100/(PI*diam[0]*dx[0]); //(mA/cm^2)
//         Ist = 0.0; // testing
          }
          else
          Ist = 0.0;

         // calculate ionic currents
          for(seg = 0; seg <= nseg-1; seg++) 
          {
		        WT_SCN1A_function ();
                cat_human_alpha1G();
                cat_human_alpha1H();
                cat_human_alpha1I();
                C4565();
                cal();
                can();
                kdr();
                ka();
                km();
                kc();
                kahp_new();
                pas();
                ca_dyn();
// debugging
/*   cout << "segment ######################################  " << seg << endl;
   cout << " The time is at ------------------------------  "<< t << endl;
   cout << v[seg] << endl;
   cout << ina[seg]<<endl;
   cout << icat[seg]<<endl;
   cout << ican[seg]<<endl;
   cout << ical[seg]<<endl;
   cout << ikdr[seg]<<endl;
   cout << ika[seg]<<endl;
   cout << ikm[seg]<<endl;
   cout << ikc[seg]<<endl;
   cout << ikahp[seg]<<endl;
   cout << ipas[seg]<<endl;
   cout << ica_pump[seg]<<endl;
   cout << cai[seg]<< endl;
*/
  
          }
    
        cable();

	// output results to a file at certain intervals
	    if(printdata >= printval) 
		{data_file();	
		printdata = 0;}
			
		printdata = printdata+1;

		cout << "t, v[0] = " << t << "     " << v[0] <<endl;

	}


	// indicate the end of the calculation: beeping
	cout << "\a\a"; 

	return (1);
}

// [Ca2+] dynamics using a simplified form

void ca_dyn()
{

  cai[0] = cai[0] +0.1*dt*(-0.15*ica[0]+(-cai[0])/5);

}

void data_file()
{
  fprintf(out_one,"%.3f\t%g\t%g\t%g\t%g\t%g\t%g\t%g\t%g\t%g\n", t, v[0],ina11[0],icat[0], ican[0], ical[0],icat_alpha1G[0],icat_alpha1H[0],icat_alpha1I[0],cai[0]);

}


// CA3 hippocampal neuron geometry (single-compartment model)
void geo() 
{
// define dimension
   double r3 = 15, L3 = 25.5; //soma

// define segment length, diameter, area assumng dx = 1, and resistances

   diam[0] = 2*r3;

   dx[0] = L3/1;

   area[0] = PI*diam[0];

// define resistances
   g[0] = 4*Ra*dx[0]/(PI*diam[0]*diam[0]);

}

//solve cable equation
void cable() 
{

  icat[0] = icat_alpha1G[0]+icat_alpha1H[0]+icat_alpha1I[0];
  ica[0] = icat[0] + ican[0] + ical[0];
  ion[0] = ina11[0]+icat[0]+ican[0]+ical[0]+ikdr[0]+ika[0]+ikm[0]+ikc[0]+ikahp[0]+ipas[0];
  current[0] = ion[0];


// inject stimulation current at soma
  current[0] = ion[0]+Ist; //(mA/cm^2)

  v[0] = v[0]-dt*current[0]/cm;

}


//**************************************************************************************
//	FUNCTION TO COMPUTE SCN1A NA CURRENT--Markov Model By Colleen Clancy @ Cornell U
//**************************************************************************************

void WT_SCN1A_function ()	
{	
	
//	Ena=(r*T/F)*log(Naout/Nain);

double a11, a12, a13, b11, b12, b13, a2, a3, a4, b4, b1, b2, b3, a5, b5, mu1, mu2, aa2, bb2, bb3, aa3;
double err1, err2, err3,err4, err5, err6, err7, err8, err9, err10, err11, err12, err13, err14; 

double  y1, z1,  x1, x2, y2, z2, e1, e2,  q1, q2,  w1, w2,  d1, d2, c1, c2, f1, j2, f2, tap1, tap2, jj1, u2, u1, p1, p2, t1, zz1, zz2;

double dny, dnz, dne1, dne, dne2, dny1, dnz1, dnx, dnx1,  dnq1, dnq,dnw, dnw1, dnd, 
 dnd1, dnc, dnc1, dnu, dnu1,  dnp, dnp1,  dntap, dntap1,  dnf1, dnj1, dnf, dnj, dnzz, dnzz1;

double qt = pow(2.0,(celsius-23.0)/10);

int i4;

 //Nain=10.0; 
 //Naout= 140.0;
//double ina;
			
int which_channels = 1;

if (which_channels ==1)   // WT rates
{
    a11= qt*2.802/(0.21*exp(-v[seg]/17.0)+0.23*exp(-v[seg]/150));
    a12= qt*2.802/(0.23*exp(-v[seg]/15.0)+0.25*exp(-v[seg]/150));
    a13= qt*2.802/(0.25*exp(-v[seg]/12.0)+0.27*exp(-v[seg]/150));
  

     b11= qt*0.4*exp(-v[seg]/20.3);
     b12= qt*0.4*exp(-(v[seg]-5)/20.3);
     b13= qt*0.4*exp(-(v[seg]-10)/20.3)/4.5;
     a3= qt*(3.7933e-7*exp(-v[seg]/7.6))*3;
     b3= qt*(0.0084+.00002*v[seg]);
    a2= qt*((9.178*exp(v[seg]/29.68))/4.5);
     b2= ((a13*a2*a3)/(b13*b3));	
     a4= (a2/100)*1.5;
      b4 = a3/5;
     a5= (a2/95000)*80.0;
      b5 = (a3/30)/10.0;

mu1= 0.0;
mu2= 0.0;


if (t<dt)
     {
     
ina11[seg] = 0.0;	
//initialize the Na+ current

//@-65
Na_MC3[seg]=0.438311;
Na_MO[seg]=2.1544e-6;
Na_MI[seg]=0.000187751;
Na_MIs[seg]=0.000501213;
Na_MC2[seg]=0.0125341;
IC3[seg]=0.527911;
IC2[seg]=0.0150964;
Na_MIs2[seg]=0.00529984;
M_C3[seg]=0;
M_C2[seg]=0;
M_C[seg]=0;
M_O[seg]=0;
M_IF[seg]=0;


   Na_MC[seg]= 1-(Na_MO[seg]+Na_MI[seg]+Na_MC2[seg]+Na_MC3[seg]+M_C3[seg]+Na_MIs[seg]
              +M_C2[seg]+M_C[seg]+M_O[seg]+IC3[seg]+IC2[seg]+Na_MIs2[seg]);
   Na_MC[seg] = 0.000155888;
        
 	}

} 

else if (which_channels ==2)   // RH rates
{
    a11= 2.802/(0.21*exp(-v[seg]/17.0)+0.23*exp(-v[seg]/150));
    a12= 2.802/(0.23*exp(-v[seg]/15.0)+0.25*exp(-v[seg]/150));
    a13= 2.802/(0.25*exp(-v[seg]/12.0)+0.27*exp(-v[seg]/150));
  

     b11= 0.4*exp(-v[seg]/20.3);
     b12= 0.4*exp(-(v[seg]-5)/20.3);
     b13= 0.4*exp(-(v[seg]-10)/20.3)/4.5;
     a3= (3.7933e-7*exp(-v[seg]/7.6))*3;
     b3= (0.0084+.00002*v[seg]);
     a2= ((9.178*exp(v[seg]/29.68))/4.5);
     b2= ((a13*a2*a3)/(b13*b3));	
     
     a4= (a2/100)*1.5;
     b4 = a3/5;
      
     a5= (a2/95000)*80.0;
     b5 = (a3/30)/10.0;

// flicker inactivation rates

aa2= ((9.178*exp(v[seg]/29.68))/4.5);
bb3= (0.0084+.00002*v[seg]);
aa3= (3.7933e-7*exp(-v[seg]/7.7))*83;

bb2= ((a13*aa2*aa3)/(b13*bb3));

mu1= 2e-5;
mu2= 2e-4;


if (t<dt)
     {
     
ina11[seg] = 0.0;				//initialize the Na+ current

//@ -65mv

Na_MC3[seg]=0.419443;
Na_MO[seg]=2.06164e-6;
Na_MI[seg]=0.000179635;
Na_MIs[seg]=0.00050059;
Na_MC2[seg]=0.0119946;
IC3[seg]=0.505087;
IC2[seg]=0.0144437;
Na_MIs2[seg]=0.00509992;
M_C3[seg]=0.0418866;
M_C2[seg]=0.00119781;
M_C[seg]=1.4897e-5;
M_O[seg]=2.05662e-7;
M_IF[seg]=7.12628e-7;


   Na_MC[seg]= 1-(Na_MO[seg]+Na_MI[seg]+Na_MC2[seg]+Na_MC3[seg]+M_C3[seg]
               +Na_MIs[seg]+M_C2[seg]+M_C[seg]+M_O[seg]+IC3[seg]+IC2[seg]+Na_MIs2[seg]+M_IF[seg]);
   
        
 	}

} 


/*  	if (t>=60000.0-dt&&Get_steady_state==1)
 	{
 	
 	cout<<"Na_MC3="<<Na_MC3<<";"<<endl;
 	cout<<"Na_MO="<<Na_MO<<";"<<endl;
 	cout<<"Na_MI="<<Na_MI<<";"<<endl;
 	cout<<"Na_MIs="<<Na_MIs<<";"<<endl;
 	cout<<"Na_MC2="<<Na_MC2<<";"<<endl;
 	cout<<"IC3="<<IC3<<";"<<endl;
 	cout<<"IC2="<<IC2<<";"<<endl;
	cout<<"Na_MIs2="<<Na_MIs2<<";"<<endl;
 	
 	cout<<"M_C3="<<M_C3<<";"<<endl;
 	cout<<"M_C2="<<M_C2<<";"<<endl;
 	cout<<"M_C="<<M_C<<";"<<endl;
 	cout<<"M_O="<<M_O<<";"<<endl;
        cout<<"M_IF="<<M_IF<<";"<<endl;
 	
 	exit (0);
 	}
        
*/  
   
 if (t>=dt)
   {	
     y1=Na_MO[seg];     z1=Na_MI[seg];	x1= Na_MC2[seg];	q1= Na_MC3[seg];	w1= Na_MIs[seg];
     d1= M_C3[seg]; 	   c1= M_C2[seg];	u1= M_C[seg];	p1= M_O[seg];	tap1= Na_MC[seg];
     jj1= IC3[seg];	   f1= IC2[seg];	e1 = Na_MIs2[seg];	zz1= M_IF[seg];
     
     
     //UPPER STATES

     Na_MC[seg]= 1-(Na_MO[seg]+Na_MI[seg]+Na_MC2[seg]+Na_MC3[seg]+M_C3[seg]
                 +Na_MIs[seg]+M_C2[seg]+M_C[seg]+M_O[seg]+IC3[seg]+IC2[seg]+Na_MIs2[seg]+M_IF[seg]);
     
     dny=((Na_MC[seg]*a13+Na_MI[seg]*b2+M_O[seg]*mu2)- (Na_MO[seg]*(b13+a2+mu1)))*dt;
     dnz=((IC2[seg]*a12+Na_MO[seg]*a2+Na_MC[seg]*b3+Na_MIs[seg]*b4)- (Na_MI[seg]*(b2+a3+a4+b12)))*dt;
     dnx= ((IC2[seg]*a3+Na_MC[seg]*b12+Na_MC3[seg]*a11+M_C2[seg]*mu2)-(Na_MC2[seg]*(a12+b11+mu1+b3)))*dt;
     dnq= (IC3[seg]*a3+Na_MC2[seg]*b11+M_C3[seg]*mu2-(Na_MC3[seg]*(a11+mu1+b3)))*dt;
     dnw= ((Na_MI[seg]*a4+Na_MIs2[seg]*b5)-Na_MIs[seg]*(b4+a5))*dt;
//     dntap= (Na_MC2[seg]*a12+Na_MO[seg]*b13+M_C[seg]*mu2+Na_MI[seg]*a3-(Na_MC[seg]*(a13+b12+mu1+b3)))*dt;
     dnj= ((Na_MC3[seg]*b3+IC2[seg]*b11)-IC3[seg]*(a11+a3))*dt;
     dnf= ((Na_MC2[seg]*b3+IC3[seg]*a11+Na_MI[seg]*b12)-IC2[seg]*(a12+b11+a3))*dt;
     dne= (Na_MIs[seg]*a5-Na_MIs2[seg]*b5)*dt;
     
     //LOWER STATES
     dnd= (Na_MC3[seg]*mu1+M_C2[seg]*b11-(M_C3[seg]*(a11+mu2)))*dt;
     dnc=  (Na_MC2[seg]*mu1+M_C3[seg]*a11+M_C[seg]*b12-(M_C2[seg]*(a12+mu2+b11)))*dt;
     dnu= (Na_MC[seg]*mu1+M_C2[seg]*a12+M_O[seg]*b13+M_IF[seg]*aa3-(M_C[seg]*(a13+mu2+b12+bb3)))*dt;
     dnp= ((M_C[seg]*a13+Na_MO[seg]*mu1+M_IF[seg]*bb2)-(M_O[seg]*(b13+mu2+aa2)))*dt;
     dnzz= ((M_O[seg]*aa2+M_C[seg]*bb3)- (M_IF[seg]*(bb2+aa3)))*dt;
     
//   cout<<dny+dnz+dnx+dnq+dnw+dntap+dnj+dnf+dnd+dnc+dnu+dnp+dne+dnzz<<endl; 		//Use this as a test that derivatives sum to zero
     
     
    y2=Na_MO[seg]+dny;
    z2=Na_MI[seg]+dnz;
    x2= Na_MC2[seg]+dnx;
    q2= Na_MC3[seg]+dnq;
    w2= Na_MIs[seg]+dnw;
//    tap2 = Na_MC[seg]+dntap;
    j2= IC3[seg]+dnj;
    f2= IC2[seg]+dnf;
    e2 = Na_MIs2[seg]+dne;
    
    d2= M_C3[seg]+dnd;
    c2= M_C2[seg]+dnc;
    u2= M_C[seg]+dnu;
    p2= M_O[seg]+dnp;
    zz2= M_IF[seg]+dnzz;
       
    
   
    
  // Na_MC[seg]= 1-(y2+z2+x2+q2+w2+d2+c2+u2+p2);

err1= y2-y1;	err2= z2-z1;	err3=x2-x1;  	err10= 0.0;
err4= q2-q1;	err5= w2-w1;	err6=d2-d1;	err11= j2-jj1;
err7= c2-c1;	err8= u2-u1;	err9= p2-p1;	err12= f2- f1;
err13= e2-e1;	err14= zz2-zz1;

	dny1= dny;
	dnz1=dnz;
	dnx1= dnx;		//UPPER STATES
	dnq1= dnq;
	dnw1=dnw;
//	dntap1=dntap;
	dnj1= dnj;
	dnf1= dnf;
	dne1= dne;
	
	
	dnd1= dnd;
	dnc1= dnc;		//LOWER STATES
	dnu1= dnu;
	dnp1= dnp;
        dnzz1= dnzz;
	
i4=0;

//cout << "testing is ok" << endl;

while (((err1>1e-5)||(err1<-1e-5)||(err2>1e-5)||(err2<-1e-5)||(err3>1e-5)
||(err3<-1e-5)||(err4>1e-5)||(err4<-1e-5)||(err5>1e-5)||(err5<-1e-5)||
(err6>1e-5)||(err6<-1e-5)||(err7>1e-5)||(err7<-1e-5)||(err8>1e-5)||
(err8<-1e-5)||(err9>1e-5)||(err9<-1e-5)||(err10>1e-5)||(err10<-1e-5)||
(err11>1e-5)||(err11<-1e-5)||(err12>1e-5)||(err12<-1e-5)||(err13>1e-5)||(err13<-1e-5)||(err14>1e-5)||(err14<-1e-5))
&&(i4<40))
 {
	y1=y2; 	 	z1=z2;	 	 x1=x2;		q1=q2;		w1= w2;
	d1= d2;		c1= c2;		u1= u2;		p1= p2;    //	tap1= tap2;
	jj1= j2;	f1= f2;		e1= e2;		zz1= zz2;

    tap1 = 1-y1-z1-x1-q1-w1-d1-c1-u1-p1-jj1-f1-e1-zz1;
	
    dny=((tap1*a13+z1*b2+p1*mu2)- (y1*(b13+a2+mu1)))*dt;
    dnz=((f1*a12+y1*a2+tap1*b3+w1*b4)- (z1*(b2+a3+a4+b12)))*dt;
    dnx=  (((f1*a3+tap1*b12+q1*a11+c1*mu2)-(x1*(a12+b11+mu1+b3)))*dt);
    dnq= ((jj1*a3+x1*b11+d1*mu2-q1*(a11+mu1+b3))*dt);
    dnw= ((z1*a4+e1*b5)-w1*(b4+a5))*dt;
//    dntap= (x1*a12+y1*b13+u1*mu2+z1*a3-(tap1*(a13+b12+mu1+b3)))*dt;
    dnj= ((q1*b3+f1*b11)-jj1*(a11+a3))*dt;
    dnf= ((x1*b3+jj1*a11+z1*b12)-f1*(a12+b11+a3))*dt;
    dne= (w1*a5- e1*b5)*dt;
    
  // cout<<dntap<<setw(15)<<t<<endl;
    
    
    dnd= (q1*mu1+c1*b11-(d1*(a11+mu2)))*dt;
    dnc= (x1*mu1+d1*a11+u1*b12-(c1*(a12+mu2+b11)))*dt;
    dnu= (tap1*mu1+c1*a12+p1*b13+zz1*aa3-(u1*(a13+mu2+b12+bb3)))*dt;
    dnp= (u1*a13+y1*mu1+zz1*bb2-(p1*(mu2+b13+aa2)))*dt;
    dnzz= (p1*aa2+u1*bb3-(zz1*(aa3+bb2)))*dt;
    
    
    
    dny= (dny+dny1)/2;
    dnz= (dnz+dnz1)/2;
    dnx= (dnx+dnx1)/2;
    dnq= (dnq+dnq1)/2;
    dnw= (dnw+dnw1)/2;
//    dntap= (dntap+dntap1)/2;
    dnj= (dnj+dnj1)/2;
    dnf= (dnf+dnf1)/2;
    dne= (dne+dne1)/2;
    
    dnd= (dnd+dnd1)/2;
    dnc= (dnc+dnc1)/2;
    dnu= (dnu+dnu1)/2;
    dnp= (dnp+dnp1)/2;
    dnzz= (dnzz+dnzz1)/2;
    
     y2= Na_MO[seg]+dny;
     z2= Na_MI[seg]+dnz;
     x2= Na_MC2[seg]+dnx;
     q2= Na_MC3[seg]+dnq;
     w2= Na_MIs[seg]+dnw;
//     tap2= Na_MC[seg]+dntap;
     j2= IC3[seg]+dnj;
     f2= IC2[seg]+dnf;
     e2= Na_MIs2[seg]+dne;
     
     d2= M_C3[seg]+dnd;
     c2= M_C2[seg]+dnc;
     u2= M_C[seg]+dnu;
     p2= M_O[seg]+dnp;
     zz2= M_IF[seg]+dnzz;
     
     dny1=dny;
     dnz1=dnz;
     dnx1= dnx;
     dnq1=dnq;
     dnw1= dnw;
//     dntap1= dntap;
     dnf1= dnf;
     dnj1= dnj;
     dne1= dne;
     
     dnd1= dnd;
     dnc1= dnc;
     dnu1= dnu;
     dnp1= dnp;
     dnzz1= dnzz;
     
    // Na_MC[seg]= 1-(y2+z2+x2+q2+w2+d2+c2+u2+p2);
     
	err1=y2-y1;
 	err2= z2-z1;
 	err3= x2-x1;
 	err4= q2-q1;
 	err5= w2-w1;
 	err6= d2-d1;
 	err7= c2-c1;
 	err8= u2-u1;
 	err9= p2-p1;
 	err10= 0.0;
 	err11= j2-jj1;
 	err12= f2-f1;
 	err13= e2-e1;
        err14= zz2-zz1;
 	
 	i4++;
	
   
 }
 
 
 if (i4<40)	
    { 
   Na_MO[seg]=y2; 	 Na_MI[seg]=z2;	 Na_MC2[seg]= x2;	Na_MC3[seg]= q2;	Na_MIs[seg]=w2;
   M_C3[seg]= d2;   	 M_C2[seg]= c2;	 M_C[seg]= u2;	M_O[seg]= p2;    //	Na_MC[seg]=tap2;
   IC3[seg]= j2;	IC2[seg]=f2;		Na_MIs2[seg]= e2; 	M_IF[seg]= zz2;

   Na_MC[seg]= 1-(Na_MO[seg]+Na_MI[seg]+Na_MC2[seg]+Na_MC3[seg]+M_C3[seg]
               +Na_MIs[seg]+M_C2[seg]+M_C[seg]+M_O[seg]+IC3[seg]+IC2[seg]+Na_MIs2[seg]+M_IF[seg]);

// cout << "Na_Mc" << Na_MC[seg] << endl;
   
  //  Na_MC[seg]=1-(Na_MO[seg]+Na_MI[seg]+ Na_MC2[seg]+Na_MIs[seg]+Na_MC3[seg]+M_C3[seg]+M_C2[seg]+M_C[seg]+M_O[seg]);
  
    }
    
   
    
  else 
  { 
//   cout<<"error_ina"<<setw(15)<<i4<<setw(15)<<t<<endl;
   exit(0);
    }
  
    
  
   
   } 
    
   M_OP[seg]= Na_MO[seg]+M_O[seg];
   M_OP[seg]= (M_OP[seg]*1.0);
   Na_OP[seg]= Na_O[seg]*0;
    
   OPs[seg]= M_OP[seg]+Na_OP[seg];
   ina11[seg]= gna11*0.035*OPs[seg]*(v[seg]-50);

if(which_channels == 1)
  ina11_wt = ina11[seg];
else if (which_channels == 2)
  ina11_mt = ina11[seg];		
    
}



// Human T-type calcium channels based on Chemin et al., J Physiol., 2002, 540, 3-14.
// CaT--alpha 1G CaV3.1
void cat_human_alpha1G() {
   double n_inf,l_inf,tau_n,tau_l,vv;
   double gbar = 0.008; // s/cm2
   double vhalfn = -49.3;  //mv
   double vhalfl = -74.2;
   double kn = 4.6;
   double kl = -5.5;
   double q10 = 2.3;
   double qt = pow(q10,(celsius-28)/10);
   double pcabar = 0.0001; // cm/s
   double z= 2;
   if(t<dt) { 
      vv = -65;
      n_inf = 1/(1+exp(-(vv-vhalfn)/kn));
      l_inf = 1/(1+exp(-(vv-vhalfl)/kl));
              
      if(vv > -56)
      tau_n = (0.8+0.025*exp(-vv/14.5))/qt;
      else    
      tau_n = (1.71*exp((vv+120)/38))/qt;
              
      tau_l = (12.3+0.12*exp(-vv/10.8))/qt;
      ntg[seg] = n_inf;
      ltg[seg] = l_inf;
   }          
              
   n_inf = 1/(1+exp(-(v[seg]-vhalfn)/kn));
   l_inf = 1/(1+exp(-(v[seg]-vhalfl)/kl));
              
   if(v[seg]>-56)
   tau_n = (0.8+0.025*exp(-v[seg]/14.5))/qt;
   else       
   tau_n = (1.71*exp((v[seg]+120)/38))/qt;
   if(v[seg]>-60)
   tau_l = (12.3+0.12*exp(-v[seg]/10.8))/qt;
   else       
   tau_l = 137/qt; 
              
   ntg[seg] = (ntg[seg] + dt*n_inf/tau_n)/(1+dt/tau_n);
   ltg[seg] = (ltg[seg] + dt*l_inf/tau_l)/(1+dt/tau_l);
              
   double w = v[seg]*0.001*z*F/(R*(celsius+273.16));
   double ghk = -0.001*z*F*(cao-cai[seg]*exp(w))*w/(exp(w)-1);
   icat_alpha1G[seg] = gcat31*pcabar*ntg[seg]*ntg[seg]*ltg[seg]*ghk;
}            
// CaT--alpha 1H  CaV3.2
// based on Vitko et al., 2005, J Neurosci, 25, 4844-4855
// provided by Dr Edward Perez-Reyes on 6/30/05
void cat_human_alpha1H() {
   double n_inf,l_inf,tau_n,tau_l,vv;
   double gbar = 0.008; // s/cm2
   double vhalfn = -48.4;  //mv
   double vhalfl = -75.6;
   double kn = 5.2;
   double kl = -6.2;
   double q10 = 2.3;
   double qt = pow(q10,(celsius-28)/10);
   double pcabar = 0.0001; // cm/s
   double z= 2;

   if(t<dt) {
      vv = -65;
      n_inf = 1/(1+exp(-(vv-vhalfn)/kn));
      l_inf = 1/(1+exp(-(vv-vhalfl)/kl));

      if(vv > -56)
      tau_n = (1.34+0.035*exp(-vv/11.8))/qt;
      else
      tau_n = (2.44*exp((vv+120)/40))/qt;

      tau_l = (18.3+0.005*exp(-vv/6.2))/qt;

//      tau_n = (1.9+1.0/(exp((vv+21.0)/11.9)+exp(-(vv+115)/21)))/qt;
//      tau_l = (18.3+(1942+exp((vv+164)/9.2))/(1+exp((vv+89.3)/3.7)))/qt;

      nth[seg] = n_inf;
      lth[seg] = l_inf;
   }

      n_inf = 1/(1+exp(-(v[seg]-vhalfn)/kn));
      l_inf = 1/(1+exp(-(v[seg]-vhalfl)/kl));

      if(v[seg] > -56)
      tau_n = (1.34+0.035*exp(-v[seg]/11.8))/qt;
      else
      tau_n = (2.44*exp((v[seg]+120)/40))/qt;
     if(v[seg]>-60)
      tau_l = (18.3+0.005*exp(-v[seg]/6.2))/qt;
     else
      tau_l = 500/qt;
//      tau_n = (1.9+1.0/(exp((v+21.0)/11.9)+exp(-(v+115)/21)))/qt;
//      tau_l = (18.3+(1942+exp((v+164)/9.2))/(1+exp((v+89.3)/3.7)))/qt;

   nth[seg] = (nth[seg] + dt*n_inf/tau_n)/(1+dt/tau_n);
   lth[seg] = (lth[seg] + dt*l_inf/tau_l)/(1+dt/tau_l);

   double w = v[seg]*0.001*z*F/(R*(celsius+273.16));
   double ghk = -0.001*z*F*(cao-cai[seg]*exp(w))*w/(exp(w)-1);
   icat_alpha1H[seg] = gcat32*pcabar*nth[seg]*nth[seg]*lth[seg]*ghk;

}

// Ca3.2 mutant channel

void C4565() {
   double nC4565_inf,lC4565_inf,tauC4565_n,tauC4565_l, carev;
   double vhalfnC4565 = -59.36;  //mv
   double vhalflC4565 = -85.5;
   double knC4565 = 7.53;
   double klC4565 = -7.18;
   double knl = 1.143;
   double q10 = 2.3;
   double qt = pow(q10,(celsius-25)/10);
   double pcabar = 0.0001; // cm/s
   double z= 2;


   if(t<dt) {
      v[seg] = vinit;
      nC4565_inf = 1/(1+exp(-(v[seg]-vhalfnC4565)/knC4565));
      lC4565_inf = 1/(1+exp(-(v[seg]-vhalflC4565)/klC4565));
      tauC4565_n = (1.52+1.0/(exp((v[seg]+40.3)/10.06)+exp(-(v[seg]+144.3)/28.6)))/qt;
      tauC4565_l = (13.7+(1943+exp((v[seg]+164)/9.2))/(1+exp((v[seg]+89.3)/3.7)))/qt;
      ntC4565[seg] = nC4565_inf;
      ltC4565[seg] = lC4565_inf;
   }

      nC4565_inf = 1/(1+exp(-(v[seg]-vhalfnC4565)/knC4565));
      lC4565_inf = 1/(1+exp(-(v[seg]-vhalflC4565)/klC4565));
      tauC4565_n = (1.52+1.0/(exp((v[seg]+40.3)/10.06)+exp(-(v[seg]+144.3)/28.6)))/qt;
      tauC4565_l = (13.7+(1943+exp((v[seg]+164)/9.2))/(1+exp((v[seg]+89.3)/3.7)))/qt;

      ntC4565[seg] = (ntC4565[seg] + dt*nC4565_inf/tauC4565_n)/(1+dt/tauC4565_n);
      ltC4565[seg] = (ltC4565[seg] + dt*lC4565_inf/tauC4565_l)/(1+dt/tauC4565_l);

//   carev = (1e3)*(R*(celsius+273.15))/(2*F)*log(cao/cai[seg]);
//   iC4565[seg] = 0.0001*gcatbar_C4565*ntC4565[seg]*ntC4565[seg]*ltC4565[seg]*(v[seg]-carev);

   double w = v[seg]*0.001*z*F/(R*(celsius+273.16));
   double ghk = -0.001*z*F*(cao-cai[seg]*exp(w))*w/(exp(w)-1);
   iC4565[seg] = gcatbar_C4565*pcabar*ntC4565[seg]*ntC4565[seg]*ltC4565[seg]*ghk;

}


// CaT--alpha 1I CaV3.3
void cat_human_alpha1I() {
   double n_inf,l_inf,tau_n,tau_l,vv;
   double gbar = 0.008; // s/cm2
   double vhalfn = -41.5;  //mv
   double vhalfl = -69.8;
   double kn = 6.2;
   double kl = -6.1;
   double q10 = 2.3;
   double qt = pow(q10,(celsius-28)/10);
   double pcabar = 0.0001; // cm/s
   double z= 2;

   if(t<dt) {
      vv = -65;
      n_inf = 1/(1+exp(-(vv-vhalfn)/kn));
      l_inf = 1/(1+exp(-(vv-vhalfl)/kl));
//      tau_n = (9.3+8.0/(exp((vv+30.0)/11.9)+exp(-(vv+120)/21)))/qt;
//      tau_l = (13.7+(1942+exp((vv+164)/9.2))/(1+exp((vv+89.3)/3.7)))/qt;

      if(vv > -60)
      tau_n = (7.2+0.02*exp(-vv/14.7))/qt;
      else
      tau_n = (0.875*exp((vv+120)/41))/qt;
      if(vv>-70)
      tau_l = (79.5+2.0*exp(-vv/9.3))/qt;
      else
      tau_l = 260/qt;

      nti[seg] = n_inf;
      lti[seg] = l_inf;
   }

      n_inf = 1/(1+exp(-(v[seg]-vhalfn)/kn));
      l_inf = 1/(1+exp(-(v[seg]-vhalfl)/kl));

//      tau_n = (9.3+8.0/(exp((v+30.0)/11.9)+exp(-(v+120)/21)))/qt;
//      tau_l = (13.7+(1942+exp((v+164)/9.2))/(1+exp((v+89.3)/3.7)))/qt;

      if(v[seg] > -60)
        tau_n = (7.2+0.02*exp(-v[seg]/14.7))/qt;
      else
        tau_n = (0.875*exp((v[seg]+120)/41))/qt;
     if(v[seg]>-60)
      tau_l = (79.5+2.0*exp(-v[seg]/9.3))/qt;
     else
      tau_l = 260/qt;
      nti[seg] = (nti[seg] + dt*n_inf/tau_n)/(1+dt/tau_n);
      lti[seg] = (lti[seg] + dt*l_inf/tau_l)/(1+dt/tau_l);

      double w = v[seg]*0.001*z*F/(R*(celsius+273.16));
      double ghk = -0.001*z*F*(cao-cai[seg]*exp(w))*w/(exp(w)-1);
      icat_alpha1I[seg] = gcat33*pcabar*nti[seg]*nti[seg]*lti[seg]*ghk;

}


// N-channel

void can() {  //N-channel from Jaffe et al. 1994, J of Neurophisiology
              double am,bm,ah,bh,tau_m,tau_h,m_inf,h_inf;


              if(t<dt) {
                        v[seg] = -65;
                        am = 0.1967*(-v[seg]+19.88)/(exp((-v[seg]+19.88)/10)-1);
                        bm = 0.046*exp(-v[seg]/20.73);
                        ah = 1.6e-4*exp(-v[seg]/48.4);
                        bh = 1/(exp((-v[seg]+39)/10)+1);
                        tau_m = 1/(am+bm);
                        tau_h = 1/(ah+bh);
                        m_inf = am*tau_m;
                        h_inf = ah*tau_h;
                        mn[seg] = m_inf;
                        hn[seg] = h_inf;
                        ican[seg] = 0.0;
              }


              am = 0.1967*(-v[seg]+19.88)/(exp((-v[seg]+19.88)/10)-1);
              bm = 0.046*exp(-v[seg]/20.73);
              ah = 1.6e-4*exp(-v[seg]/48.4);
              bh = 1/(exp((-v[seg]+39)/10)+1);
              tau_m = 1/(am+bm);
              tau_h = 1/(ah+bh);
              m_inf = am*tau_m;
              h_inf = ah*tau_h;

              mn[seg] = (mn[seg] + dt*m_inf/tau_m)/(1+dt/tau_m);
              hn[seg] = (hn[seg] + dt*h_inf/tau_h)/(1+dt/tau_h);

             double ki = 0.001;
             double KTF = 25/293.15*(celsius+273.15);
             double f = KTF/2;
             double nu = v[seg]/f;
             double efun;
             if(fabs(nu)<1e-4)
             efun = 1-nu/2;
             else
             efun = nu/(exp(nu)-1);

             double ghk = -f*(1-cai[seg]/cao*exp(nu))*efun;
//             double gcanbar = 0.0003;
   
             double gcan = gcanbar*mn[seg]*mn[seg]*hn[seg]*ki/(ki+cai[seg]);

             ican[seg] = gcan*ghk;

//             return ican;

}

// L-channel

void cal() {  //L-channel Jaffe et al. 1994. J Neurophisiology
              double am,bm,tau_m,m_inf; 
              if(t<dt) {
              v[seg] = - 65;
              am = 15.69*(-v[seg]+81.5)/(exp((-v[seg]+81.5)/10)-1);
              bm = 0.29*exp(-v[seg]/10.86);
              tau_m = 1/(am+bm);
              m_inf = am*tau_m;
              ml[seg] = m_inf;
              ical[seg] = 0.0;
              }

              am = 15.69*(-v[seg]+81.5)/(exp((-v[seg]+81.5)/10)-1);
              bm = 0.29*exp(-v[seg]/10.86);
              tau_m = 1/(am+bm);
              m_inf = am*tau_m;

              ml[seg]= (ml[seg] + dt*m_inf/tau_m)/(1+dt/tau_m);

              double ki = 0.001;
//              double gcalbar = 0.003;
              double gcal = gcalbar*ml[seg]*ml[seg]*(ki/(ki+cai[seg]));
              double KTF = 25/293.15*(celsius+273.15);
              double f = KTF/2;
              double nu = v[seg]/f;
              double efun;
              if(fabs(nu) < 1e-4)
              efun = 1-nu/2;
              else
              efun = nu/(exp(nu)-1); 
     
              double ghk = -f*(1-cai[seg]/cao*exp(nu))*efun;

              ical[seg] = gcal*ghk;

//              return ical;

}

// K channels

// Kdr channel (Ca-independent k channel)

void kdr() { // From Migliore et al. 1995. J Neurophysiology

              double an,bn,al,bl,tau_n,tau_l,n_inf,l_inf;

              if(t<dt) {
                        v[seg] = -65;
                        an = 0.03*exp(1.e-3*5*0.4*(v[seg]+32)*F/(R*T));
                        bn = 0.03*exp(-1.e-3*5*0.6*(v[seg]+32)*F/(R*T));
                        al = 0.001*exp(-1.e-3*2*(v[seg]+61)*F/(R*T));
                        bl = 0.001;
                        tau_n = 1/(an+bn);
                        tau_l = 1/(al+bl);
                        n_inf = an*tau_n;
                        l_inf = al*tau_l;
                        nkdr[seg] = n_inf;
                        lkdr[seg] = l_inf;
                        ikdr[seg] = 0.0;
              }


              an = 0.03*exp(1.e-3*5*0.4*(v[seg]+32)*F/(R*T));
              bn = 0.03*exp(-1.e-3*5*0.6*(v[seg]+32)*F/(R*T));
              al = 0.001*exp(-1.e-3*2*(v[seg]+61)*F/(R*T));
              bl = 0.001;
              tau_n = 1/(an+bn);
              tau_l = 1/(al+bl);
              n_inf = an*tau_n;
              l_inf = al*tau_l;

              nkdr[seg] = (nkdr[seg] + dt*n_inf/tau_n)/(1+dt/tau_n);
              lkdr[seg] = (lkdr[seg] + dt*l_inf/tau_l)/(1+dt/tau_l);

              ikdr[seg] = gkdr*nkdr[seg]*nkdr[seg]*nkdr[seg]*lkdr[seg]*(v[seg]+91);

//             return ikdr;


}


// Ka channel (Ca-independent k channel)

void ka() { // From Migliore et al. 1995. J Neurophysiology

              double an,bn,al,bl,tau_n,tau_l,n_inf,l_inf;

              if(t<dt) {
                        v[seg] = -65;
                        an = 0.02*exp(1.e-3*3*0.6*(v[seg]+33.6)*F/(R*T));
                        bn = 0.02*exp(-1.e-3*1.2*(v[seg]+33.6)*F/(R*T));
                        al = 0.08*exp(-1.e-3*4*(v[seg]+83)*F/(R*T));
                        bl = 0.08;
                        tau_n = 1/(an+bn);
                        tau_l = 1/(al+bl);
                        n_inf = an*tau_n;
                        l_inf = al*tau_l;
                        nka[seg] = n_inf;
                        lka[seg] = l_inf;
                        ika[seg] = 0.0;
              }


              an = 0.02*exp(1.e-3*3*0.6*(v[seg]+33.6)*F/(R*T));
              bn = 0.02*exp(-1.e-3*1.2*(v[seg]+33.6)*F/(R*T));
              al = 0.08*exp(-1.e-3*4*(v[seg]+83)*F/(R*T));
              bl = 0.08;
              tau_n = 1/(an+bn);
              tau_l = 1/(al+bl);
              n_inf = an*tau_n;
              l_inf = al*tau_l;

              nka[seg] = (nka[seg] + dt*n_inf/tau_n)/(1+dt/tau_n);
              lka[seg] = (lka[seg] + dt*l_inf/tau_l)/(1+dt/tau_l);

              ika[seg] = gka*nka[seg]*lka[seg]*(v[seg]+91.0);

//             return ika;

}


// Km channel (Ca-independent k channel)

void km() { // From Migliore et al. 1995. J Neurophysiology

              double am,bm,tau_m,m_inf,q10;
              q10 = pow(5.0,((celsius-23)/10));

              if(t<dt) {
                        v[seg] = -65;
                        am = q10*0.006*exp(1.e-3*10*0.06*(v[seg]+55)*F/(R*T));
                        bm = q10*0.006*exp(-1.e-3*9.4*(v[seg]+55)*F/(R*T));
                        tau_m = 1/(am+bm);
                        m_inf = am*tau_m;
                        mkm[seg] = m_inf;
                        ikm[seg] = 0.0;
              }


              am = q10*0.006*exp(1.e-3*10*0.06*(v[seg]+55)*F/(R*T));
              bm = q10*0.006*exp(-1.e-3*9.4*(v[seg]+55)*F/(R*T));
              tau_m = 1/(am+bm);
              m_inf = am*tau_m;

              mkm[seg] = (mkm[seg] + dt*m_inf/tau_m)/(1+dt/tau_m);

              ikm[seg] = gkm*mkm[seg]*(v[seg]+91.0);

//             return ikm;


}


// Kc channel (Ca-dependent k channel)

void kc() { // From Migliore et al. 1995. J Neurophysiology
              // Based on data by Moczydlowski & Latorre, J Gen Physiol. 1983, 82:511-542 
              double ao,bo,tau_o,o_inf;

              double k1 = 0.48e-3;
              double k2 = 0.13e-6;
              double F1 = 96.4853;

              if(t<dt) {
                        v[seg] = -65;
                        ao = cai[seg]*0.28/(cai[seg]+k1*exp(-2*0.84*v[seg]*F1/(R*T)));
                        bo = 0.48/(1+cai[seg]/(k2*exp(-2*1.0*v[seg]*F1/(R*T))));
                        tau_o = 1/(ao+bo);
                        o_inf = ao*tau_o;
                        okc[seg] = o_inf;
                        ikc[seg] = 0.0;
              }


              ao = cai[seg]*0.28/(cai[seg]+k1*exp(-2*0.84*v[seg]*F1/(R*T)));
              bo = 0.48/(1+cai[seg]/(k2*exp(-2*1.0*v[seg]*F1/(R*T))));
              tau_o = 1/(ao+bo);
              o_inf = ao*tau_o;

              okc[seg] = (okc[seg] + dt*o_inf/tau_o)/(1+dt/tau_o);

//             double gkc = 0.01;
               ikc[seg] = gkc*okc[seg]*(v[seg]+91);

//             return ikc;

}


// Kahp channel (Ca-dependent k channel)

void kahp_new() { // From Lyle J. Borg-Graham, 1999  formulation. parameters are from Gold et al., 2006.
        
              double aw,bw,tau_w,w_inf,tau0;
 
              if(t<dt) {
                        cai[seg]  = 5e-5;  
                        aw = 1e31;
                        bw = 0.2;
                        tau0 = 5.0;
                        tau_w = 1/(aw*(pow(cai[seg],2))+bw)+tau0;
                        w_inf = aw*pow(cai[seg],2)/(aw*(pow(cai[seg],2))+bw);
                        wkahp[seg] = w_inf;
                        ikahp[seg] = 0.0;
              }
         
         
              aw = 1e31;
              bw = 0.2;
              tau_w = 1/(aw*pow(cai[seg],2)+bw)+tau0;
              w_inf = aw*pow(cai[seg],2)/(aw*pow(cai[seg],2)+bw);
     
              wkahp[seg] = (wkahp[seg] + dt*w_inf/tau_w)/(1+dt/tau_w);
 
              ikahp[seg] = gkahp*pow(wkahp[seg],8)*(v[seg]+91);
}


// leak current or passive current

void pas() {

//       double gpas = 1/Rm;
       ipas[seg] = gpas*(v[seg]-vrest);
//       return ipas;

}



Loading data, please wait...