HH-Modelica

This project contains a modular implementation of the Hodgkin-Huxley equations that has two aims:

1. Increase the understandability of the Hodgkin-Huxley model for novices by keeping the amount of variables and equations per component very low and bridging the gap between biological meaning and electrical analogy.
2. Serve as a unifying basis for extensions by experts and therefore facilitate the creation of more complex Hogdkin-Huxley-type models.

Full modular model

model HHmodular "'flat' version of the modular model (no membrane container)"
  extends PotentialAdapter(v = l2.v);
  HHmodelica.Components.PotassiumChannel c_pot annotation(
    Placement(visible = true, transformation(origin = {-33, 3}, extent = {{-17, -17}, {17, 17}}, rotation = 0)));
  HHmodelica.Components.SodiumChannel c_sod annotation(
    Placement(visible = true, transformation(origin = {1, 3}, extent = {{-17, -17}, {17, 17}}, rotation = 0)));
  HHmodelica.Components.LeakChannel c_leak annotation(
    Placement(visible = true, transformation(origin = {35, 3}, extent = {{-17, -17}, {17, 17}}, rotation = 0)));
  HHmodelica.Components.LipidBilayer l2 annotation(
    Placement(visible = true, transformation(origin = {-67, 3}, extent = {{-17, -17}, {17, 17}}, rotation = 0)));
  HHmodelica.Components.CurrentClamp clamp annotation(
    Placement(visible = true, transformation(origin = {69, 3}, extent = {{-17, -17}, {17, 17}}, rotation = 0)));
equation
  connect(l2.p, c_pot.p) annotation(
    Line(points = {{-66, 20}, {-66, 40}, {-33, 40}, {-33, 20}}, color = {0, 0, 255}));
  connect(c_pot.p, c_sod.p) annotation(
    Line(points = {{-33, 20}, {-32, 20}, {-32, 40}, {0, 40}, {0, 20}, {2, 20}}, color = {0, 0, 255}));
  connect(c_sod.p, c_leak.p) annotation(
    Line(points = {{2, 20}, {2, 20}, {2, 40}, {34, 40}, {34, 20}, {36, 20}}, color = {0, 0, 255}));
  connect(c_leak.p, clamp.p) annotation(
    Line(points = {{36, 20}, {36, 20}, {36, 40}, {68, 40}, {68, 20}, {70, 20}}, color = {0, 0, 255}));
  connect(clamp.n, c_leak.n) annotation(
    Line(points = {{70, -14}, {68, -14}, {68, -40}, {36, -40}, {36, -14}, {36, -14}}, color = {0, 0, 255}));
  connect(c_leak.n, c_sod.n) annotation(
    Line(points = {{36, -14}, {34, -14}, {34, -40}, {2, -40}, {2, -14}, {2, -14}}, color = {0, 0, 255}));
  connect(c_sod.n, c_pot.n) annotation(
    Line(points = {{2, -14}, {0, -14}, {0, -40}, {-32, -40}, {-32, -14}, {-33, -14}}, color = {0, 0, 255}));
  connect(c_pot.n, l2.n) annotation(
    Line(points = {{-33, -14}, {-33, -40}, {-66, -40}, {-66, -14}}, color = {0, 0, 255}));
  connect(l2.temp, c_pot.temp) annotation(
    Line(points = {{-58, 12}, {-58, 16}, {-42, 16}, {-42, 13}}, color = {255, 0, 0}));
  connect(c_pot.temp, c_sod.temp) annotation(
    Line(points = {{-42, 13}, {-40, 13}, {-40, 16}, {-8, 16}, {-8, 12}}, color = {255, 0, 0}));
  annotation(
    experiment(StartTime = 0, StopTime = 0.03, Tolerance = 1e-6, Interval = 1e-05),
    __OpenModelica_simulationFlags(s = "dassl"),
    __MoST_experiment(variableFilter = "v_m|clamp\\.(v|i)|c_pot\\.(g|gate_act\\.n)|c_sod\\.(g|gate_act\\.n|gate_inact\\.n)"));
end HHmodular;

1. $\mathrm{clamp.p.i}+\mathrm{l2.i}+\mathrm{c\_leak.i}+\mathrm{c\_sod.i}+\mathrm{c\_pot.i}\equiv 0.0$
2. $(((-\mathrm{l2.i}-\mathrm{clamp.cur.i\_const})-\mathrm{c\_leak.i})-\mathrm{c\_pot.i})-\mathrm{c\_sod.i}\equiv 0.0$
3. $\mathrm{clamp.g.p.i}+(\mathrm{clamp.cur.i\_const}-\mathrm{clamp.p.i})\equiv 0.0$
4. $\mathrm{c\_leak.i}\equiv \mathrm{c\_leak.g\_max}\cdot (\mathrm{l2.v}-\mathrm{c\_leak.v\_eq})$
5. ${\mathrm{l2.v}}^{\prime }\equiv 1000.0\cdot \mathrm{l2.i}/\mathrm{l2.c}$
6. $\mathrm{v\_m}\equiv \mathrm{e\_r}-\mathrm{l2.v}$
7. Within group c_pot (prefix _ indicates shortened variable name)
   1. ${\mathrm{\_gate\_act.n}}^{\prime }\equiv \mathrm{\_gate\_act.phi}\cdot (\mathrm{fopen}(\mathrm{l2.v},-10.0,0.1,100.0)\cdot (1.0-\mathrm{\_gate\_act.n})-\mathrm{gate\_act.fclose}(\mathrm{l2.v},0.0125,125.0)\cdot \mathrm{\_gate\_act.n})$
   2. $\mathrm{\_g}\equiv \mathrm{\_g\_max}\cdot {\mathrm{\_gate\_act.n}}^{4.0}$
   3. $\mathrm{\_i}\equiv \mathrm{\_g}\cdot (\mathrm{l2.v}-\mathrm{\_v\_eq})$
8. Within group c_sod (prefix _ indicates shortened variable name)
   1. ${\mathrm{\_gate\_act.n}}^{\prime }\equiv \mathrm{\_gate\_act.phi}\cdot (\mathrm{fopen}(\mathrm{l2.v},-25.0,0.1,1000.0)\cdot (1.0-\mathrm{\_gate\_act.n})-\mathrm{gate\_act.fclose}(\mathrm{l2.v},0.05555555555555555,4000.0)\cdot \mathrm{\_gate\_act.n})$
   2. ${\mathrm{\_gate\_inact.n}}^{\prime }\equiv \mathrm{\_gate\_inact.phi}\cdot (\mathrm{gate\_act.fclose}(\mathrm{l2.v},0.05,70.0)\cdot (1.0-\mathrm{\_gate\_inact.n})-\mathrm{gate\_inact.fclose}(\mathrm{l2.v},-30.0,-0.1,1000.0)\cdot \mathrm{\_gate\_inact.n})$
   3. $\mathrm{\_g}\equiv \mathrm{\_g\_max}\cdot {\mathrm{\_gate\_act.n}}^{3.0}\cdot \mathrm{\_gate\_inact.n}$
   4. $\mathrm{\_i}\equiv \mathrm{\_g}\cdot (\mathrm{l2.v}-\mathrm{\_v\_eq})$

Functions:

function fopen
  input Real x "membrane potential (as displacement from resting potential)";
  input Real x0 = -10.0 "offset for x (fitting parameter)";
  input Real sx = 0.1 "scaling factor for x (fitting parameter)";
  input Real sy = 100.0 "scaling factor for y (fitting parameter)";
  output Real y "rate of change of the gating variable at given V=x";
  protected Real x_adj "adjusted x with offset and scaling factor";
algorithm
  x_adj := sx * (x - x0);
  if abs(x - x0) < 1e-06 then
    y := sy;
  else
    y := sy * x_adj / (exp(x_adj) - 1.0);
  end if;
end fopen;

function gate_inact.fclose
  input Real x "input value";
  input Real x0 = -30.0 "x-value of sigmoid midpoint (fitting parameter)";
  input Real sx = -0.1 "growth rate/steepness (fitting parameter)";
  input Real y_max = 1000.0 "maximum value";
  output Real y "result";
  protected Real x_adj "adjusted x with offset and scaling factor";
algorithm
  x_adj := sx * (x - x0);
  y := y_max / (exp(-x_adj) + 1.0);
end gate_inact.fclose;

function gate_act.fclose
  input Real x "input value";
  input Real sx = 0.0125 "scaling factor for x axis (fitting parameter)";
  input Real sy = 125.0 "scaling factor for y axis (fitting parameter)";
  output Real y "result";
algorithm
  y := sy * exp(sx * x);
end gate_act.fclose;

Monolithic reference model

model HHmono "monolithic version of the Hodgkin-Huxley model"
  extends PotentialAdapter;
  parameter Real e_r(unit = "mV") = -75 "resting potential";
  Real v_m(unit = "mV") = e_r - v "absolute membrane potential (v_in - v_out)";
  parameter Real Cm(unit = "uF/cm2") = 1;
  parameter Real gbarNa(unit = "mmho/cm2") = 120 "max sodium conductance";
  parameter Real gbarK(unit = "mmho/cm2") = 36 "max potassium conductance";
  parameter Real gbar0(unit = "mmho/cm2") = 0.3;
  parameter Real VNa(unit = "mV") = -115;
  parameter Real VK(unit = "mV") = 12;
  parameter Real Vl(unit = "mV") = -10.613;
  parameter Real Temp = 6.3;
  parameter Real phi = 3 ^ ((Temp - 6.3) / 10);
  parameter Real Vdepolar(unit = "mV") = -90;
  parameter Real Vnorm(unit = "mV") = 1 "for non-dimensionalizing v in function expressions, i.e. exp(v/Vnorm) replaces exp(v).";
  parameter Real msecm1(unit = "1/msec") = 1 "for adding units to alpha and beta variables";
  parameter Real alphan0(unit = "1/msec") = 0.1 / (exp(1) - 1) "always use v=0 to calculate i.c.";
  parameter Real betan0(unit = "1/msec") = 0.125;
  parameter Real alpham0(unit = "1/msec") = 2.5 / (exp(2.5) - 1);
  parameter Real betam0(unit = "1/msec") = 4;
  parameter Real alphah0(unit = "1/msec") = 0.07;
  parameter Real betah0(unit = "1/msec") = 1 / (exp(3) + 1);
  parameter Real minusI(unit = "nA/cm2") = 40;
  input Real Vclamp(unit = "mV");
  parameter Boolean clamp_0no_1yes = false;
  //Variables for all the algebraic equations
  Real INa(unit = "nA/cm2") "Ionic currents";
  Real IK(unit = "nA/cm2") "Ionic currents";
  Real Il(unit = "nA/cm2") "Ionic currents";
  Real alphan(unit = "1/msec") "rate constant of particles from out to in";
  Real betan(unit = "1/msec") "rate constant from in to out";
  Real alpham(unit = "1/msec") "rate constant of activating molecules from out to in";
  Real betam(unit = "1/msec") "rate constant of activating molecules from in to out";
  Real alphah(unit = "1/msec") "rate constant of inactivating molecules from out to in";
  Real betah(unit = "1/msec") "rate constant of inactivating molecules from in to out";
  Real gNa(unit = "mmho/cm2") "Sodium conductance";
  Real gK(unit = "mmho/cm2") "potassium conductance";
  Real Iion(unit = "nA/cm2");
  //State variables for all the ODEs
  Real VV(unit = "mV");
  Real v(unit = "mV") "displacement of the membrane potential from its resting value (depolarization negative)";
  Real n "proportion of the particles in a certain position";
  Real m "proportion of activating molecules on the inside";
  Real h "proportion of inactivating molecules on the outside";
protected
  Real ninf;
  Real minf;
  Real hinf;
  Real taun(unit = "msec");
  Real taum(unit = "msec");
  Real tauh(unit = "msec");
initial equation
  VV = if clamp_0no_1yes then Vclamp else Vdepolar;
  n = alphan0 / (alphan0 + betan0);
  m = alpham0 / (alpham0 + betam0);
  h = alphah0 / (alphah0 + betah0);
equation
//if v/Vnorm == -10 then
//  alphan = 0.1;
//else
  alphan = 0.01 * (v / Vnorm + 10) / (exp((v / Vnorm + 10) / 10) - 1);
//end if;
  betan = msecm1 * (0.125 * exp(v / Vnorm / 80));
//if v/Vnorm == -25 then
//  alpham = 1;
//else
  alpham = 0.1 * (v / Vnorm + 25) / (exp((v / Vnorm + 25) / 10) - 1);
//end if;
  betam = msecm1 * (4 * exp(v / Vnorm / 18));
  alphah = msecm1 * (0.07 * exp(v / Vnorm / 20));
  betah = msecm1 * (1 / (exp((v / Vnorm + 30) / 10) + 1));
  minf = alpham / (alpham + betam);
  ninf = alphan / (alphan + betan);
  hinf = alphah / (alphah + betah);
  taun = 1 / (alphan + betan);
  tauh = 1 / (alphah + betah);
  taum = 1 / (alpham + betam);
  gNa = gbarNa * m ^ 3 * h;
  gK = gbarK * n ^ 4;
  INa = gNa * (v - VNa);
  IK = gK * (v - VK);
  Il = gbar0 * (v - Vl);
  Iion = INa + IK + Il;
  if clamp_0no_1yes then
    der(VV) = 0;
    v = Vclamp;
  else
    der(VV) = ((-minusI) - INa - IK - Il) / Cm;
    v = VV;
  end if;
  der(n) = phi * (alphan * (1 - n) - betan * n);
  der(m) = phi * (alpham * (1 - m) - betam * m);
  der(h) = phi * (alphah * (1 - h) - betah * h);
  annotation(
    experiment(StartTime = 0, StopTime = 30, Tolerance = 1e-6, Interval = 0.01),
    __OpenModelica_simulationFlags(s = "dassl"),
    __MoST_experiment(variableFilter = "v_m|v|gK|gNa|n|m|h"),
    documentation(info = "
    This is a 1:1 translation of the JSim implementation
    available at
    https://www.physiome.org/jsim/models/webmodel/NSR/Hodgkin_Huxley1952/.
  "));
end HHmono;

1. $\mathrm{v\_m}\equiv \mathrm{e\_r}-\mathrm{VV}$
2. $h^{\prime }\equiv \mathrm{phi}\cdot (\mathrm{alphah}\cdot (1.0-h)-\mathrm{betah}\cdot h)$
3. $m^{\prime }\equiv \mathrm{phi}\cdot (\mathrm{alpham}\cdot (1.0-m)-\mathrm{betam}\cdot m)$
4. $n^{\prime }\equiv \mathrm{phi}\cdot (\mathrm{alphan}\cdot (1.0-n)-\mathrm{betan}\cdot n)$
5. ${\mathrm{VV}}^{\prime }\equiv (((-\mathrm{IK}-\mathrm{Il})-\mathrm{INa})-\mathrm{minusI})/\mathrm{Cm}$
6. $\mathrm{Iion}\equiv \mathrm{INa}+\mathrm{IK}+\mathrm{Il}$
7. $\mathrm{Il}\equiv \mathrm{gbar0}\cdot (\mathrm{VV}-\mathrm{Vl})$
8. $\mathrm{IK}\equiv \mathrm{gK}\cdot (\mathrm{VV}-\mathrm{VK})$
9. $\mathrm{INa}\equiv \mathrm{gNa}\cdot (\mathrm{VV}-\mathrm{VNa})$
10. $\mathrm{gK}\equiv \mathrm{gbarK}\cdot n^{4.0}$
11. $\mathrm{gNa}\equiv \mathrm{gbarNa}\cdot m^{3.0}\cdot h$
12. $\mathrm{taum}\equiv 1.0/(\mathrm{alpham}+\mathrm{betam})$
13. $\mathrm{tauh}\equiv 1.0/(\mathrm{alphah}+\mathrm{betah})$
14. $\mathrm{taun}\equiv 1.0/(\mathrm{alphan}+\mathrm{betan})$
15. $\mathrm{hinf}\equiv \mathrm{alphah}/(\mathrm{alphah}+\mathrm{betah})$
16. $\mathrm{ninf}\equiv \mathrm{alphan}/(\mathrm{alphan}+\mathrm{betan})$
17. $\mathrm{minf}\equiv \mathrm{alpham}/(\mathrm{alpham}+\mathrm{betam})$
18. $\mathrm{betah}\equiv \mathrm{msecm1}/(1.0+e^{\mathrm{VV}/\mathrm{Vnorm}\times 10.0+3.0})$
19. $\mathrm{alphah}\equiv 0.07000000000000001\cdot \mathrm{msecm1}\cdot e^{\mathrm{VV}/\mathrm{Vnorm}\times 20.0}$
20. $\mathrm{betam}\equiv 4.0\cdot \mathrm{msecm1}\cdot e^{\mathrm{VV}/\mathrm{Vnorm}\times 18.0}$
21. $\mathrm{alpham}\equiv 0.1\cdot (25.0+\mathrm{VV}/\mathrm{Vnorm})/(-1+e^{\mathrm{VV}/\mathrm{Vnorm}\times 10.0+2.5})$
22. $\mathrm{betan}\equiv 0.125\cdot \mathrm{msecm1}\cdot e^{\mathrm{VV}/\mathrm{Vnorm}\times 80.0}$
23. $\mathrm{alphan}\equiv 0.01\cdot (10.0+\mathrm{VV}/\mathrm{Vnorm})/(-1+e^{\mathrm{VV}/\mathrm{Vnorm}\times 10.0+1.0})$

Functions:

Model structure

The structure of the modular implementation follows both the physiology of the squid giant axon and the electrical analogy that Hodgkin and Huxley used to derive the equations.

Interfaces

The lipid bilayer, the ion channels, and the current clamp used for the experiment protocol each are electrical components with a basic interface.

The voltage-dependent gating molecules of the sodium and potassium channels have an additional interface since they are temperature dependent.

connector TemperatureInput = input Real(unit = "degC") "membrane temperature" annotation(
  Icon(coordinateSystem(preserveAspectRatio = true, extent = {{-100, -100}, {100, 100}}), graphics = {Ellipse(extent = {{-100, 100}, {100, -100}}, lineColor = {255, 0, 0}, fillColor = {255, 255, 255}, fillPattern = FillPattern.Solid)}));

connector TemperatureOutput = output Real(unit = "degC") "membrane temperature" annotation(
  Icon(coordinateSystem(preserveAspectRatio = true, extent = {{-100, -100}, {100, 100}}), graphics = {Ellipse(extent = {{-100, 100}, {100, -100}}, lineColor = {255, 0, 0}, fillColor = {255, 0, 0}, fillPattern = FillPattern.Solid)}));

Base components

The simplest component of the model is the lipid bilayer, which just acts as a capacitor.

model LipidBilayer "lipid bilayer separating external and internal potential (i.e. acting as a capacitor)"
  extends TwoPinComponent;
  extends HHmodelica.Icons.LipidBilayer;
  TemperatureOutput temp = temp_m annotation(
    Placement(transformation(extent = {{40, 48}, {60, 68}})));
  parameter Real temp_m(unit = "degC") = 6.3 "constant membrane temperature";
  parameter Real c(unit = "uF/cm2") = 1 "membrane capacitance";
  parameter Real v_init(unit = "mV") = -90 "short initial stimulation";
initial equation
  v = v_init;
equation
  der(v) = 1000 * i / c "multiply with 1000 to get mV/s instead of v/s";
end LipidBilayer;

The central components that regulate the shape of the action potential curve are the voltage-dependent ion channels. These channels all follow the same basic structure

partial model IonChannel "ionic current through the membrane"
  extends TwoPinComponent;
  extends HHmodelica.Icons.IonChannel;
  Real g(unit = "mmho/cm2") "ion conductance, needs to be defined in subclasses";
  parameter Real v_eq(unit = "mV") "equilibrium potential (as displacement from resting potential)";
  parameter Real g_max(unit = "mmho/cm2") "maximum conductance";
equation
  i = g * (v - v_eq);
end IonChannel;

partial model GatedIonChannel "ion channel that has voltage-dependent gates"
  extends IonChannel;
  TemperatureInput temp "membrane temperature to determine reaction coefficient" annotation(
    Placement(transformation(extent = {{-40, 48}, {-60, 68}})));
end GatedIonChannel;

The specific ion channels that are built upon those base classes are the fast sodium channel, the (rapid delayed rectifier) potassium channel, and a leak channel.

model SodiumChannel "channel selective for Na+ cations"
  extends GatedIonChannel(g_max = 120, v_eq = -115);
  extends HHmodelica.Icons.Activatable;
  extends HHmodelica.Icons.Inactivatable;
  Gate gate_act(redeclare function fopen = goldmanFit(x0 = -25, sy = 1000, sx = 0.1), redeclare function fclose = expFit(sx = 1 / 18, sy = 4000), v = v, temp = temp) "activation gate";
  Gate gate_inact(redeclare function fopen = expFit(sx = 1 / 20, sy = 70), redeclare function fclose = logisticFit(x0 = -30, sx = -0.1, y_max = 1000), v = v, temp = temp) "inactivation gate";
equation
  g = g_max * gate_act.n ^ 3 * gate_inact.n;
end SodiumChannel;

model PotassiumChannel "channel selective for K+ cations"
  extends GatedIonChannel(g_max = 36, v_eq = 12);
  extends HHmodelica.Icons.Activatable;
  Gate gate_act(redeclare function fopen = goldmanFit(x0 = -10, sy = 100, sx = 0.1), redeclare function fclose = expFit(sx = 1 / 80, sy = 125), v = v, temp = temp) "activation gate";
equation
  g = g_max * gate_act.n ^ 4;
end PotassiumChannel;

model LeakChannel "constant leakage current of ions through membrane"
  extends IonChannel(g_max = 0.3, v_eq = -10.613);
  extends HHmodelica.Icons.OpenChannel;
equation
  g = g_max;
end LeakChannel;

The behavior of the sodium and potassium channels is defined by gating molecules that have an open and a closed conformation. All gates in the model follow the same structure.

Fitting functions

Different fitting functions are used to determine the behavior of fclose and fopen in the

function expFit "exponential function with scaling parameters for x and y axis"
  input Real x "input value";
  input Real sx "scaling factor for x axis (fitting parameter)";
  input Real sy "scaling factor for y axis (fitting parameter)";
  output Real y "result";
algorithm
  y := sy * exp(sx * x);
end expFit;

function logisticFit "logistic function with sigmoidal shape"
  input Real x "input value";
  input Real x0 "x-value of sigmoid midpoint (fitting parameter)";
  input Real sx "growth rate/steepness (fitting parameter)";
  input Real y_max "maximum value";
  output Real y "result";
protected
  Real x_adj "adjusted x with offset and scaling factor";
algorithm
  x_adj := sx * (x - x0);
  y := y_max / (exp(-x_adj) + 1);
end logisticFit;

Hodgkin and Huxley state that this formula was (in part) used because it "bears a close resemblance to the equation derived by Goldman (1943) for the movements of a charged particle in a constant field".

We suppose that this statement refers to equation 11 of Goldman (1943), which is also called the Goldman-Hodgkin-Katz flux equation:

j_i = u_i * F / a * dV * (n'_i * exp(-z_i * beta * dV - n0_i)) / exp(-z_i * beta * dV)

Factoring out n0_i from the denominator, substituting n'_i/n0_i = exp(V_0 * beta * z_i) and grouping and renaming variables, the GHK flux equation can be written as

y := sy * sx * x * (exp((x - x0) * sx) - 1) / (exp(x * sx) - 1)

with sx = -z_i * beta, x = dV, x0 = V_0, and sy = n0_i * u_i * F / a * 1 / sx.

With this notation, the similarity becomes apparent, as omitting the denominator (exp((x - x0) * sx) - 1) and using x-x0 instead of x in the rest of the formula gives exactly the goldmanFit used by Hodgkin and Huxley.

function goldmanFit "fitting function related to Goldmans formula for the movement of a charged particle in a constant electrical field"
  input Real x "membrane potential (as displacement from resting potential)";
  input Real x0 "offset for x (fitting parameter)";
  input Real sx "scaling factor for x (fitting parameter)";
  input Real sy "scaling factor for y (fitting parameter)";
  output Real y "rate of change of the gating variable at given V=x";
protected
  Real x_adj "adjusted x with offset and scaling factor";
algorithm
  x_adj := sx * (x - x0);
  if abs(x - x0) < 1e-6 then
    y := sy;
  else
    y := sy * x_adj / (exp(x_adj) - 1);
  end if;
// using L'Hôpital to find limit for x_adj->0
  annotation(
    Documentation(info = "
        <html>
          <p>Hodgkin and Huxley state that this formula was (in part) used
          because it &quot;bears a close resemblance to the equation derived
          by Goldman (1943) for the movements of a charged particle in a constant
          field&quot;.</p>
          <p>We suppose that this statement refers to equation 11 of Goldman
          (1943), which is also called the Goldman-Hodgkin-Katz flux equation:</p>
          <blockquote>
            j_i = u_i * F / a * dV * (n'_i * exp(-z_i * beta * dV - n0_i)) / exp(-z_i * beta * dV)
          </blockquote>
          <p>Factoring out n0_i from the denominator, substituting
          n'_i/n0_i = exp(V_0 * beta * z_i) and grouping and renaming
          variables, the GHK flux equation can be written as</p>
          <blockquote>
            y := sy * sx * x * (exp((x - x0) * sx) - 1) / (exp(x * sx) - 1)
          </blockquote>
          <p>with sx = -z_i * beta, x = dV, x0 = V_0, and
          sy = n0_i * u_i * F / a * 1 / sx.</p>
          <p>With this notation, the similarity becomes apparent, as omitting
          the denominator (exp((x - x0) * sx) - 1) and using x-x0 instead of
          x in the rest of the formula gives exactly the goldmanFit used by
          Hodgkin and Huxley.</p>
        </html>
      "));
end goldmanFit;