start work on mechanism descriptions

data/hh.mod

+TITLE hh.mod   squid sodium, potassium, and leak channels
+
+COMMENT
+ This is the original Hodgkin-Huxley treatment for the set of sodium, 
+  potassium, and leakage channels found in the squid giant axon membrane.
+  ("A quantitative description of membrane current and its application 
+  conduction and excitation in nerve" J.Physiol. (Lond.) 117:500-544 (1952).)
+ Membrane voltage is in absolute mV and has been reversed in polarity
+  from the original HH convention and shifted to reflect a resting potential
+  of -65 mV.
+ Remember to set celsius=6.3 (or whatever) in your HOC file.
+ See squid.hoc for an example of a simulation using this model.
+ SW Jaslove  6 March, 1992
+ENDCOMMENT
+
+UNITS {
+    (mA) = (milliamp)
+    (mV) = (millivolt)
+    (S) = (siemens)
+}
+
+NEURON {
+    SUFFIX hh
+    USEION na READ ena WRITE ina
+    USEION k READ ek WRITE ik
+    NONSPECIFIC_CURRENT il
+    RANGE gnabar, gkbar, gl, el, gna, gk
+    GLOBAL minf, hinf, ninf, mtau, htau, ntau
+    THREADSAFE : assigned GLOBALs will be per thread
+}
+
+PARAMETER {
+    gnabar = .12 (S/cm2)    <0,1e9>
+    gkbar = .036 (S/cm2)    <0,1e9>
+    gl = .0003 (S/cm2)  <0,1e9>
+    el = -54.3 (mV)
+}
+
+STATE {
+    m h n
+}
+
+ASSIGNED {
+    v (mV)
+    celsius (degC)
+    ena (mV)
+    ek (mV)
+
+    gna (S/cm2)
+    gk (S/cm2)
+    ina (mA/cm2)
+    ik (mA/cm2)
+    il (mA/cm2)
+    minf hinf ninf
+    mtau (ms) htau (ms) ntau (ms)
+}
+
+? currents
+BREAKPOINT {
+    SOLVE states METHOD cnexp
+    gna = gnabar*m*m*m*h
+    ina = gna*(v - ena)
+    gk = gkbar*n*n*n*n
+    ik = gk*(v - ek)
+    il = gl*(v - el)
+}
+
+INITIAL {
+    rates(v)
+    m = minf
+    h = hinf
+    n = ninf
+}
+
+DERIVATIVE states
+{
+    rates(v)
+    m' =  (minf-m)/mtau
+    h' = (hinf-h)/htau
+    n' = (ninf-n)/ntau
+}
+
+PROCEDURE rates(v(mV))
+{
+    :Computes rate and other constants at current v.
+    :Call once from HOC to initialize inf at resting v.
+    LOCAL  alpha, beta, sum, q10
+
+    q10 = 3^((celsius - 6.3)/10)
+
+    :"m" sodium activation system
+    alpha = .1 * vtrap(-(v+40),10)
+    beta =  4 * exp(-(v+65)/18)
+    sum = alpha + beta
+    mtau = 1/(q10*sum)
+    minf = alpha/sum
+
+    :"h" sodium inactivation system
+    alpha = .07 * exp(-(v+65)/20)
+    beta = 1 / (exp(-(v+35)/10) + 1)
+    sum = alpha + beta
+    htau = 1/(q10*sum)
+    hinf = alpha/sum
+
+    :"n" potassium activation system
+    alpha = .01*vtrap(-(v+55),10) 
+    beta = .125*exp(-(v+65)/80)
+    sum = alpha + beta
+    ntau = 1/(q10*sum)
+    ninf = alpha/sum
+}
+
+FUNCTION vtrap(x,y) {  :Traps for 0 in denominator of rate eqns.
+    if (fabs(x/y) < 1e-6) {
+        vtrap = y*(1 - x/y/2)
+    }else{
+        vtrap = x/(exp(x/y) - 1)
+    }
+}
+

src/cell.cpp

@@ -2,11 +2,22 @@
 
 namespace nestmc {
 
+cell::cell()
+{
+    // insert a placeholder segment for the soma
+    segments_.push_back(make_segment<placeholder_segment>());
+    parents_.push_back(0);
+}
+
 int cell::num_segments() const
 {
     return segments_.size();
 }
 
+//
+// note: I think that we have to enforce that the soma is the first
+//       segment that is added
+//
 void cell::add_soma(value_type radius, point_type center)
 {
     if(has_soma()) {
@@ -15,24 +26,15 @@ void cell::add_soma(value_type radius, point_type center)
         );
     }
 
-    // soma has intself as its own parent
-    soma_ = num_segments();
-    parents_.push_back(num_segments());
-
     // add segment for the soma
     if(center.is_set()) {
-        segments_.push_back(
-            make_segment<soma_segment>(radius, center)
-        );
+        segments_[0] = make_segment<soma_segment>(radius, center);
     }
     else {
-        segments_.push_back(
-            make_segment<soma_segment>(radius)
-        );
+        segments_[0] = make_segment<soma_segment>(radius);
     }
 }
 
-// add a cable that is provided by the user as a segment_ptr
 void cell::add_cable(cell::index_type parent, segment_ptr&& cable)
 {
     // check for a valid parent id
@@ -52,16 +54,44 @@ void cell::add_cable(cell::index_type parent, segment_ptr&& cable)
     parents_.push_back(parent);
 }
 
+segment* cell::segment(int index)
+{
+    if(index<0 || index>=num_segments()) {
+        throw std::out_of_range(
+            "attempt to access a segment with invalid index"
+        );
+    }
+    return segments_[index].get();
+}
+
+segment const* cell::segment(int index) const
+{
+    if(index<0 || index>=num_segments()) {
+        throw std::out_of_range(
+            "attempt to access a segment with invalid index"
+        );
+    }
+    return segments_[index].get();
+}
+
 
 bool cell::has_soma() const
 {
-    return soma_ > -1;
+    return !segment(0)->is_placeholder();
 }
 
 soma_segment* cell::soma()
 {
     if(has_soma()) {
-        return segments_[soma_].get()->as_soma();
+        return segment(0)->as_soma();
+    }
+    return nullptr;
+}
+
+cable_segment* cell::cable(int index)
+{
+    if(index>0 && index<num_segments()) {
+        return segment(index)->as_cable();
     }
     return nullptr;
 }
@@ -95,7 +125,7 @@ std::vector<segment_ptr> const& cell::segments() const
     return segments_;
 }
 
-void cell::construct()
+void cell::construct() const
 {
     if(num_segments()) {
         tree_ = cell_tree(parents_);
@@ -104,6 +134,7 @@ void cell::construct()
 
 cell_tree const& cell::graph() const
 {
+    construct();
     return tree_;
 }
 

src/cell.hpp

 #pragma once
 
-#include <vector>
+#include <mutex>
 #include <stdexcept>
+#include <thread>
+#include <vector>
 
 #include "segment.hpp"
 #include "cell_tree.hpp"
@@ -17,6 +19,9 @@ namespace nestmc {
         using value_type = double;
         using point_type = point<value_type>;
 
+        // constructor
+        cell();
+
         /// add a soma to the cell
         /// radius must be specified
         void add_soma(value_type radius, point_type center=point_type());
@@ -37,9 +42,18 @@ namespace nestmc {
 
         bool has_soma() const;
 
+        class segment* segment(int index);
+        class segment const* segment(int index) const;
+
         /// access pointer to the soma
+        /// returns nullptr if the cell has no soma
         soma_segment* soma();
 
+        /// access pointer to a cable segment
+        /// will throw an std::out_of_range exception if
+        /// the cable index is not valid
+        cable_segment* cable(int index);
+
         /// the volume of the cell
         value_type volume() const;
 
@@ -48,10 +62,6 @@ namespace nestmc {
 
         std::vector<segment_ptr> const& segments() const;
 
-        /// generate the internal representation of the connectivity
-        /// graph for the cell segments
-        void construct();
-
         /// the connectivity graph for the cell segments
         cell_tree const& graph() const;
 
@@ -61,6 +71,11 @@ namespace nestmc {
 
         private:
 
+        /// generate the internal representation of the connectivity
+        /// graph for the cell segments
+        void construct() const;
+
+
         //
         // the local description of the cell which can be modified by the user
         // in a ad-hoc manner (adding segments, modifying segments, etc)
@@ -70,15 +85,15 @@ namespace nestmc {
         std::vector<index_type> parents_;
         // the segments
         std::vector<segment_ptr> segments_;
-        // index of the soma
-        int soma_ = -1;
 
         //
         // fixed cell description, which is computed from the layout description
-        // above
+        // this computed whenever a call to the graph() method is made
+        // the graph method is const, so tree_ is mutable
         //
+        // note: this isn't remotely thread safe!
 
-        cell_tree tree_;
+        mutable cell_tree tree_;
     };
 
     // create a cable by forwarding cable construction parameters provided by the user

src/segment.hpp

@@ -20,7 +20,8 @@ struct segment_properties {
 enum class segmentKind {
     soma,
     dendrite,
-    axon
+    axon,
+    none
 };
 
 // forward declarations of segment specializations
@@ -68,6 +69,11 @@ class segment {
         return nullptr;
     }
 
+    virtual bool is_placeholder() const
+    {
+        return false;
+    }
+
     segment_properties<value_type> properties;
 
     protected:
@@ -75,6 +81,35 @@ class segment {
     segmentKind kind_;
 };
 
+class placeholder_segment : public segment
+{
+    public:
+
+    using base = segment;
+    using base::kind_;
+    using base::value_type;
+
+    placeholder_segment()
+    {
+        kind_ = segmentKind::none;
+    }
+
+    value_type volume() const override
+    {
+        return std::numeric_limits<value_type>::quiet_NaN();
+    }
+
+    value_type area() const override
+    {
+        return std::numeric_limits<value_type>::quiet_NaN();
+    }
+
+    bool is_placeholder() const override
+    {
+        return true;
+    }
+};
+
 class soma_segment : public segment
 {
     public :

src/stimulus.hpp

+#pragma once
+
+namespace nestmc {
+
+class i_clamp {
+    public:
+
+    using value_type = double;
+
+    i_clamp(value_type del, value_type dur, value_type amp)
+    :   delay_(del),
+        duration_(dur),
+        amplitude_(amp)
+    {}
+
+    value_type delay() const {
+        return delay_;
+    }
+    value_type duration() const {
+        return duration_;
+    }
+    value_type amplitude() const {
+        return amplitude_;
+    }
+
+    void set_delay(value_type d) {
+        delay_ = d;
+    }
+    void set_duration(value_type d) {
+        duration_ = d;
+    }
+    void set_amplitude(value_type a) {
+        amplitude_ = a;
+    }
+
+    // current is set to amplitude for time in the half open interval:
+    //      t \in [delay, delay+duration)
+    value_type amplitude(double t) {
+        if(t>=delay_ && t<(delay_+duration_)) {
+            return amplitude_;
+        }
+        return 0;
+    }
+
+    private:
+
+    value_type delay_     = 0;
+    value_type duration_  = 0;
+    value_type amplitude_ = 0;
+};
+
+} // namespace nestmc

tests/CMakeLists.txt

@@ -16,6 +16,7 @@ set(TEST_SOURCES
     test_cell.cpp
     test_point.cpp
     test_segment.cpp
+    test_stimulus.cpp
     test_swcio.cpp
     test_tree.cpp
     test_run.cpp

tests/test_cell.cpp

@@ -145,7 +145,6 @@ TEST(cell_type, multiple_cables)
         EXPECT_EQ(c.area(), 8. + math::area_sphere(soma_radius));
 
         // construct the graph
-        c.construct();
         auto const& con = c.graph();
 
         EXPECT_EQ(con.num_segments(), 5u);

tests/test_run.cpp

@@ -8,10 +8,32 @@ TEST(run, init)
 
     nestmc::cell cell;
 
-    cell.add_soma(18.8);
-    auto& props = cell.soma()->properties;
+    cell.add_soma(12.6157/2.0);
+    //auto& props = cell.soma()->properties;
 
-    cell.construct();
+    cell.add_cable(0, segmentKind::dendrite, 0.5, 0.5, 200);
 
-    EXPECT_EQ(cell.graph().num_segments(), 1u);
+    EXPECT_EQ(cell.graph().num_segments(), 2u);
+
+    for(auto &s : cell.segments()) {
+        std::cout << "volume : " << s->volume()
+                  << " area : " << s->area()
+                  << " ratio : " << s->volume()/s->area() << std::endl;
+    }
+
+#ifdef example_1
+
+    // in this context (i.e. attached to a segment on a high-level cell)
+    // a mechanism is essentially a set of parameters
+    // - the only "state" is that used to define parameters
+    cell.soma()->add_mechanism("hh");
+
+    auto& soma_hh = cell.soma()->mechanisms("hh");
+    soma_hh.set("gnabar", 0.12);
+    soma_hh.set("gkbar", 0.036);
+    soma_hh.set("gl", 0.0003);
+    soma_hh.set("el", -54.3);
+
+    fvm_cell fv(cell);
+#endif
 }

tests/test_stimulus.cpp

+#include "gtest.h"
+
+#include "../src/stimulus.hpp"
+
+TEST(stimulus, i_clamp)
+{
+    using namespace nestmc;
+
+    // stimulus with delay 2, duration 0.5, amplitude 6.0
+    i_clamp stim(2.0, 0.5, 6.0);
+
+    EXPECT_EQ(stim.delay(), 2.0);
+    EXPECT_EQ(stim.duration(), 0.5);
+    EXPECT_EQ(stim.amplitude(), 6.0);
+
+    // test that current only turned on in the half open interval
+    // t \in [2, 2.5)
+    EXPECT_EQ(stim.amplitude(0.0), 0.0);
+    EXPECT_EQ(stim.amplitude(1.0), 0.0);
+    EXPECT_EQ(stim.amplitude(2.0), 6.0);
+    EXPECT_EQ(stim.amplitude(2.4999), 6.0);
+    EXPECT_EQ(stim.amplitude(2.5), 0.0);
+
+    // update: delay 1.0, duration 1.5, amplitude 3.0
+    stim.set_delay(1.0);
+    stim.set_duration(1.5);
+    stim.set_amplitude(3.0);
+
+    EXPECT_EQ(stim.delay(), 1.0);
+    EXPECT_EQ(stim.duration(), 1.5);
+    EXPECT_EQ(stim.amplitude(), 3.0);
+}
+