diff --git a/data/hh.mod b/data/hh.mod
new file mode 100755
index 0000000000000000000000000000000000000000..6020be0c1f6d99ab90a3639916eb32cdb1fab155
--- /dev/null
+++ b/data/hh.mod
@@ -0,0 +1,120 @@
+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)
+ }
+}
+
diff --git a/src/cell.cpp b/src/cell.cpp
index 5d9cec1b22e43d265f89ec3998f3c33e566f6ebc..c0d51bc32741c17b73da8215d7076b5325b21cb9 100644
--- a/src/cell.cpp
+++ b/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_;
}
diff --git a/src/cell.hpp b/src/cell.hpp
index bc3d49e6289099b59e24e813c1a91e6b386c16c5..214a5b3c9ba090a2c38e007b50f1594512d53268 100644
--- a/src/cell.hpp
+++ b/src/cell.hpp
@@ -1,7 +1,9 @@
#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
diff --git a/src/segment.hpp b/src/segment.hpp
index 7710ed9714aa112dea731cf27bfd648027bcca16..a2429032f28d97aa61f87eb31f79b9783eca4246 100644
--- a/src/segment.hpp
+++ b/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 :
diff --git a/src/stimulus.hpp b/src/stimulus.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..94a4dddfddf1cab80272de74bd3a94f9130989fd
--- /dev/null
+++ b/src/stimulus.hpp
@@ -0,0 +1,52 @@
+#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
diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt
index 3bfe22c59526798ea1108ad1862a7e418a734182..53db3ff016a2765cecef53d807b20d09f5c8eaa1 100644
--- a/tests/CMakeLists.txt
+++ b/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
diff --git a/tests/test_cell.cpp b/tests/test_cell.cpp
index d9cb96b8e816d1942438d98ff13a5ee3e1bba940..68748512f594702d84a68d89c3eb72ec8010af25 100644
--- a/tests/test_cell.cpp
+++ b/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);
diff --git a/tests/test_run.cpp b/tests/test_run.cpp
index 9ac528f7a5768c08b8043d019f9207b7339506e7..f80cb41f612d2fd53db9162ec6b8400634cf91bc 100644
--- a/tests/test_run.cpp
+++ b/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
}
diff --git a/tests/test_stimulus.cpp b/tests/test_stimulus.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..43be56f158e07e409c94b21cca7c293026b8d830
--- /dev/null
+++ b/tests/test_stimulus.cpp
@@ -0,0 +1,33 @@
+#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);
+}
+