hhModel.cpp

#include <iostream>
#include <cstdlib>
#include <stdio.h>
#include <math.h>
#include <array>
using namespace std;

int main()
{
    double dt = 0.01;
    int numSamples = 100000;

    // create an rectangular external stimulus
    double I[numSamples];
    for(int i = 0; i <= numSamples; i++)
        I[i] = 0;
    for(int i = 25000; i <= 75000; i++)
            I[i] = 0.01;

    // init constants
    double Vinit = -65;
    double Vref = -65;

    double Smemb = 4000;  // [um^2] surface area of the membrane
    double Cmemb = 1.0;   // [uF/cm^2] membrane capacitance density
    double Cm = Cmemb * Smemb * 1e-8;   // [uF] membrane capacitance

    double GNa = 120;
    double GK = 36;
    double GL = 0.3;
    double ENa = 125;
    double EK = -55;
    double EL = -25;

    double v[numSamples]; // [mV] membrane potential vector

    double m, h, n, aM, bM, aH, bH, aN, bN, dv_dt, dn_dt, dm_dt, dh_dt, IL, IK, INa;

    double gNa = GNa * Smemb * 1e-8;// Na conductance [mS]
    double gK = GK * Smemb * 1e-8;  // K conductance [mS]
    double gL = GL * Smemb * 1e-8;  // leak conductance [mS]

    // initial values
    v[0] = Vinit;           // initial membrane potential
    m = 0;
    h = 0;
    n = 0;

    // To compare the code to Matlab I run the numerical integration 1000 times
    for (int nn = 0; nn <= 1000; nn++){

    // Numerical integration
    for (int j = 0; j <= numSamples; j++){
        // ionic currents
        INa = gNa * m*m*m * h * (ENa - v[j]);
        IK = gK * n*n*n*n * (EK - v[j]);
        IL = gL * (EL - v[j]);

        aM = 0.1 * (v[j] - Vref -25) / ( 1 - exp(-(v[j]-Vref-25)/10));
        bM = 4 * exp(-(v[j]-Vref)/18);

        aH = 0.07 * exp(-(v[j]-Vref)/20);
        bH = 1 / ( 1 + exp( -(v[j]-Vref-30)/10 ) );

        aN = 0.01 * (v[j]-Vref-10) / ( 1 - exp(-(v[j]-Vref-10)/10) );
        bN = 0.125 * exp(-(v[j]-Vref)/80);

        // derivatives
        dv_dt = ( INa + IK + IL + I[j] ) / Cm;
        dm_dt = (1.0-m) * aM - m * bM;
        dh_dt = (1.0-h) * aH - h * bH;
        dn_dt = (1.0-n) * aN - n * bN;

        // calculate next step
        v[j+1] = v[j] + dv_dt * dt;
        m = m + dm_dt * dt;
        h = h + dh_dt * dt;
        n = n + dn_dt * dt;
    }
    }
}