Gen_mark_mod.html

/* C <--> O <--> I

Function to compute current for simple three-state model
Colleen Clancy
Compile with cc (C++)
The function should be called from Main() which contains the
time loop and value for dt.

C --> O = a // rate constants
O --> C = b // can be adjusted for voltage or ion dependence
O --> I = aa
I --> O = bb

Global Variables (or variables passed from Main())
double C, O, I;

double Comp_current() // current function

{
double Eion=(r*T/F)*log(IONout/IONin); //channel reversal potential

double Gion = 10.0; //maximum membrane conductance

double a, aa, b, bb, err1, err2, y1, z1, y2, z2, dny, dnz, dny1, dnz1;
int i4;

   a= 10;
aa= 100;   // rate constants
    b = 10;   // can be adjusted for voltage or ion dependence
bb = 100;

   if (t==0)    // initial values
     {
    C= 1.0;
I = 0.0;
    O = 1.0- C I;
}

   else

   {
     y1=C;     z1= I;

dny= (O*b- C*a)*dt;
dnz= (O*aa- I*bb)*dt;

     y2= C+dny;
     z2= I+dnz;

    O= 1.0-y2-z2;

err1= y2-y1; err2= z2-z1;

dny1= dny;
dnz1=dnz;

i4=0;
while (((err1>1e-5)||(err1<-1e-5)||(err2>1e-5)||(err2<-1e-5))&&(i4<40))
{
y1=y2;   z1=z2;

    dny= (O*b- y1*a)*dt;
    dnz= (O*aa- z1*bb)*dt;

    dny= (dny+dny1)/2;
    dnz= (dnz+dnz1)/2;

     y2= C+dny;
     z2= I+dnz;

     dny1=dny;
     dnz1=dnz;

     O= 1.0-y2-z2;

err1=y2-y1;
err2= z2-z1;

i4++;

}

if (i4<40)
    {
   C=y2; I=z2;
    }

else
{
   cout<<"error_ina"<<setw(15)<<i4<<setw(15)<<t<<endl;
   exit(0);
    }
    O= 1-C -I;

    }

Gen_current = Gion*(O)*(v-Eion);

return Gen_current;
}