Add some example mechanisms to language proposal

doc/specification/examples/ca_dynamics.al

@@ -0,0 +1,12 @@
+# Adapted from Allen Institute Ca_dynamics.mod, in turn based on model of Destexhe et al. 1994.
+
+interface concentration "CaDynamics" {
+    export parameter gamma = 0.05;    # Proportion of unbuffered calcium.
+    export parameter decay = 80 ms;   # Calcium removal time constant.
+    export parameter minCai = 1e-4 mM;
+    export parameter depth = 0.1 µm;  # Depth of shell.
+
+    bind flux = molar flux "ca";
+    bind cai = internal concentration "ca";
+    effect molar flux rate "ca" = -gamma*flux - (cai - minCai)/decay;
+}

doc/specification/examples/expsyn_stdp.al

@@ -0,0 +1,64 @@
+# Based on Arbor default/expsyn_stdp.mod
+
+interface discrete "expsyn_stdp" {
+    # A scaling factor for incoming (pre-synaptic) events is required, as the
+    # weight of an event is dimensionelss.
+
+    export parameter A     =  1 μS;    # pre-synaptic event contribution factor
+    export parameter Apre  =  0.01 μS; # pre-synaptic event plasticity contribution
+    export parameter Apost = -0.01 μS; # post-synaptic event plasticity contribution
+
+    export parameter τ      = 2 ms; # synaptic time constant
+    export parameter τpre  = 10 ms; # pre-synaptic plasticity contribution time constant
+    export parameter τpost = 10 ms; # pre-synaptic plasticity contribution time constant
+
+    export parameter gmax  = 10 μS; # maximum synaptic conductance
+    export parameter e = 0 mV;      # reversal potential
+
+    initial state = {
+        g = 0 μS;
+        apre = 0 μS;
+        apost = 0 μS;
+        w_plastic = 0 μS;
+    };
+
+    # Expression below could have been written without the qualifying 'S.'
+    # by using a 'with S;' or 'with state;' instead.
+
+    bind S = state;
+
+    evolve state' = {
+        g' = -S.g/τ;
+        apre' = -S.apre/τpre;
+        apost'= -S.apost/τpost;
+    };
+
+    bind v = membrane potential;
+    effect current = S.g*(v - e);
+
+    on w = event; state = {
+        # With proposed clamp syntax, this could be:
+        # g = (S.g + S.w_plastic + w·A) ↓ [0 S, gmax];
+
+        g = let g = S.g + S.w_plastic + w·A;
+            | g < 0 S   → 0 S
+            | g > gmax  → gmax
+            | otherwise → g;
+
+        apre = S.apre + Apre;
+        w_plastic = S.w_platic + S.apost;
+    };
+
+    on δt = post; state = {
+        apost = S.apost + Apost;
+        w_plastic = S.w_plastic + S.apre;
+    };
+
+    # The 'δt' above is ignored; it could be incorporated for an update that accounts for the
+    # delay between the post-synaptic event and its triggering of the `on` clause. e.g.
+    #
+    # on δt = post; state = {
+    #     apost = S.apost + Apost·exp(-δt/τpost);
+    #     w_plastic = S.w_plastic + S.apre·exp(δt/τpre);
+    # };
+}

doc/specification/examples/kd.al

@@ -0,0 +1,79 @@
+# Adapted from Allen Institute Kd.mod, in turn based on Kd model of Foust et al. 2011.
+#
+# Three versions, one which takes the reversal potential as a parameter, and one where it
+# is computed dynamically from ion concentrations.
+
+module Kd {
+    parameter gbar = 10^-5 S/cm^2;
+
+    def mInf = fn (v: voltage) ->
+        1 - 1/(1 + exp((v + 43 mV)/8 mV));
+
+    def hInf = fn (v: voltage) ->
+        1/(1 + exp((v + 67 mV)/7.3 mV));
+
+    type state = {
+        m: real;
+        h: real;
+    };
+
+    def state0 = fn (v: voltage) ->
+        {
+            m = mInf(v);
+            h = hInf(v);
+        };
+
+    def rate = fn (s: state, v: voltage) ->
+        with s; {
+            m' = (m - mInf(v))/1 ms;
+            h' = (h - hInf(v))/1500 ms;
+        };
+
+    def current = fn (s: state, v_minus_ek: voltage) ->
+        with s;
+        gbar*m*h*v_minus_ek: current/area;
+}
+
+interface density "Kd" {
+    import Kd;
+    export density parameter Kd.gbar as gbar;
+    export parameter ek = -77 mV;
+
+    bind v = membrane potential;
+
+    initial state = Kd.state0(v);
+    evolve state' = Kd.rate(state, v);
+    effect current density "k" = Kd.current(state, v - ek);
+}
+
+interface density "Kd_nernst" {
+    import Kd;
+    export density parameter Kd.gbar as gbar;
+
+    bind v = membrane potential;
+    bind ek = nernst potential "k";
+
+    initial state = Kd.state0(v);
+    evolve state' = Kd.rate(state, v);
+    effect current density "k" = Kd.current(state, v - ek);
+}
+
+# Kd_nernst is equivalent to:
+#
+# interface density "Kd_nernst" {
+#     import Kd;
+#     export density parameter Kd.gbar as gbar;
+#
+#     bind v = membrane potential;
+#     bind T = temperature;
+#     bind ki = internal concentration "k";
+#     bind ko = external concentration "k";
+#     bind kz = charge "k";
+#
+#     initial state = Kd.state0(v);
+#     evolve state' = Kd.rate(state, v);
+#     effect current density "k" =
+#         let ek = nernst(kz, T, ki, ko);
+#         Kd.current(state, v - ek);
+# }
+

doc/specification/examples/na_ts.al

@@ -0,0 +1,56 @@
+# Adapted from Allen Institute NaTs.mod, in turn based on model of Colbert and Pan 2002.
+
+interface density "NaTs" {
+    type rate = real/time;
+
+    export density parameter gbar    = 10^-5 S/cm^2;
+
+    export parameter mαF: rate       = 0.182 ms⁻¹;
+    export parameter mβF: rate       = 0.124 ms⁻¹;
+    export parameter mv½: voltage    = -40 mV;
+    export parameter mk: voltage     = 6 mV;
+
+    export parameter hαF: rate       = 0.015 ms⁻¹;
+    export parameter hβF: rate       = 0.015 ms⁻¹;
+    export parameter hv½: voltage    = -66 mV;
+    export parameter hk: voltage     = 6 mV;
+
+    def rates = fn (v: voltage, T: temperature) →
+        let qt = 2.3^((T-296.15 K)/10 K);
+
+        let mα = mαF·mk/exprel((mv½-v)/mk);
+        let mβ = mβF·mk/exprel((v-mv½)/mk);
+
+        let hα = hαF·hk/exprel((v-hv½)/hk);
+        let hβ = hβF·hk/exprel((hv½-v)/hk);
+
+        {
+            mi = mα/(mα + mβ);
+            mτ = 1/(mα + mβ)/qt;
+
+            hi = hα/(hα + hβ);
+            hτ = 1/(hα + hβ)/qt;
+        };
+
+    bind T = temperature;
+    bind v = membrane potential;
+    bind ena = nernst potential "na";
+
+    initial state =
+        with rates(v, T); {
+            m = mi;
+            h = hi;
+        };
+
+    evolve state' =
+        with state;
+        with rates(v, T); {
+            m' = (mi - m)/mτ;
+            h' = (hi - h)/hτ;
+        };
+
+    effect current density "na" =
+        with state;
+        gbar·m³·h·(v - ena);
+}
+

doc/specification/examples/nav.al

@@ -0,0 +1,85 @@
+# Adapted from Allen Institute NaV.mod, based on kinetics
+# from Carter et al. 2012 doi:10.1016/j.neuron.2012.08.033 (see Figure 7A).
+
+interface density "NaV" {
+    export density parameter g̅ = 0.015 S/cm²;
+
+    export parameter Con  =  0.01 ms⁻¹; # closed C1 → inactivated I1 transition
+    export parameter Coff = 40    ms⁻¹; # inactivated I1 → closed C1 transitions
+    export parameter Oon  =  8    ms⁻¹; # open O → inactivated I6 transition
+    export parameter Ooff =  0.05 ms⁻¹; # inactivated I6 → open O transition
+    export parameter α =   400    ms⁻¹; # closed Cx right transitions (activation)
+    export parameter β =    12    ms⁻¹; # closed Cx left transitions (deactivation)
+    export parameter γ =   250    ms⁻¹; # closed → open transition
+    export parameter δ =    60    ms⁻¹; # open → closed transition
+
+    export parameter a = 2.51;          # factor for right Ix transitions
+    export parameter b = 5.32;          # inverse factor for left Ix transitions
+
+    export parameter αvdep =  24 mV;    # Vdep of activation
+    export parameter βvdep = -24 mV;    # Vdep of deactivation
+
+    # With the proposed range-limited extension, fields below would be defined
+    #     C1: [0, 1];
+    # etc.
+    type state = {
+        C1: real;
+        C2: real;
+        C3: real;
+        C4: real;
+        C5: real;
+        O:  real;
+        I1: real;
+        I2: real;
+        I3: real;
+        I4: real;
+        I5: real;
+        I6: real;
+    };
+
+    def kinetics = fn(s: state, v: voltage, Q: real) →
+        # scale rates by Q and voltage dependencies
+        let Con  = Q·Con;
+        let Coff = Q·Coff;
+        let Oon  = Q·Oon;
+        let Ooff = Q·Ooff;
+        let α    = Q·α·exp(v/αvdep);
+        let β    = Q·β·exp(v/βvdep);
+        let γ    = Q·γ;
+        let δ    = Q·δ;
+
+        with s; {
+            C1 ⇄ C2  ( 4·α,   β );
+            C2 ⇄ C3  ( 3·α, 2·β );
+            C3 ⇄ C4  ( 2·α, 3·β );
+            C4 ⇄ C5  (   α, 4·β );
+            C5 ⇄ O   (   γ,   δ );
+
+            I1 ⇄ I2  ( 4·a·α,   β/b );
+            I2 ⇄ I2  ( 3·a·α, 2·β/b );
+            I3 ⇄ I4  ( 2·a·α, 3·β/b );
+            I4 ⇄ I5  (   a·α, 4·β/b );
+            I5 ⇄ I6  (     γ,   δ   );
+
+            C1 ⇄ I1  (     Con, Coff    );
+            C2 ⇄ I2  (   a·Con, Coff/b  );
+            C3 ⇄ I3  (  a²·Con, Coff/b² );
+            C4 ⇄ I4  (  a³·Con, Coff/b³ );
+            C5 ⇄ I5  (  a⁴·Con, Coff/b⁴ );
+            O  ⇄ I6  (     Oon, Ooff    );
+        }: state'
+
+    bind v = membrane potential;
+    bind T = temperature;
+    bind ena = nernst potential "na";
+
+    def qt(T: temperature) → 2.3^((T - 310.15 K)/10 K);
+
+    initial steady state from state =
+        { C1 = 1; C2 = 0; C3 = 0; C4 = 0; C5 = 0; O = 0;
+          I1 = 0; I2 = 0; I3 = 0; I4 = 0; I5 = 0; I6 = 0; };
+
+    evolve state' = kinetics(state, v, qt(T));
+
+    effect current density "na" = with state; g̅·O·(v - ena);
+}

doc/specification/examples/sk.al

@@ -0,0 +1,18 @@
+# Adapted from Allen Institute SK.mod, in turn based on model of Kohler et al. 1996.
+# Takes fixed ek as parameter.
+
+interface density "SK" {
+    export density parameter gbar = 10^-5 S/cm^2;
+    export parameter zTau = 1 ms;
+
+    def zInf(v: voltage; c: concentration) ->
+        | c==0 -> 0
+        | otherwise -> 1/(1 + (0.00043 mM/c)^4.8);
+
+    bind v = voltage;
+    bind cai = internal concentration "ca";
+
+    initial state = zInf(v, cai);
+    evolve state' = (zInf(v, cai) - state)/zTau;
+    effect current density "k" = gbar*state*(v - ek);
+}