hhModel-v2.cpp

#include <iostream>
#include <cmath>
#include <chrono>

uint32_t num_calls = 0; // to fool optimizer

namespace constants
{
    constexpr const auto time_step = 0.01;
    constexpr const auto num_samples = 100000;

    constexpr const auto initial_voltage = -65.0;
    constexpr const auto reference_voltage = -65.0;

    constexpr const auto membrane_surface = 4000.0;  // [um^2] surface area of the membrane
    constexpr const auto membrane_capacitance_density = 1.0;   // [uF/cm^2] membrane capacitance density
    constexpr const auto membrane_capacitance = membrane_capacitance_density * membrane_surface * 1e-8;   // [uF] membrane capacitance

    constexpr const auto GNa = 120.0;
    constexpr const auto GK = 36.0;
    constexpr const auto GL = 0.3;
    constexpr const auto ENa = 125.0;
    constexpr const auto EK = -55.0;
    constexpr const auto EL = -25.0;

    constexpr const auto sodium_conductance = GNa * membrane_surface * 1e-8;// Na conductance [mS]
    constexpr const auto kalium_conductance = GK * membrane_surface * 1e-8;  // K conductance [mS]
    constexpr const auto leak_conductance = GL * membrane_surface * 1e-8;  // leak conductance [mS]
}

class simulation
{
    double membrane_potential;
    double m;
    double h;
    double n;

    double step_m(const double delta_time) const
    {
        const auto voltage_difference = membrane_potential - constants::reference_voltage;

        const auto aM = 0.1 * (voltage_difference - 25.0) / (1.0 - exp(-(voltage_difference - 25.0) / 10.0));
        const auto bM = 4.0 * exp(-voltage_difference / 18.0);

        auto const derivative_m = (1.0 - m) * aM - m * bM;

        return (m + derivative_m * delta_time);
    }

    double step_h(const double delta_time) const
    {
        const auto voltage_difference = membrane_potential - constants::reference_voltage;

        const auto aH = 0.07 * exp(-voltage_difference / 20.0);
        const auto bH = 1.0 / (1.0 + exp(-(voltage_difference - 30.0) / 10.0));

        auto const derivative_h = (1.0 - h) * aH - h * bH;

        return h + derivative_h * delta_time;
    }

    double step_n(const double delta_time) const
    {
        const auto voltage_difference = membrane_potential - constants::reference_voltage;

        const auto aN = 0.01 * (voltage_difference - 10.0) / (1.0 - exp(-(voltage_difference - 10.0) / 10.0));
        const auto bN = 0.125 * exp(-voltage_difference / 80.0);

        auto const derivative_n = (1.0 - n) * aN - n * bN;

        return n + derivative_n * delta_time;
    }

    double step_membrane_potential(const double time, const double stimulus) const
    {
        const auto sodium_current = constants::sodium_conductance * m*m*m * h * (constants::ENa - membrane_potential);
        const auto kalium_current = constants::kalium_conductance * n*n*n*n * (constants::EK - membrane_potential);
        const auto leak_current = constants::leak_conductance * (constants::EL - membrane_potential);

        auto const derivative_potential = (sodium_current + kalium_current + leak_current + stimulus) / constants::membrane_capacitance;

        return membrane_potential + derivative_potential * time;
    }
public:
    simulation(const double initial_potential, const double initial_m, const double initial_h, const double initial_n) : membrane_potential(initial_potential), m(initial_m), h(initial_h), n(initial_n) {}

    void step(const double delta_time, const double current_stimulus)
    {
        membrane_potential = step_membrane_potential(delta_time, current_stimulus);
        m = step_m(delta_time);
        h = step_h(delta_time);
        n = step_n(delta_time);

        ++num_calls; // to fool optimizer
    }

    double latest_potential() const
    {
        return membrane_potential;
    }
};

void run_simulation()
{
    auto sim = simulation{ constants::initial_voltage, 0, 0, 0 };

    for (auto j = 0; j < constants::num_samples; j++) {
        auto stimulus = 0.0;
        if (j >= 25000 && j <= 75000) stimulus = 0.01;

        sim.step(constants::time_step, stimulus);
    }
}

int main(int argc, char *argv[])
{
    using seconds = std::chrono::duration<double, std::ratio<1, 1>>;

    if(argc < 2)
    {
        std::cerr << "requires number of benchmarking loops as command line argument #1";
        return 1;
    }

    const auto num_loops = atoi(argv[1]);

    const auto start_time = std::chrono::high_resolution_clock::now();

    for (auto loop = 0; loop < num_loops; ++loop) {
        run_simulation();
    }

    const auto end_time = std::chrono::high_resolution_clock::now();
    const auto time_taken = std::chrono::duration_cast<seconds>(end_time - start_time);

    std::cout << "time taken: " << time_taken.count() << "s\n";
    std::cout << "number of loops: " << num_loops << "\n";
    std::cout << "number of simulation steps: " << (num_calls / num_loops) << "\n"; // to fool optimizer
}