Edit on GitHub

Ion channels

Note

This documentation is work in progress. Currently, the extension of Documenter.jl in my package MoST.jl is still experimental. As the package evolves further, this documentation will increase in readability.

Interfaces

InaMo.Currents.Interfaces

The current-defining models in InaMo are compatible to Modelica.Electric and sometimes use other models from this library. This is achived with the interfaces InaMo.Currents.Interfaces.TwoPinCell, InaMo.Currents.Interfaces.TwoPinVertical and InaMo.Currents.Interfaces.OnePortVertical, which are all copies of interfaces from Modelica.Electric.Interfaces, which only change the connector placement.

Membrane components always have a positive pin at the top representing the extracellular potential and a negative pin at the bottom representing the intracellular potential. InaMo.Currents.Interfaces.TwoPinCell places the connectors more freely to visually put them on the outside and inside of the cell in the icon.

InaMo.Currents.Interfaces.TwoPinVertical

partial model TwoPinVertical "Copy of Modelica.Electrical.Interfaes.TwoPin with vertical connector placement"
  SI.Voltage v "Voltage drop of the two pins (= p.v - n.v)";
  Modelica.Electrical.Analog.Interfaces.PositivePin p(v.nominal = 1e-3, i.nominal = 1e-12) "Positive electrical pin" annotation(
    Placement(transformation(extent = {{-10, 110}, {10, 90}})));
  Modelica.Electrical.Analog.Interfaces.NegativePin n(v.nominal = 1e-3, i.nominal = 1e-12) "Negative electrical pin" annotation(
    Placement(transformation(extent = {{-10, -110}, {10, -90}})));
equation
  v = p.v - n.v;
end TwoPinVertical;

InaMo.Currents.Interfaces.TwoPinCell

partial model TwoPinCell "Copy of Modelica.Electrical.Interfaes.TwoPin with adjusted connector placement"
  SI.Voltage v "Voltage drop of the two pins (= p.v - n.v)";
  Modelica.Electrical.Analog.Interfaces.PositivePin p(v.nominal = 1e-3, i.nominal = 1e-12) "Positive electrical pin" annotation(
    Placement(transformation(extent = {{-60, 110}, {-40, 90}})));
  Modelica.Electrical.Analog.Interfaces.NegativePin n(v.nominal = 1e-3, i.nominal = 1e-12) "Negative electrical pin" annotation(
    Placement(transformation(extent = {{-60, -10}, {-40, 10}})));
equation
  v = p.v - n.v;
end TwoPinCell;

InaMo.Currents.Interfaces.OnePortVertical

partial model OnePortVertical "Copy of Modelica.Electrical.Interfaes.OnePort with vertical connector placement"
  SI.Voltage v "Voltage drop of the two pins (= p.v - n.v)";
  SI.Current i "Current flowing from pin p to pin n";
  Modelica.Electrical.Analog.Interfaces.PositivePin p(v.nominal = 1e-3, i.nominal = 1e-12) "Positive electrical pin" annotation(
    Placement(transformation(extent = {{-10, 110}, {10, 90}})));
  Modelica.Electrical.Analog.Interfaces.NegativePin n(v.nominal = 1e-3, i.nominal = 1e-12) "Negative electrical pin" annotation(
    Placement(transformation(extent = {{-10, -110}, {10, -90}})));
equation
  v = p.v - n.v;
  0 = p.i + n.i;
  i = p.i;
end OnePortVertical;

InaMo.Currents.Interfaces.GatedIonChannel

This is the base model for all ion channel components in InaMo.

It defines the inner variables used in other base components and leaves the two variables open_ratio and i_open open for definition in subclasses.

Typically this interface is not used on its own in favor of the more specific interfaces InaMo.Currents.Interfaces.IonChannelElectric and InaMo.Currents.Interfaces.IonChannelGHK.

partial model GatedIonChannel "ion channel with voltage dependent gates"
  extends OnePortVertical;
  extends InaMo.Icons.InsideBottomOutsideTop;
  extends InaMo.Icons.LipidBilayerWithGap;
  inner SI.Voltage v_gate = v "voltage used for activation/inactivation gates";
  inner SI.Current i_ion = i "current used for ion flux";
  Real open_ratio "ratio between 0 (fully closed) and 1 (fully open)";
  Real i_open(nominal = 1e-12) "i if open_ratio = 1";
equation
  i = open_ratio * i_open;
  annotation(
    Documentation(info = "<html>
  <p>This is the base model for all ion channel components in InaMo.</p>
  <p>It defines the <code>inner</code> variables used in other base components
  and leaves the two variables <code>open_ratio</code> and
  <code>i_open</code> open for definition in subclasses.</p>
  <p>Typically this interface is not used on its own in favor of the more
  specific interfaces InaMo.Currents.Interfaces.IonChannelElectric and
  InaMo.Currents.Interfaces.IonChannelGHK.</p>
</html>"));
end GatedIonChannel;

InaMo.Currents.Interfaces.IonChannelElectric

This model is the base class for most ion channels in InaMo.

It provides the electrical equations for i_open and therefore only requires to define open_ratio based on the current state of the gates in the channel.

partial model IonChannelElectric "ion channel based on electrical analog (voltage source + conductor)"
  extends GatedIonChannel;
  parameter SI.ElectricPotential v_eq "equilibrium potential";
  parameter SI.Conductance g_max "maximum conductance";
  SI.Conductance g(nominal = 1e-9) = open_ratio * g_max "ion conductance";
equation
  i_open = g_max * (v - v_eq);
  annotation(
    Documentation(info = "<html>
  <p>This model is the base class for most ion channels in InaMo.</p>
  <p>It provides the electrical equations for <code>i_open</code> and
  therefore only requires to define <code>open_ratio</code> based on the
  current state of the gates in the channel.</p>
</html>"));
end IonChannelElectric;

InaMo.Currents.Interfaces.IonChannelGHK

This base class is an alternative to the usual formulation using InaMo.Currents.Interfaces.IonChannelElectric. It uses the Goldman-Hodgkin-Katz (GHK) flux equation to define the current based on ion concnetrations and membrane permeability.

partial model IonChannelGHK "ion channel with Goldman-Hodgkin-Katz (GHK) behavior"
  extends GatedIonChannel;
  SI.Concentration ion_in "intracellular concentration of ion";
  parameter SI.Concentration ion_ex "extracellular concentration of ion";
  parameter PermeabilityFM ion_p "permeability of ion";
  parameter Integer ion_z = 1 "valence of ion";
  outer parameter SI.Temperature temp "membrane temperature";
equation
  i_open = ghkFlux(v, temp, ion_in, ion_ex, ion_p, ion_z) * unitArea "multiply with unit area to preserve correct units";
  annotation(
    Documentation(info = "<html>
  <p>This base class is an alternative to the usual formulation using
  InaMo.Currents.Interfaces.IonChannelElectric.
  It uses the Goldman-Hodgkin-Katz (GHK) flux equation to define the
  current based on ion concnetrations and membrane permeability.</p>
</html>"));
end IonChannelGHK;

Base classes

InaMo.Currents.Basic.BackgroundChannel

model BackgroundChannel "generic background channel, which is always fully open"
  extends InaMo.Currents.Interfaces.IonChannelElectric;
  extends InaMo.Icons.OpenChannel;
  extends InaMo.Icons.Current(current_name = "I_b");
equation
  open_ratio = 1;
end BackgroundChannel;

InaMo.Currents.Basic.GateAB

This model reprensents a gate in a Hodgkin-Huxley-style ion channel.

To define a specific gate, the fitting functions falpha and fbeta have to be redeclared in the ion channel model where the gate is used. This can be achieved using the redeclare mechanism in Modelica as follows: 

    GateAB activation(
      redeclare function falpha = myFancyFittingFunction(p1=1, p2=2, ...),
      redeclare function fbeta = myFancyFittingFunction(p1=3, p2=4, ...)
    );
    inner SI.ElectricPotential v_gate = myVoltageVariable;

Note that you only have to supply v_gate once, since it is automatically connected to any components which define a outer variable with the same name. If you use the base class InaMo.Currents.Interfaces.GatedIonChannel, you even do not have to supply this variable at all, as it is already defined in the base class.

model GateAB "gating molecule with an open and a closed conformation governed by two functions alpha and beta"
  extends InaMo.Icons.Gate;
  import InaMo.Functions.Fitting.*;
  replaceable function falpha = goldman(x0 = 0, sx = 1, sy = 1) "rate of transfer from closed to open conformation";
  replaceable function fbeta = expFit(x0 = 0, sx = 1, sy = 1) "rate of transfer from open to closed conformation";
  Real n(start = falpha(0) / (falpha(0) + fbeta(0)), fixed = true) "ratio of molecules in open conformation";
  outer SI.ElectricPotential v_gate "membrane potential of enclosing component";
  Real alpha(unit = "1") = falpha(v_gate) "rate of transfer from closed to open conformation";
  Real beta(unit = "1") = fbeta(v_gate) "rate of transfer from open to closed conformation";
  Real steady(unit = "1") = alpha / (alpha + beta) "steady state achieved if current voltage is held constant";
  SI.Duration tau = 1 / (alpha + beta) "time constant for obtaining steady state (i.e. time until difference between n and steady has reduced by a factor of 1/e)";
equation
  der(n) = alpha * (1 - n) - beta * n;
  annotation(
    Documentation(info = "<html>
    <p>This model reprensents a gate in a Hodgkin-Huxley-style ion channel.</p>
    <p>To define a specific gate, the fitting functions <code>falpha</code>
    and <code>fbeta</code> have to be redeclared in the ion channel model where
    the gate is used.
    This can be achieved using the redeclare mechanism in Modelica as follows:
    </p>
    <pre>
    GateAB activation(
      redeclare function falpha = myFancyFittingFunction(p1=1, p2=2, ...),
      redeclare function fbeta = myFancyFittingFunction(p1=3, p2=4, ...)
    );
    inner SI.ElectricPotential v_gate = myVoltageVariable;
    </pre>
    <p>Note that you only have to supply <code>v_gate</code> once, since
    it is automatically connected to any components which define a
    <code>outer</code> variable with the same name.
    If you use the base class InaMo.Currents.Interfaces.GatedIonChannel,
    you even do not have to supply this variable at all, as it is already
    defined in the base class.</p>
  </html>"),
    Icon(graphics = {Text(origin = {-1, -41}, extent = {{-29, 31}, {29, -31}}, textString = "α/β")}));
end GateAB;

InaMo.Currents.Basic.GateTS

This model reprensents a gate in a Hodgkin-Huxley-style ion channel.

To define a specific gate, the fitting functions fsteady and ftau have to be redeclared in the ion channel model where the gate is used. This can be achieved using the redeclare mechanism in Modelica as follows: 

    GateAB activation(
      redeclare function fsteady = myFancyFittingFunction(p1=1, p2=2, ...),
      redeclare function ftau = myFancyFittingFunction(p1=3, p2=4, ...)
    );
    inner SI.ElectricPotential v_gate = myVoltageVariable;

Note that you only have to supply v_gate once, since it is automatically connected to any components which define a outer variable with the same name. If you use the base class InaMo.Currents.Interfaces.GatedIonChannel, you even do not have to supply this variable at all, as it is already defined in the base class.

model GateTS "gating molecule with an open and a closed conformation governed by two functions tau and steady"
  extends InaMo.Icons.Gate;
  import InaMo.Functions.Fitting.*;
  replaceable function ftau = genLogistic "time constant for obtaining steady state (i.e. time until difference between n and steady has reduced by a factor of 1/e)";
  replaceable function fsteady = genLogistic "value that n would reach if v_gate was held constant";
  Real n(start = fsteady(0), fixed = true) "ratio of molecules in open conformation";
  outer SI.ElectricPotential v_gate "membrane potential of enclosing component";
  Real steady = fsteady(v_gate) "value that n would reach if v_gate was held constant";
  SI.Duration tau = ftau(v_gate) "time constant for obtaining steady state (i.e. time until difference between n and steady has reduced by a factor of 1/e)";
equation
  der(n) = (steady - n) / tau;
  annotation(
    Documentation(info = "<html>
    <p>This model reprensents a gate in a Hodgkin-Huxley-style ion channel.</p>
    <p>To define a specific gate, the fitting functions <code>fsteady</code>
    and <code>ftau</code> have to be redeclared in the ion channel model where
    the gate is used.
    This can be achieved using the redeclare mechanism in Modelica as follows:
    </p>
    <pre>
    GateAB activation(
      redeclare function fsteady = myFancyFittingFunction(p1=1, p2=2, ...),
      redeclare function ftau = myFancyFittingFunction(p1=3, p2=4, ...)
    );
    inner SI.ElectricPotential v_gate = myVoltageVariable;
    </pre>
    <p>Note that you only have to supply <code>v_gate</code> once, since
    it is automatically connected to any components which define a
    <code>outer</code> variable with the same name.
    If you use the base class InaMo.Currents.Interfaces.GatedIonChannel,
    you even do not have to supply this variable at all, as it is already
    defined in the base class.</p>
  </html>"),
    Icon(graphics = {Text(origin = {-1, -41}, extent = {{-29, 31}, {29, -31}}, textString = "τ/∞")}));
end GateTS;

InaMo.Currents.Basic.GateTSShift

model GateTSShift "like GateTS but with an additional variable that shifts the steady state curve along the x axis"
  extends GateTS(steady = fsteady(v_gate - shift));
  parameter SI.ElectricPotential shift;
end GateTSShift;

InaMo.Currents.Basic.InstantGate

This model can be used if gating (or similar) voltage-dependent behavior has to be described, but the kinetics are so fast that they can be assumed to instantaneously arrive at their steady state.

model InstantGate "gate with instantaneous behavior"
  extends InaMo.Icons.Gate;
  outer parameter Real FoRT;
  // FIXME: this is only here to acoid a bug in OMC regarding lookup of outer variables
  replaceable function fn = expFit(x0 = 0, sx = 1, sy = 1) "function that defines behavior of gating variable n";
  Real n = fn(v_gate) "ratio of molecules in open conformation";
  outer SI.ElectricPotential v_gate "membrane potential";
  annotation(
    Documentation(info = "<html>
  <p>This model can be used if gating (or similar) voltage-dependent behavior
  has to be described, but the kinetics are so fast that they can be assumed
  to instantaneously arrive at their steady state.</p>
</html>"));
end InstantGate;

Ion Channels

InaMo.Currents.Basic.SodiumChannelBase

model SodiumChannelBase "base model for sodium channels (agnostic about whether na_in is constant or variable)"
  extends InaMo.Currents.Interfaces.IonChannelGHK(ion_in = 0, ion_ex = na_ex, ion_p = na_p, ion_z = 1);
  // ion in must be declared in subclasses
  extends InaMo.Icons.Activatable;
  extends InaMo.Icons.Inactivatable;
  extends InaMo.Icons.Current(current_name = "I_Na");
  outer parameter SI.Concentration na_ex "extracellular sodium concentration";
  outer parameter PermeabilityFM na_p "membrane permeability to Na+ ions";
  // Note: mV -> V by setting x0 /= 1000 and sx *= 1000
  // Note: time scale is already in seconds => no futher changes required
  GateAB act(redeclare function falpha = goldman(x0 = -0.0444, sx = -1000 / 12.673, sy = 460 * 12.673), redeclare function fbeta = expFit(x0 = -0.0444, sx = -1000 / 12.673, sy = 18400)) "voltage-dependent activation";
  // Note: mv -> V by setting x0 /= 1000 and sx *= 1000
  // Note: time scale is already in seconds => no further changes required
  function inact_steady = pseudoABSteady(redeclare function falpha = expFit(x0 = -0.0669, sx = -1000 / 5.57, sy = 44.9), redeclare function fbeta = genLogistic(y_min = 0, y_max = 1491, x0 = -0.0946, sx = 1000 / 12.9, se = 323.3)) "function used for steady state for inactivation gates";
  GateTS inact_fast(redeclare function fsteady = inact_steady, redeclare function ftau = genLogistic(y_min = 0.00035, y_max = 0.03 + 0.00035, x0 = -0.040, sx = -1000 / 6.0, se = 1)) "inactivation gate for fast sodium channels (type1/h1)";
  GateTS inact_slow(redeclare function fsteady = inact_steady, redeclare function ftau = genLogistic(y_min = 0.00295, y_max = 0.12 + 0.00295, x0 = -0.060, sx = -1000 / 2.0, se = 1)) "inactivation gate for slow sodium channels (type2/h2)";
  Real inact_total = 0.635 * inact_fast.n + 0.365 * inact_slow.n "total inactivation resulting from fast and slow inactivation terms";
equation
  open_ratio = act.n ^ 3 * inact_total;
end SodiumChannelBase;

InaMo.Currents.Atrioventricular.SodiumChannel

This model only defines that the intracellular sodium concentration is assumed to be constant. All other channel properties are defined in the base class InaMo.Currents.Basic.SodiumChannelBase.

model SodiumChannel "fast sodium channel (I_Na) as used by Inada 2009"
  extends InaMo.Currents.Basic.SodiumChannelBase(ion_in = na_in);
  outer parameter SI.Concentration na_in "intracellular sodium concentration";
  annotation(
    Documentation(info = "<html>
  <p>This model only defines that the intracellular sodium concentration is
  assumed to be constant.
  All other channel properties are defined in the base class
  InaMo.Currents.Basic.SodiumChannelBase.</p>
</html>"));
end SodiumChannel;

InaMo.Currents.Atrioventricular.RapidDelayedRectifierChannel

model RapidDelayedRectifierChannel "Rapid delayed rectifier potassium channel (I_K,r)"
  extends InaMo.Currents.Interfaces.IonChannelElectric(g_max = 1.5e-9);
  extends InaMo.Icons.Activatable;
  extends InaMo.Icons.Inactivatable;
  extends InaMo.Icons.Current(current_name = "I_K,r");
  function act_steady = genLogistic(x0 = -10.22e-3, sx = 1000 / 8.5) "function for steady state of activation gates";
  GateTS act_fast(redeclare function ftau = pseudoABTau(redeclare function falpha = expFit(sx = 0.0398e3, sy = 17), redeclare function fbeta = expFit(sx = -0.051e3, sy = 0.211)), redeclare function fsteady = act_steady) "voltage-dependent fast activation";
  GateTS act_slow(redeclare function ftau = gaussianAmp(y_min = 0.33581, y_max = 0.90673 + 0.33581, x0 = -10e-3, sigma = sqrt(988.05 / 2) / 1000), redeclare function fsteady = act_steady) "voltage-dependent slow activation";
  GateTS inact(redeclare function ftau = pseudoABTau(redeclare function falpha = expFit(sx = 0.00942e3, sy = 603.6), redeclare function fbeta = expFit(sx = -0.0183e3, sy = 92.01)), redeclare function fsteady = fprod(redeclare function fa = genLogistic(x0 = -4.9e-3, sx = -1000 / 15.14), redeclare function fb = gaussianAmp(y_min = 1, y_max = (-0.3) + 1, sigma = sqrt(500 / 2) / 1000))) "voltage-dependent inactivation";
  Real act_total = 0.9 * act_fast.n + 0.1 * act_slow.n "total activation due to slow and fast activation terms";
equation
  open_ratio = act_total * inact.n;
end RapidDelayedRectifierChannel;

InaMo.Currents.Atrioventricular.LTypeCalciumChannel

Additional to the equations in Inada 2009, this model also includes an acetylcholine-dependent term found in the C++ code, which inhibits the current up to a percentage given by k_ach for high acetylcholine concentrations. However, the parameter settings for this term cannot be found in the C++ code. We therefore assume that the term was not used in the simulations presented in Inada 2009.

model LTypeCalciumChannel "L-type calcium channel (I_Ca,L)"
  extends InaMo.Currents.Interfaces.IonChannelElectric(g_max = 18.5e-9, v_eq = 62.1e-3);
  extends InaMo.Concentrations.Interfaces.TransmembraneCaFlow(n_ca = 1);
  extends InaMo.Icons.Activatable;
  extends InaMo.Icons.Inactivatable;
  extends InaMo.Icons.Current(current_name = "I_Ca,L");
  outer parameter Boolean use_ach "model ACh dependence or not";
  parameter Real k_ach(unit = "1", min = 0, max = 1) = 0 "ratio of maximum channel inhibition by acetylcholine";
  outer parameter SI.Concentration ach "extracellular(?) acetylcholine concentration";
  GateTS act(redeclare function ftau = pseudoABTau(redeclare function falpha = fsum(redeclare function fa = goldman(x0 = -35e-3, sx = -1000 / 2.5, sy = 26.12 * 2.5), redeclare function fb = goldman(sx = -0.208e3, sy = 78.11 / 0.208)), redeclare function fbeta = goldman(x0 = 5e-3, sx = 0.4e3, sy = 10.52 / 0.4)), redeclare function fsteady = genLogistic(x0 = -3.2e-3, sx = 1000 / 6.61)) "voltage-dependent activation";
  // parameters for AN node
  function inact_steady = genLogistic(x0 = -29e-3, sx = -1000 / 6.31);
  GateTS inact_slow(redeclare function ftau = gaussianAmp(y_min = 0.06, y_max = 1.08171 + 0.06, x0 = -40e-3, sigma = sqrt(138.04 / 2) / 1000), redeclare function fsteady = inact_steady) "voltage-dependent slow inactivation";
  GateTS inact_fast(redeclare function ftau = gaussianAmp(y_min = 0.01, y_max = 0.1539 + 0.01, x0 = -40e-3, sigma = sqrt(185.67 / 2) / 1000), redeclare function fsteady = inact_steady) "voltage-dependent fast inactivation";
  Real inact_total = 0.675 * inact_fast.n + 0.325 * inact_slow.n;
protected
  // TODO check order of magnitude for MM-constant
  Real ach_factor = if use_ach then 1 - k_ach * michaelisMenten(ach, 0.9e-4) else 1;
equation
  open_ratio = act.n * inact_total * ach_factor;
  annotation(
    Documentation(info = "<html>
  <p>Additional to the equations in Inada 2009, this model also includes
  an acetylcholine-dependent term found in the C++ code, which inhibits the
  current up to a percentage given by <code>k_ach</code> for high acetylcholine
  concentrations.
  However, the parameter settings for this term cannot be found in the C++
  code.
  We therefore assume that the term was not used in the simulations presented
  in Inada 2009.</p>
</html>"));
end LTypeCalciumChannel;

InaMo.Currents.Atrioventricular.LTypeCalciumChannelN

This is a slight variation of I_Ca,L used for N-type cells. It only exchanges the equation for the steady state of the activation gate. 

model LTypeCalciumChannelN "L-type calcium channel (I_Ca,L) for N-type cells"
  extends LTypeCalciumChannel(act(redeclare function fsteady = genLogistic(x0 = -18.2e-3, sx = 1000 / 5)));
  annotation(
    Documentation(info = "<html>
  <p>This is a slight variation of I_Ca,L used for N-type cells.
  It only exchanges the equation for the steady state of the activation gate.
  </p>
</html>"));
end LTypeCalciumChannelN;

InaMo.Currents.Atrioventricular.InwardRectifier

This model is an extension of the inward rectifier in Lindblad 1996 by Inada et al.. It includes an additional voltage-dependent activation term voltage_act, which decreases the current to 50% for low voltages and then increases it back to its full value for v > 30 mV, which gives the I-V curve a dual mode shape. We are not sure if this change is sensible, but we include it since it is both present in the article and the C++ code of Inada 2009.

model InwardRectifier "inward rectifying potassium channel (I_K1)"
  extends InaMo.Currents.Interfaces.IonChannelElectric(g_max = 12.5e-9, v_eq = -81.9e-3);
  extends InaMo.Icons.Activatable;
  extends InaMo.Icons.Inactivatable;
  extends InaMo.Icons.Current(current_name = "I_K1");
  outer parameter SI.Temperature temp "membrane temperature of enclosing component";
  // FIXME: FoRT should not need the "inner" keyword, this is just done to
  // help OMC with name lookup of variables due to a bug
  inner parameter Real FoRT = Modelica.Constants.F / temp / Modelica.Constants.R "helper variable to simplyfiy equations";
  parameter Boolean use_vact = true "use voltage-dependent activation gate? (only Inada 2009)";
  outer parameter SI.Concentration k_ex "extracellular potassium concentration";
  Real n_pot = michaelisMenten(k_ex, 0.59) "[K+]_ex-dependent gating variable";
  // Note: R in mJ/(mol * K) -> R in J/(mol * K) by setting sx /= 1000
  // Note: mv -> V by setting x0 /= 1000 and sx *= 1000
  // Note: sx /= 1000 and sx *= 1000 cancel each other out => no change in sx
  InstantGate voltage_inact(redeclare function fn = genLogistic(x0 = v_eq - 3.6e-3, sx = -1.393 * FoRT)) "voltage-dependent inactivation (Lindblad 1996)";
  // Note: mv -> V by setting x0 /= 1000 and sx *= 1000
  InstantGate voltage_act(redeclare function fn = genLogistic(y_min = 0.5, x0 = -30e-3, sx = 1000 / 5)) "voltage-dependent activation (only Inada 2009)";
equation
  if use_vact then
    open_ratio = n_pot ^ 3 * voltage_inact.n * voltage_act.n;
  else
    open_ratio = n_pot ^ 3 * voltage_inact.n;
  end if;
  annotation(
    Documentation(info = "<html>
  <p>This model is an extension of the inward rectifier in Lindblad 1996 by
  Inada et al..
  It includes an additional voltage-dependent activation term
  <code>voltage_act</code>, which decreases the current to 50% for low
  voltages and then increases it back to its full value for v > 30 mV, which
  gives the I-V curve a dual mode shape.
  We are not sure if this change is sensible, but we include it since it is
  both present in the article and the C++ code of Inada 2009.</p>
</html>"));
end InwardRectifier;

InaMo.Currents.Atrioventricular.TransientOutwardChannel

This model implements the equations for I_to in Inada 2009.

NOTE: inact_fast.ftau.y_min is given as 0.1266 (126.6 ms) in Inada 2009, but this is inconsistent with plot S2C in the same article. We therefore suspect that this is an order of magnitude error and the value should be 0.01266 (12.66 ms) instead, which is also consistent with the C++ code.

NOTE: inact_slow.ftau.y_min is given as 0.1 s and this value is also used in the C++ code. However, Figure S2D in Inada 2009 clearly shows a value around 0.2 s. We stick with 0.1 s, because we assume that Figure S2D may just use an old setting and we hope that the C++ code is the most accurate source with regard to actual simulation behavior.

model TransientOutwardChannel "transient outward potassium current (I_to)"
  extends InaMo.Currents.Interfaces.IonChannelElectric(g_max = 20e-9);
  extends InaMo.Icons.Activatable;
  extends InaMo.Icons.Inactivatable;
  extends InaMo.Icons.Current(current_name = "I_to");
  // v_eq ~= -0.08696 V (calculated with nernst)
  GateTS act(redeclare function ftau = pseudoABTau(redeclare function falpha = expFit(x0 = -30.61e-3, sx = 0.09e3, sy = 1.037 / 3.188e-3), redeclare function fbeta = expFit(x0 = -23.84e-3, sx = -0.12e3, sy = 0.396 / 3.188e-3), off = 0.596e-3), redeclare function fsteady = genLogistic(x0 = 7.44e-3, sx = 1000 / 16.4)) "voltage-dependent activation";
  function inact_steady = genLogistic(x0 = -33.8e-3, sx = -1000 / 6.12);
  GateTS inact_slow(redeclare function ftau = gaussianAmp(y_min = 0.1, y_max = 4 + 0.1, x0 = -65e-3, sigma = sqrt(500 / 2) / 1000), redeclare function fsteady = inact_steady) "voltage-dependent slow inactivation";
  GateTS inact_fast(redeclare function ftau = genLogistic(y_min = 0.01266, y_max = 4.72716 + 0.01266, x0 = -154.5e-3, sx = -1000 / 23.96), redeclare function fsteady = inact_steady) "voltage-dependent fast inactivation";
  Real inact_total = 0.55 * inact_slow.n + 0.45 * inact_fast.n "total inactivation resulting from fast and slow inactivation gates";
equation
  open_ratio = act.n * inact_total;
  annotation(
    Documentation(info = "<html>
  <p>This model implements the equations for I_to in Inada 2009.</p>
  <p>NOTE: inact_fast.ftau.y_min is given as 0.1266 (126.6 ms) in Inada 2009,
  but this is inconsistent with plot S2C in the same article.
  We therefore suspect that this is an order of magnitude error and the
  value should be 0.01266 (12.66 ms) instead, which is also consistent with
  the C++ code.</p>
  <p>NOTE: inact_slow.ftau.y_min is given as 0.1 s and this value is also used
  in the C++ code.
  However, Figure S2D in Inada 2009 clearly shows a value around 0.2 s.
  We stick with 0.1 s, because we assume that Figure S2D may just use an old
  setting and we hope that the C++ code is the most accurate source with
  regard to actual simulation behavior.</p>
</html>"));
end TransientOutwardChannel;

InaMo.Currents.Atrioventricular.HyperpolarizationActivatedChannel

model HyperpolarizationActivatedChannel "hyperpolarization-activated channel responsible for \"funny\" current (I_f, HCN4)"
  extends InaMo.Currents.Interfaces.IonChannelElectric(g_max = 1e-9, v_eq = -30e-3);
  extends InaMo.Icons.Activatable;
  extends InaMo.Icons.Current(current_name = "I_f");
  outer parameter Boolean use_ach;
  outer parameter SI.Concentration ach;
  parameter Real act_shift = if use_ach then -7.2 * hillLangmuir(ach, 1.26e-5, 0.69) else 0;
  GateTSShift act(redeclare function ftau = gaussianAmp(y_min = 0.25, y_max = 2 + 0.25, x0 = -70e-3, sigma = sqrt(500 / 2) / 1000), redeclare function fsteady = genLogistic(x0 = -83.19e-3, sx = -1000 / 13.56), shift = act_shift);
equation
  open_ratio = act.n;
end HyperpolarizationActivatedChannel;

InaMo.Currents.Atrioventricular.SustainedInwardChannel

This model was originally developed by Kurata et al. (Kurata 2002). It was slightly changed by Inada et al. by adjusting parameter settings and the equation for act.fsteady.

NOTE: v_eq is not given in Inada 2009 (where it is called E_st). We therefore assume the value given in Kurata 2002. The CellML implementation of Inada 2009 uses -37.4 mV instead of 37.4 mV, but this seems to be an error.

NOTE: Since Kurata 2002 gives tau in ms, but inada uses s as unit, we needed to multiply both alpha and beta for act and inact by a factor of 1000.

NOTE: Inada 2009 uses act.ftau.falpha.fbeta.sx = 1000/700, but Kurata 2002 uses -1000/700, which is more plausible. We therefore use the value from Kurata 2002.

model SustainedInwardChannel "sustained inward current (I_st)"
  extends InaMo.Currents.Interfaces.IonChannelElectric(g_max = 0.1e-9, v_eq = 37.4e-3);
  extends InaMo.Icons.Activatable;
  extends InaMo.Icons.Inactivatable;
  extends InaMo.Icons.Current(current_name = "I_st");
  GateTS act(redeclare function ftau = pseudoABTau(redeclare function falpha = pseudoABTau(redeclare function falpha = expFit(sy = 0.15e-3, sx = -1000 / 11), redeclare function fbeta = expFit(sy = 0.2e-3, sx = -1000 / 700)), redeclare function fbeta = pseudoABTau(redeclare function falpha = expFit(sy = 16e-3, sx = 1000 / 8), redeclare function fbeta = expFit(sy = 15e-3, sx = 1000 / 50))), redeclare function fsteady = genLogistic(x0 = -49.1e-3, sx = 1000 / 8.98)) "voltage-dependent activation";
  GateAB inact(redeclare function falpha = pseudoABTau(redeclare function falpha = expFit(sy = 3100 / 0.1504e3, sx = 1000 / 13), redeclare function fbeta = expFit(sy = 700 / 0.1504e3, sx = 1000 / 70)), redeclare function fbeta = fsum(redeclare function fa = pseudoABTau(redeclare function falpha = expFit(sy = 95 / 0.1504e3, sx = -1000 / 10), redeclare function fbeta = expFit(sy = 50 / 0.1504e3, sx = -1000 / 700)), redeclare function fb = genLogistic(y_max = 0.000229e3, sx = 1000 / 5))) "voltage-dependent inactivation";
equation
  open_ratio = act.n * inact.n;
  annotation(
    Documentation(info = "<html>
  <p>This model was originally developed by Kurata et al. (Kurata 2002).
  It was slightly changed by Inada et al. by adjusting parameter settings and
  the equation for act.fsteady.</p>
  <p>NOTE: v_eq is not given in Inada 2009 (where it is called E_st).
  We therefore assume the value given in Kurata 2002.
  The CellML implementation of Inada 2009 uses -37.4 mV instead of 37.4 mV,
  but this seems to be an error.</p>
  <p>NOTE: Since Kurata 2002 gives tau in ms, but inada uses s as unit, we
  needed to multiply both alpha and beta for act and inact by a factor of
  1000.</p>
  <p>NOTE: Inada 2009 uses act.ftau.falpha.fbeta.sx = 1000/700, but Kurata
  2002 uses -1000/700, which is more plausible.
  We therefore use the value from Kurata 2002.</p>
</html>"));
end SustainedInwardChannel;

InaMo.Currents.Atrioventricular.AcetylcholineSensitiveChannel

This channel formulation is found in the C++ implementation by Inada et al., but not in the article. We are not aware of any scientific publication that describes I_ACh with these precise equations and therefore there is also no reference plot available.

We can only speculate that Inada et al. probably experimented with including this current but ultimately chose not to. Our models therefore also by default disable this current.

model AcetylcholineSensitiveChannel "acetylcholine sensitive potassium channel (I_ACh) as found in the C++ implementation of Inada 2009"
  extends InaMo.Currents.Interfaces.IonChannelElectric(g_max = g_ach * g_k);
  extends InaMo.Icons.Activatable;
  extends InaMo.Icons.Inactivatable;
  extends InaMo.Icons.Current(current_name = "I_ACh");

  function constValue
    input Real x;
    input Real c = 0;
    output Real y;
  algorithm
    y := c;
  end constValue;

  parameter SI.Concentration k_ach = 3.5e-7;
  outer parameter SI.Temperature temp;
  outer parameter SI.Concentration k_ex;
  outer parameter SI.Concentration ach;
  parameter Real g_ach = hillLangmuir(ach, k_ach, 1.5);
  parameter Real g_k = michaelisMenten(k_ex, 10);
  inner parameter Real FoRT = Modelica.Constants.F / (Modelica.Constants.R * temp);
  GateAB inact_fast(redeclare function falpha = constValue(c = 73.1), redeclare function fbeta = genLogistic(y_max = 120, x0 = -50e-3, sx = 1000 / 15));
  GateAB inact_slow(redeclare function falpha = constValue(c = 3.7), redeclare function fbeta = genLogistic(y_max = 5.82, x0 = -50e-3, sx = 1000 / 15));
  InstantGate act(redeclare function fn = genLogistic(x0 = v_eq + 140e-3, sx = FoRT / 2.5));
equation
  open_ratio = act.n * inact_slow.n * inact_fast.n;
  annotation(
    Documentation(info = "<html>
  <p>This channel formulation is found in the C++ implementation by
  Inada et al., but not in the article.
  We are not aware of any scientific publication that describes I_ACh with
  these precise equations and therefore there is also no reference plot
  available.</p>
  <p>We can only speculate that Inada et al. probably experimented with
  including this current but ultimately chose not to.
  Our models therefore also by default disable this current.</p>
</html>"));
end AcetylcholineSensitiveChannel;