diff --git a/src/cell.cpp b/src/cell.cpp
index 3c867c83fcd622001d0f24c3ac8ab22fd20f72ab..72bad02a72af98aa9ee63e0c5cd1dbd3239a040a 100644
--- a/src/cell.cpp
+++ b/src/cell.cpp
@@ -4,6 +4,20 @@
namespace nest {
namespace mc {
+int find_compartment_index(
+ segment_location const& location,
+ compartment_model const& graph
+) {
+ EXPECTS(location.segment<graph.segment_index.size());
+ const auto& si = graph.segment_index;
+ const auto seg = location.segment;
+
+ auto first = si[seg];
+ auto n = si[seg+1] - first;
+ auto index = std::floor(n*location.position);
+ return index<n ? first+index : first+n-1;
+}
+
cell::cell()
{
// insert a placeholder segment for the soma
@@ -162,6 +176,20 @@ compartment_model cell::model() const
return m;
}
+
+void cell::add_stimulus( segment_location loc, i_clamp stim)
+{
+ if(!(loc.segment<num_segments())) {
+ throw std::out_of_range(
+ util::pprintf(
+ "can't insert stimulus in segment % of a cell with % segments",
+ loc.segment, num_segments()
+ )
+ );
+ }
+ stimulii_.push_back({loc, std::move(stim)});
+}
+
std::vector<int> const& cell::segment_parents() const
{
return parents_;
diff --git a/src/cell.hpp b/src/cell.hpp
index 1c1a43c45e9f929e7aa2b2c141bb5584d3172f6c..6e5d95f90578c7afda672f40c0653b73613c18d7 100644
--- a/src/cell.hpp
+++ b/src/cell.hpp
@@ -5,8 +5,9 @@
#include <thread>
#include <vector>
-#include "segment.hpp"
-#include "cell_tree.hpp"
+#include <segment.hpp>
+#include <cell_tree.hpp>
+#include <stimulus.hpp>
namespace nest {
namespace mc {
@@ -19,6 +20,21 @@ struct compartment_model {
std::vector<int> segment_index;
};
+struct segment_location {
+ segment_location(int s, double l)
+ : segment(s), position(l)
+ {
+ EXPECTS(position>=0. && position<=1.);
+ }
+ int segment;
+ double position;
+};
+
+int find_compartment_index(
+ segment_location const& location,
+ compartment_model const& graph
+);
+
/// high-level abstract representation of a cell and its segments
class cell {
public:
@@ -83,12 +99,28 @@ class cell {
compartment_model model() const;
+ void add_stimulus(segment_location loc, i_clamp stim);
+
+ std::vector<std::pair<segment_location, i_clamp>>&
+ stimulii() {
+ return stimulii_;
+ }
+
+ const std::vector<std::pair<segment_location, i_clamp>>&
+ stimulii() const {
+ return stimulii_;
+ }
+
private:
// storage for connections
std::vector<index_type> parents_;
+
// the segments
std::vector<segment_ptr> segments_;
+
+ // the stimulii
+ std::vector<std::pair<segment_location, i_clamp>> stimulii_;
};
// create a cable by forwarding cable construction parameters provided by the user
diff --git a/src/fvm.hpp b/src/fvm.hpp
index 33439dd95ae23431eb8d605c37eb2e441cc3ffe8..c69f76289a07b6061c2586eb37a84115507575ce 100644
--- a/src/fvm.hpp
+++ b/src/fvm.hpp
@@ -14,6 +14,7 @@
#include <mechanism_interface.hpp>
#include <util.hpp>
#include <segment.hpp>
+#include <stimulus.hpp>
#include <vector/include/Vector.hpp>
@@ -165,6 +166,8 @@ class fvm_cell {
/// the ion species
std::map<mechanisms::ionKind, ion_type> ions_;
+
+ std::vector<std::pair<uint32_t, i_clamp>> stimulii_;
};
////////////////////////////////////////////////////////////////////////////////
@@ -232,7 +235,6 @@ fvm_cell<T, I>::fvm_cell(nest::mc::cell const& cell)
auto radius_center = math::mean(c.radius);
auto area_face = math::area_circle( radius_center );
face_alpha_[i] = area_face / (c_m * r_L * c.length);
- cv_capacitance_[i] = c_m;
auto halflen = c.length/2;
@@ -371,6 +373,13 @@ fvm_cell<T, I>::fvm_cell(nest::mc::cell const& cell)
ion_ca().reversal_potential()(all) = 12.5 * std::log(2.0/5e-5);// mV
ion_ca().internal_concentration()(all) = 5e-5; // mM
ion_ca().external_concentration()(all) = 2.0; // mM
+
+ // add the stimulii
+ for(const auto& stim : cell.stimulii()) {
+ auto idx = find_compartment_index(stim.first, graph);
+ std::cout << "adding stimulus at compartment " << idx << "\n";
+ stimulii_.push_back( {idx, stim.second} );
+ }
}
template <typename T, typename I>
@@ -402,7 +411,7 @@ void fvm_cell<T, I>::setup_matrix(T dt)
for(auto i=1u; i<d.size(); ++i) {
// TODO get this right
// probably requires scaling a by cv_areas_[i] and cv_areas_[p[i]]
- auto a = 1e7*dt * face_alpha_[i];
+ auto a = 1e5*dt * face_alpha_[i];
d[i] += a;
l[i] = -a;
@@ -438,19 +447,20 @@ void fvm_cell<T, I>::advance(T dt)
current_(all) = 0.;
- // update currents
+ // update currents from ion channels
for(auto& m : mechanisms_) {
m->set_params(t_, dt);
m->nrn_current();
}
- // the factor scales the injected current to 10^2.nA
- auto ie_factor = 100.;
- auto ie = 0.1;
- auto loc = size()-1;
- //auto loc = 0;
- if(t_>=5. && t_<8.)
- current_[loc] -= ie_factor*ie/cv_areas_[loc];
+ // add current contributions from stimulii
+ for(auto& stim : stimulii_) {
+ auto ie = stim.second.amplitude(t_);
+ auto loc = stim.first;
+
+ // the factor of 100 scales the injected current to 10^2.nA
+ current_[loc] -= 100.*ie/cv_areas_[loc];
+ }
//std::cout << "t " << t_ << " current " << current_;
diff --git a/src/stimulus.hpp b/src/stimulus.hpp
index 94a4dddfddf1cab80272de74bd3a94f9130989fd..f7c7538cf398a7154588126620d97b7c054953a4 100644
--- a/src/stimulus.hpp
+++ b/src/stimulus.hpp
@@ -1,6 +1,7 @@
#pragma once
-namespace nestmc {
+namespace nest {
+namespace mc {
class i_clamp {
public:
@@ -49,4 +50,5 @@ class i_clamp {
value_type amplitude_ = 0;
};
-} // namespace nestmc
+} // namespace mc
+} // namespace nest
diff --git a/tests/test_run.cpp b/tests/test_run.cpp
index 12012f0f8fb885f50152edefdf87e05fd5010ef8..82ac6fd452039bfb7f3d703c4145f5bc8601a931 100644
--- a/tests/test_run.cpp
+++ b/tests/test_run.cpp
@@ -1,8 +1,35 @@
+#include <fstream>
+
#include "gtest.h"
#include "../src/cell.hpp"
#include "../src/fvm.hpp"
+void write_vis_file(const std::string& fname, std::vector<std::vector<double>> values)
+{
+ auto m = values.size();
+ if(!m) return;
+
+ std::ofstream fid(fname);
+ if(!fid.is_open()) return;
+
+ auto n = values[0].size();
+ for(const auto& v : values) {
+ std::cout << n << " ----- " << v.size() << "\n";
+ if(n!=v.size()) {
+ std::cerr << "all output arrays must have the same length\n";
+ return;
+ }
+ }
+
+ for(auto i=0u; i<n; ++i) {
+ for(auto j=0u; j<m; ++j) {
+ fid << " " << values[j][i];
+ }
+ fid << "\n";
+ }
+}
+
// based on hh/Neuron/steps_A.py
TEST(run, single_compartment)
{
@@ -18,25 +45,20 @@ TEST(run, single_compartment)
soma->mechanism("membrane").set("r_L", 123); // no effect for single compartment cell
soma->add_mechanism(hh_parameters());
- std::cout << soma->mechanism("membrane");
+ // add stimulus to the soma
+ cell.add_stimulus({0,0}, {5., 3., 0.1});
// make the lowered finite volume cell
fvm::fvm_cell<double, int> model(cell);
- //std::cout << "CV areas " << model.cv_areas() << "\n";
-
- //std::cout << "-----------------------------\n";
-
// set initial conditions
using memory::all;
model.voltage()(all) = -65.;
model.initialize(); // have to do this _after_ initial conditions are set
- //std::cout << "-----------------------------\n";
-
// run the simulation
- auto dt = 0.02; // ms
- auto tfinal = 30.; // ms
+ auto dt = 0.025; // ms
+ auto tfinal = 25.; // ms
int nt = tfinal/dt;
std::vector<double> result;
result.push_back(model.voltage()[0]);
@@ -44,7 +66,6 @@ TEST(run, single_compartment)
model.advance(dt);
result.push_back(model.voltage()[0]);
}
- std::cout << "took " << nt << " time steps" << std::endl;
{
std::ofstream fid("v.dat");
@@ -60,6 +81,8 @@ TEST(run, ball_and_stick)
{
using namespace nest::mc;
+ std::vector<std::vector<double>> results(4);
+
nest::mc::cell cell;
// setup global state for the mechanisms
@@ -71,9 +94,13 @@ TEST(run, ball_and_stick)
// add dendrite of length 200 um and diameter 1 um with passive channel
auto dendrite = cell.add_cable(0, segmentKind::dendrite, 0.5, 0.5, 200);
- dendrite->set_compartments(5); // 5 compartments
+ dendrite->set_compartments(81); // 5 compartments
dendrite->add_mechanism(pas_parameters());
+ // add stimulus
+ //cell.add_stimulus({0,0}, {5., 3., 0.1}); // soma
+ cell.add_stimulus({1,1}, {5., 30., 0.1});
+
// make the lowered finite volume cell
fvm::fvm_cell<double, int> model(cell);
@@ -86,26 +113,21 @@ TEST(run, ball_and_stick)
// run the simulation
auto dt = 0.02; // ms
- auto tfinal = 20.; // ms
+ auto tfinal = 25.; // ms
int nt = tfinal/dt;
- std::vector<double> result;
- result.push_back(model.voltage()[0]);
+ results[0].push_back(0.);
+ results[1].push_back(model.voltage()[0]);
+ results[2].push_back(model.voltage()[10]);
+ results[3].push_back(model.voltage()[20]);
for(auto i=0; i<nt; ++i) {
model.advance(dt);
- result.push_back(model.voltage()[0]);
+ results[0].push_back((i+1)*dt);
+ results[1].push_back(model.voltage()[0]);
+ results[2].push_back(model.voltage()[10]);
+ results[3].push_back(model.voltage()[20]);
}
- std::cout << "took " << nt << " time steps" << std::endl;
- std::cout << model.voltage() << "\n";
-
- {
- std::ofstream fid("v.dat");
- auto t = 0.;
- for(auto v:result) {
- fid << t << " " << v << "\n";
- t += dt;
- }
- }
+ write_vis_file("v.dat", results);
}
TEST(run, cable)
diff --git a/tests/test_stimulus.cpp b/tests/test_stimulus.cpp
index 43be56f158e07e409c94b21cca7c293026b8d830..8fd9307b61a780ad93f430338f7d76c941e20bc0 100644
--- a/tests/test_stimulus.cpp
+++ b/tests/test_stimulus.cpp
@@ -4,7 +4,7 @@
TEST(stimulus, i_clamp)
{
- using namespace nestmc;
+ using namespace nest::mc;
// stimulus with delay 2, duration 0.5, amplitude 6.0
i_clamp stim(2.0, 0.5, 6.0);
diff --git a/tests/vis.py b/tests/vis.py
index 9e6290e137c56b552fc25e22b3d28217f5526de5..2a1834819da3f6d1c6e0865ec194c15ce637bb35 100644
--- a/tests/vis.py
+++ b/tests/vis.py
@@ -3,11 +3,20 @@ import numpy as np
data = np.loadtxt('v.dat')
t = data[:,0]
-v = data[:,1]
+nfields = data.shape[1]
+
+for i in range(1,nfields) :
+ v = data[:,i]
+ pyplot.plot(t, v)
-pyplot.plot(t, v)
-pyplot.ylim(-80, 50)
pyplot.xlabel('time (ms)')
pyplot.ylabel('mV')
+pyplot.grid()
pyplot.show()
+
+#mi = min(v)
+#ma = max(v)
+#rng = (ma-mi)/2 * 1.1
+#mn = 0.5*(mi+ma)
+#pyplot.ylim(mn-rng, mn+rng)