function HodgkinHuxleySim(modelParams, simParams) {
    // Init
    this.model = modelParams;
    this.sim = simParams;
    this.var = {};

    // Run one stop of the model
    this.step = function() {
        let n = {
            alpha: this.computeAlphaN(this.model.V - this.model.Eref),
            beta:  this.computeBetaN(this.model.V - this.model.Eref)
        }
        let m = {
            alpha: this.computeAlphaM(this.model.V - this.model.Eref),
            beta:  this.computeBetaM(this.model.V - this.model.Eref)
        }
        let h = {
            alpha: this.computeAlphaH(this.model.V - this.model.Eref),
            beta:  this.computeBetaH(this.model.V - this.model.Eref)
        }

        n.tau = 1 / (n.alpha + n.beta);
        m.tau = 1 / (m.alpha + m.beta);
        h.tau = 1 / (h.alpha + h.beta);

        n.inf = n.alpha * n.tau;
        m.inf = m.alpha * m.tau;
        h.inf = h.alpha * h.tau;

        this.model.n        += this.sim.dt / n.tau * (n.inf - this.model.n);
        this.model.m        += this.sim.dt / m.tau * (m.inf - this.model.m);
        this.model.h        += this.sim.dt / h.tau * (h.inf - this.model.h);

        this.model.Na.G     = this.model.Na.GMax * Math.pow(this.model.m, 3) * this.model.h;
        this.model.K.G      = this.model.K.GMax * Math.pow(this.model.n, 4);

        this.model.Na.I     = this.model.Na.G * (this.model.Na.E - this.model.V);
        this.model.K.I      = this.model.K.G * (this.model.K.E - this.model.V);
        this.model.L.I      = this.model.L.G * (this.model.L.E - this.model.V);

        this.model.V        += this.sim.dt / this.model.C * (this.model.Na.I + this.model.K.I + this.model.L.I + this.sim.Iinj);
    }

    // Run the whole model
    this.run = function() {
        this.model.V  = this.model.V0;
        this.var.Vs   = new Array(Math.floor(this.sim.dur/this.sim.dt));
        this.var.Iinj = new Array(Math.floor(this.sim.dur/this.sim.dt));
        this.var.IK   = new Array(Math.floor(this.sim.dur/this.sim.dt));
        this.var.GK   = new Array(Math.floor(this.sim.dur/this.sim.dt));
        this.var.INa  = new Array(Math.floor(this.sim.dur/this.sim.dt));
        this.var.GNa  = new Array(Math.floor(this.sim.dur/this.sim.dt));

        // Sample rate of simulation data to generate plotting data
        let plotSampleRate = Math.ceil(this.sim.dur/4000/this.sim.dt);

        // For each step of the simulation
        let i = 0;
        for (let t = 0; t < this.sim.dur; t += this.sim.dt) {
            // Determine the current state of injection
            let rawInj = 0;
            if (t > this.sim.inj1Start && t < this.sim.inj1Start+this.sim.inj1Dur) {
                this.sim.inj1Amp = this.sim.inj1Amp ? this.sim.inj1Amp : 0;
                rawInj += this.sim.inj1Amp;
            }
            if (t > this.sim.inj2Start && t < this.sim.inj2Start+this.sim.inj2Dur) {
                this.sim.inj2Amp = this.sim.inj2Amp ? this.sim.inj2Amp : 0;
                rawInj += this.sim.inj2Amp;
            }
            this.sim.Iinj = rawInj;
            
            // Advance the simulation by one step
            this.step();

            if (i == plotSampleRate) {
                this.var.Vs.push([t, this.model.V]);
                this.var.Iinj.push([t, this.sim.Iinj]);
                this.var.IK.push([t, -this.model.K.I]);
                this.var.GK.push([t, this.model.K.G]);
                this.var.INa.push([t, -this.model.Na.I]);
                this.var.GNa.push([t, this.model.Na.G]);
                i = 0;
            }

            i++;
        }
        
        this.var.EKdata = new Array;
        this.var.EKdata[0] = [0, this.model.K.E];
        this.var.EKdata[1] = [this.sim.dur, this.model.K.E];
        this.var.ENadata = new Array;
        this.var.ENadata[0] = [0, this.model.Na.E];
        this.var.ENadata[1] = [this.sim.dur, this.model.Na.E];

        let output = {};
        output.model = this.model;
        output.sim = this.sim;
        output.var = this.var;
        return output;
    }

    // Simulation functions
    this.computeAlphaN = function(V) {
        return (V === 10) ? this.computeAlphaN(V+0.001) : (10-V) / (100 * (Math.exp((10-V)/10)-1));
    };
    this.computeBetaN = function(V) {
        return 0.125 * Math.exp(-V/80);
    };
    this.computeAlphaM = function(V) {
        return (V === 25) ? this.computeAlphaM(V+0.001) : (25-V) / (10 * (Math.exp((25-V)/10)-1));
    };
    this.computeBetaM = function(V) {
        return 4 * Math.exp(-V/18)
    };
    this.computeAlphaH = function(V) {
        return 0.07*Math.exp(-V/20);
    };
    this.computeBetaH = function(V) {
        return 1 / (Math.exp((30-V)/10)+1);
    };
};