Bughunt! Now hiring! (Hodgkin Huxley Neuron)

[Horst]

Hi!

I'm quite new to C++ and as it goes I have problem of which I think is quite trivial. However I'm completely stuck.

What I want to do is create a simulation of a Hodgkin & Huxley Neuron. For all non-biologitsts, this is a very old and classic model derived from empirical data.
It consists of 4 differential equations and 6 regular equations.
I've already searched for examples, and actually found some using nearly the same values but mine doesn't work. This is quite frustrating since I don't know where to look anymore. The given values are all correct, they were taken from a working model. The equations were all triple checked. So my guess is I'm making a stupid beginners mistake on which I can't quite put my finger....
Any help/suggestion would be greatly appreciated. If you need any further information just ask.
So what the model should do is, come to a resting potential around -65mV that is Vi should become -65. However I'm stuck at around -7.

#include <math.h>
#include <fstream>
#include <iostream>
using namespace std;


//define functions --------------------------
double i_Na();
double i_K();
double i_L();

double val_n(double ni);
double val_m(double mi);
double val_h(double hi);

double val_V(double Vi);


//begin variables ------------------------------
//time
double ti,t0,tn;

//iteration
int i;
int n;

double s;

//potential
double V0,Vi;

//membrane permeability
double n0,ni;
double m0,mi;
double h0,hi;

//---permeability

double alpha_n;
double beta_n;

double alpha_m;
double beta_m;

double alpha_h;
double beta_h;

//------constants
double C;

double VNa;
double VK;
double VL;

double gNa;
double gK;
double gL;


int main(){

//begin initialising variables ----------------------

t0 = 0;               //Intervall - START
tn = 400;            //Intervall - STOP
n =  200000;            //SAMPLES
i = 0;                //Sample-count

s = (tn - t0) / n;    //step range (increment)


//constants declaration

C = 1;          //Capacity in uF*cm^-2

VNa = 50;       //Sodium equilibrium potential in mV 50
VK  = -77;      //Potassium equilibrium potential in mV -77
VL  = -54.387;  //Leak potential in mV -54.387


gNa = 120;      //Sodium conductance in mS * cm^-2
gK = 36;        //Potassium conductance in mS * cm^-2
gL = 0.3;       //Leakage 'conductance' in mS * cm^-2


//initialising starting variables

V0 = 0;          //0

n0 = 0.31768;  //0.31768
h0 = 0.59612;    //0.59612
m0 = 0.052932;   //0.052932

Vi = V0;

ni = n0;
mi = m0;
hi = h0;


//end initialising variables --------------------


ofstream daten;
daten.open ("werte.txt");
	

for (i; i <= n; i++)
{

ti = t0 + s * i;


//---data-log

daten << ti;
daten << "\t";
daten << Vi;
daten << "\n";
  

//mi ------------

mi = val_m(mi);

//ni ------------

ni = val_n(ni);

//hi ------------

hi = val_h(hi);



//constrain m,n,h ------------
if(mi <= 0.000001){
  mi = 0.000001;
}
if(mi >= 1){
  mi = 1;
}
if(ni <= 0.000001){
  ni = 0.000001;
}
if(ni >= 1){
  ni = 1;
}
if(hi <= 0.000001){
  hi = 0.000001;
}
if(hi >= 1){
  hi = 1;
}


//Membrane Potential --------------


Vi = val_V(Vi);

//--------------------------
}

	daten.close();
	cout << "done!";
	return 0;

}

//begin formulae -----------------------

//individual currents ******************

//i_Na sodium current 
double i_Na()
{
  return (gNa  * pow(mi,3) * hi * (Vi - VNa));
}
  
//i_K potassium current
double i_K()
{
  return (gK * pow(ni,4) * (Vi - VK));
}
  
//i_L leakage current
double i_L()
{
  return (gL * (Vi - VL));
}


//conductance activation and inhibition ****************************
//potassium channel activation factor
double val_n(double ni)
{
  alpha_n = (0.01 * (Vi + 10)) / exp(((Vi +10) / 10) -1);
  beta_n = 0.125 * exp(Vi / 80);
  
  return (ni + s * (alpha_n * (1 - ni) - beta_n * ni));
}

//sodium channel activation factor
double val_m(double mi)
{
  alpha_m = (0.1 * (Vi + 25)) / exp(((Vi + 25) / 10) -1);
  beta_m = 4 * exp(Vi / 18);
  
  return (mi + s * (alpha_m * (1 - mi) - beta_m * mi));
}

//sodium channel inhibition factor
double val_h(double hi)
{
  alpha_h = 0.07 * exp(Vi / 20);
  beta_h = 1 / (exp((Vi + 30) / 10) +1);
  
  return (hi + s * (alpha_h * (1 - hi) - beta_h * hi));
}



//membrane potential **************************
double val_V(double Vi)
{
return (Vi + s * ((0 - (i_Na() + i_K() + i_L())) / C));
}

All the equations and values can be found here:
http://www.ebi.ac.uk/biomodels-main/BIOMD0000000020

I know it looks quite complicated, it might be easyiest to start looking at the values of m,n and h because they're the first three differential equations depending on Vi (Volt at step i), from there on Vi directly depends on m,n,h the rest is all constants.

cheers!
Horst
(and my eternal gratitude should someone be able to help me :-))

[guestgulkan]

Fist thing first:
daten << ri; //what is this variable??

[Horst]

Whoops!
This is a variable that belongs to a function stimulating the neuron.
Basically it sets a peridoic input of current that looks like a shark fin and should cause the neuron to depolarize.
A much simpler model would look something like this: ____|____
but since there's no point in stimulating until it's stable I removed it. In my code it's also commented out.
[
return (Vi + s * ((0 - (i_Na() + i_K() + i_L())) / C));
Where now a zero is, the function would return the stimulation.
0 for no stim. around 15mV for stim.
]

[Horst]

some further info:

Variables m, n, and h determine how many channels of the current ion are active.

n^4 for example represents active potassium channels. if 1 all channels are open and the membrane potential is changing.

m^3 * h represents the same for sodium channels, with the exception that m is activation and h is inhibition.

If you want to take a look at a working model, this is how it should look like (java simulation):
http://www.afodor.net/HHModel.htm

And here you can see the behavior of h,m and n respectively to V:
http://icwww.epfl.ch/~gerstner/SPNM/node14.html