Compute the HH52 model 1,000 times

hh52.cpp

#include <math.h>
#include <stddef.h>
#include <stdlib.h>

#include <iostream>

typedef enum {
    VARIABLE_OF_INTEGRATION,
    STATE,
    CONSTANT,
    COMPUTED_CONSTANT,
    ALGEBRAIC
} VariableType;

typedef struct {
    char name[8];
    char units[16];
    char component[25];
    VariableType type;
} VariableInfo;

const size_t STATE_COUNT = 4;
const size_t VARIABLE_COUNT = 18;

const VariableInfo STATE_INFO[] = {
    {"V", "millivolt", "membrane", STATE},
    {"h", "dimensionless", "sodium_channel_h_gate", STATE},
    {"m", "dimensionless", "sodium_channel_m_gate", STATE},
    {"n", "dimensionless", "potassium_channel_n_gate", STATE}
};

const VariableInfo VARIABLE_INFO[] = {
    {"i_Stim", "microA_per_cm2", "membrane", ALGEBRAIC},
    {"i_L", "microA_per_cm2", "leakage_current", ALGEBRAIC},
    {"i_K", "microA_per_cm2", "potassium_channel", ALGEBRAIC},
    {"i_Na", "microA_per_cm2", "sodium_channel", ALGEBRAIC},
    {"Cm", "microF_per_cm2", "membrane", CONSTANT},
    {"E_R", "millivolt", "membrane", CONSTANT},
    {"E_L", "millivolt", "leakage_current", COMPUTED_CONSTANT},
    {"g_L", "milliS_per_cm2", "leakage_current", CONSTANT},
    {"E_Na", "millivolt", "sodium_channel", COMPUTED_CONSTANT},
    {"g_Na", "milliS_per_cm2", "sodium_channel", CONSTANT},
    {"alpha_m", "per_millisecond", "sodium_channel_m_gate", ALGEBRAIC},
    {"beta_m", "per_millisecond", "sodium_channel_m_gate", ALGEBRAIC},
    {"alpha_h", "per_millisecond", "sodium_channel_h_gate", ALGEBRAIC},
    {"beta_h", "per_millisecond", "sodium_channel_h_gate", ALGEBRAIC},
    {"E_K", "millivolt", "potassium_channel", COMPUTED_CONSTANT},
    {"g_K", "milliS_per_cm2", "potassium_channel", CONSTANT},
    {"alpha_n", "per_millisecond", "potassium_channel_n_gate", ALGEBRAIC},
    {"beta_n", "per_millisecond", "potassium_channel_n_gate", ALGEBRAIC}
};

double * createStatesArray()
{
    double *res = (double *) malloc(STATE_COUNT*sizeof(double));

    for (size_t i = 0; i < STATE_COUNT; ++i) {
        res[i] = NAN;
    }

    return res;
}

double * createVariablesArray()
{
    double *res = (double *) malloc(VARIABLE_COUNT*sizeof(double));

    for (size_t i = 0; i < VARIABLE_COUNT; ++i) {
        res[i] = NAN;
    }

    return res;
}

void deleteArray(double *array)
{
    free(array);
}

void initialiseVariables(double *states, double *rates, double *variables)
{
    variables[4] = 1.0;
    variables[5] = 0.0;
    variables[7] = 0.3;
    variables[9] = 120.0;
    variables[15] = 36.0;
    states[0] = 0.0;
    states[1] = 0.6;
    states[2] = 0.05;
    states[3] = 0.325;
}

void computeComputedConstants(double *variables)
{
    variables[6] = variables[5]-10.613;
    variables[8] = variables[5]-115.0;
    variables[14] = variables[5]+12.0;
}

void computeRates(double voi, double *states, double *rates, double *variables)
{
    variables[0] = ((voi >= 10.0) && (voi <= 10.5))?-20.0:0.0;
    variables[1] = variables[7]*(states[0]-variables[6]);
    variables[2] = variables[15]*pow(states[3], 4.0)*(states[0]-variables[14]);
    variables[3] = variables[9]*pow(states[2], 3.0)*states[1]*(states[0]-variables[8]);
    rates[0] = -(-variables[0]+variables[3]+variables[2]+variables[1])/variables[4];
    variables[10] = 0.1*(states[0]+25.0)/(exp((states[0]+25.0)/10.0)-1.0);
    variables[11] = 4.0*exp(states[0]/18.0);
    rates[2] = variables[10]*(1.0-states[2])-variables[11]*states[2];
    variables[12] = 0.07*exp(states[0]/20.0);
    variables[13] = 1.0/(exp((states[0]+30.0)/10.0)+1.0);
    rates[1] = variables[12]*(1.0-states[1])-variables[13]*states[1];
    variables[16] = 0.01*(states[0]+10.0)/(exp((states[0]+10.0)/10.0)-1.0);
    variables[17] = 0.125*exp(states[0]/80.0);
    rates[3] = variables[16]*(1.0-states[3])-variables[17]*states[3];
}

void computeVariables(double voi, double *states, double *rates, double *variables)
{
    variables[1] = variables[7]*(states[0]-variables[6]);
    variables[3] = variables[9]*pow(states[2], 3.0)*states[1]*(states[0]-variables[8]);
    variables[10] = 0.1*(states[0]+25.0)/(exp((states[0]+25.0)/10.0)-1.0);
    variables[11] = 4.0*exp(states[0]/18.0);
    variables[12] = 0.07*exp(states[0]/20.0);
    variables[13] = 1.0/(exp((states[0]+30.0)/10.0)+1.0);
    variables[2] = variables[15]*pow(states[3], 4.0)*(states[0]-variables[14]);
    variables[16] = 0.01*(states[0]+10.0)/(exp((states[0]+10.0)/10.0)-1.0);
    variables[17] = 0.125*exp(states[0]/80.0);
}

void printHeaders()
{
    std::cout << "voi";

    for (size_t i = 0; i < STATE_COUNT; ++i) {
        std::cout << "," << STATE_INFO[i].name;
    }

    for (size_t i = 0; i < VARIABLE_COUNT; ++i) {
        std::cout << "," << VARIABLE_INFO[i].name;
    }

    std::cout << std::endl;
}

void printValues(double voi, const double *states, const double *variables)
{
    std::cout << voi;

    for (size_t i = 0; i < STATE_COUNT; ++i) {
        std::cout << "," << states[i];
    }

    for (size_t i = 0; i < VARIABLE_COUNT; ++i) {
        std::cout << "," << variables[i];
    }

    std::cout << std::endl;
}

void runModel(bool output = true)
{
    // Create our various arrays.

    double voi = 0.0;
    double *states = createStatesArray();
    double *rates = createStatesArray();
    double *variables = createVariablesArray();

    // Initialise our states and variables, and compute our computed constants, and output the initial value/guess of
    // our states and variables.

    initialiseVariables(states, rates, variables);
    computeComputedConstants(variables);
    computeRates(voi, states, rates, variables);
    computeVariables(voi, states, rates, variables);

    if (output) {
        printHeaders();
        printValues(voi, states, variables);
    }

    // Run our model.

    double step = 0.001;
    size_t iMax = 50.0 / step;
    size_t iModulo = 0.1 / step;

    for (size_t i = 1; i <= 50000; ++i) {
        voi = i * step;

        // Integrate our model.

        computeRates(voi, states, rates, variables);

        for (size_t j = 0; j < STATE_COUNT; ++j) {
            states[j] += step * rates[j];
        }

        // Compute our variables.

        computeVariables(voi, states, rates, variables);

        // Output the value of our states and variables.

        if (output && (i % iModulo == 0)) {
            printValues(voi, states, variables);
        }
    }

    // Clean up after ourselves.

    deleteArray(states);
    deleteArray(rates);
    deleteArray(variables);
}

int main()
{
    // runModel();
    for (size_t i = 0; i < 1000; ++i) {
        runModel(false);
    }

    return 0;
}