From 4018fbf3c8c2bcae11b02a81b994e884dc652c90 Mon Sep 17 00:00:00 2001
From: bcumming <bcumming@cscs.ch>
Date: Wed, 4 May 2016 16:26:37 +0200
Subject: [PATCH] very close to working prototype
---
.gitignore | 1 +
mechanisms/mod/hh.mod | 4 +-
mechanisms/mod/pas.mod | 2 +-
src/cell.cpp | 8 ++-
src/cell.hpp | 10 ++--
src/fvm.hpp | 97 +++++++++++++++++++++++++++++++------
src/matrix.hpp | 55 +++++++++++++++------
src/mechanism.hpp | 27 ++++++++---
src/mechanism_interface.hpp | 5 +-
src/parameter_list.hpp | 19 +++++---
tests/test_mechanisms.cpp | 9 ++--
tests/test_run.cpp | 88 +++++++++++++++++++++++++++++++++
tests/vis.py | 13 +++++
13 files changed, 279 insertions(+), 59 deletions(-)
create mode 100644 tests/vis.py
diff --git a/.gitignore b/.gitignore
index cd1f0915..f60f1b41 100644
--- a/.gitignore
+++ b/.gitignore
@@ -32,6 +32,7 @@
*.dot
*.pdf
*.jpg
+*.dat
# latex output
*.aux
diff --git a/mechanisms/mod/hh.mod b/mechanisms/mod/hh.mod
index 975d601d..c96280fd 100644
--- a/mechanisms/mod/hh.mod
+++ b/mechanisms/mod/hh.mod
@@ -29,11 +29,13 @@ NEURON {
}
PARAMETER {
+ : neuron uses S/cm2, i have to scale by a factor of 10
gnabar = .12 (S/cm2) : <0,1e9>
gkbar = .036 (S/cm2) : <0,1e9>
gl = .0003 (S/cm2) : <0,1e9>
el = -54.3 (mV)
- celsius = 37 (degC)
+ celsius = 6.3 (degC)
+ :celsius = 37 (degC)
}
STATE {
diff --git a/mechanisms/mod/pas.mod b/mechanisms/mod/pas.mod
index 413bd250..b7dc2ad6 100644
--- a/mechanisms/mod/pas.mod
+++ b/mechanisms/mod/pas.mod
@@ -16,7 +16,7 @@ INITIAL {}
PARAMETER {
g = .001 (S/cm2) :<0,1e9>
- e = -70 (mV)
+ e = -65 (mV) : we use -65 for the ball and stick model, instead of Neuron default of -70
}
ASSIGNED {
diff --git a/src/cell.cpp b/src/cell.cpp
index ab01007a..3c867c83 100644
--- a/src/cell.cpp
+++ b/src/cell.cpp
@@ -20,7 +20,7 @@ int cell::num_segments() const
// 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)
+soma_segment* cell::add_soma(value_type radius, point_type center)
{
if(has_soma()) {
throw std::domain_error(
@@ -35,9 +35,11 @@ void cell::add_soma(value_type radius, point_type center)
else {
segments_[0] = make_segment<soma_segment>(radius);
}
+
+ return segments_[0]->as_soma();
}
-void cell::add_cable(cell::index_type parent, segment_ptr&& cable)
+cable_segment* cell::add_cable(cell::index_type parent, segment_ptr&& cable)
{
// check for a valid parent id
if(cable->is_soma()) {
@@ -54,6 +56,8 @@ void cell::add_cable(cell::index_type parent, segment_ptr&& cable)
}
segments_.push_back(std::move(cable));
parents_.push_back(parent);
+
+ return segments_.back()->as_cable();
}
segment* cell::segment(int index)
diff --git a/src/cell.hpp b/src/cell.hpp
index 68bc4c5a..1c1a43c4 100644
--- a/src/cell.hpp
+++ b/src/cell.hpp
@@ -33,18 +33,18 @@ class cell {
/// add a soma to the cell
/// radius must be specified
- void add_soma(value_type radius, point_type center=point_type());
+ soma_segment* add_soma(value_type radius, point_type center=point_type());
/// add a cable
/// parent is the index of the parent segment for the cable section
/// cable is the segment that will be moved into the cell
- void add_cable(index_type parent, segment_ptr&& cable);
+ cable_segment* add_cable(index_type parent, segment_ptr&& cable);
/// add a cable by constructing it in place
/// parent is the index of the parent segment for the cable section
/// args are the arguments to be used to consruct the new cable
template <typename... Args>
- void add_cable(index_type parent, Args ...args);
+ cable_segment* add_cable(index_type parent, Args ...args);
/// the number of segments in the cell
int num_segments() const;
@@ -93,7 +93,7 @@ class cell {
// create a cable by forwarding cable construction parameters provided by the user
template <typename... Args>
-void cell::add_cable(cell::index_type parent, Args ...args)
+cable_segment* cell::add_cable(cell::index_type parent, Args ...args)
{
// check for a valid parent id
if(parent>=num_segments()) {
@@ -103,6 +103,8 @@ void cell::add_cable(cell::index_type parent, Args ...args)
}
segments_.push_back(make_segment<cable_segment>(std::forward<Args>(args)...));
parents_.push_back(parent);
+
+ return segments_.back()->as_cable();
}
} // namespace mc
diff --git a/src/fvm.hpp b/src/fvm.hpp
index 674c7d0e..d680edc8 100644
--- a/src/fvm.hpp
+++ b/src/fvm.hpp
@@ -79,6 +79,11 @@ class fvm_cell {
return cv_capacitance_;
}
+ /// return the voltage in each CV
+ vector_view voltage() {
+ return voltage_;
+ }
+
std::size_t size() const {
return matrix_.size();
}
@@ -121,8 +126,17 @@ class fvm_cell {
return ions_[mechanisms::ionKind::k];
}
+ /// make a time step
+ void advance(value_type dt);
+
+ /// set initial states
+ void initialize();
+
private:
+ /// current time
+ value_type t_ = value_type{0};
+
/// the linear system for implicit time stepping of cell state
matrix_type matrix_;
@@ -241,6 +255,8 @@ fvm_cell<T, I>::fvm_cell(nest::mc::cell const& cell)
cv_capacitance_[i] /= cv_areas_[i];
}
+ std::cout << "capacitance " << cv_capacitance_ << "\n";
+
/////////////////////////////////////////////
// create mechanisms
/////////////////////////////////////////////
@@ -290,12 +306,13 @@ fvm_cell<T, I>::fvm_cell(nest::mc::cell const& cell)
}
// instantiate the mechanism
+ index_view node_index(compartment_index.data(), compartment_index.size());
mechanisms_.push_back(
- helper->new_mechanism(
- &matrix_,
- index_view(compartment_index.data(), compartment_index.size())
- )
+ helper->new_mechanism(voltage_, current_, node_index)
);
+ std::cout << "created mech " << mech.first
+ << " with size " << mechanisms_.back()->size()
+ << " and indexes " << node_index << "\n";
}
/////////////////////////////////////////////
@@ -340,20 +357,23 @@ fvm_cell<T, I>::fvm_cell(nest::mc::cell const& cell)
// FIXME: Hard code parameters for now.
// Take defaults for reversal potential of sodium and potassium from
// the default values in Neuron.
- // We don't use the other parameters for the HH model, so I leave the NaN.
- // To set them I would have to go spelunking in the Neuron source.
+ // Neuron's defaults are defined in the file
+ // nrn/src/nrnoc/membdef.h
using memory::all;
- ion_ca().reversal_potential()(all) = std::numeric_limits<value_type>::quiet_NaN();
- ion_ca().internal_concentration()(all) = std::numeric_limits<value_type>::quiet_NaN();
- ion_ca().external_concentration()(all) = std::numeric_limits<value_type>::quiet_NaN();
- ion_na().reversal_potential()(all) = -50.0;
- ion_na().internal_concentration()(all) = std::numeric_limits<value_type>::quiet_NaN();
- ion_na().external_concentration()(all) = std::numeric_limits<value_type>::quiet_NaN();
+ constexpr value_type DEF_vrest = -65.0; // same name as #define in Neuron
+
+ ion_na().reversal_potential()(all) = 115+DEF_vrest; // mV
+ ion_na().internal_concentration()(all) = 10.0; // mM
+ ion_na().external_concentration()(all) = 140.0; // mM
+
+ ion_k().reversal_potential()(all) = -12.0+DEF_vrest;// mV
+ ion_k().internal_concentration()(all) = 54.4; // mM
+ ion_k().external_concentration()(all) = 2.5; // mM
- ion_k().reversal_potential()(all) = -77.0;
- ion_k().internal_concentration()(all) = std::numeric_limits<value_type>::quiet_NaN();
- ion_k().external_concentration()(all) = std::numeric_limits<value_type>::quiet_NaN();
+ 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
}
template <typename T, typename I>
@@ -399,6 +419,53 @@ void fvm_cell<T, I>::setup_matrix(T dt)
}
}
+template <typename T, typename I>
+void fvm_cell<T, I>::initialize()
+{
+ // initialize mechanism states
+ for(auto& m : mechanisms_) {
+ m->nrn_init();
+ }
+}
+
+template <typename T, typename I>
+void fvm_cell<T, I>::advance(T dt)
+{
+ using memory::all;
+
+ current_(all) = 0.;
+
+ // update currents
+ for(auto& m : mechanisms_) {
+ m->set_params(t_, dt);
+ m->nrn_current();
+ }
+
+ //if(t_>=10.) {
+ current_[0] -= 0.1;
+ //}
+ //std::cout << "t " << t_ << " current " << current_;
+
+ // set matrix diagonals and rhs
+ setup_matrix(dt);
+
+ //std::cout << " rhs " << matrix_.rhs() << " d " << matrix_.d();
+
+ // solve the linear system
+ matrix_.solve();
+
+ voltage_(all) = matrix_.rhs();
+
+ //std::cout << " v " << voltage_ << "\n";
+
+ // update states
+ for(auto& m : mechanisms_) {
+ m->nrn_state();
+ }
+
+ t_ += dt;
+}
+
} // namespace fvm
} // namespace mc
} // namespace nest
diff --git a/src/matrix.hpp b/src/matrix.hpp
index e6d406ee..db9d2df7 100644
--- a/src/matrix.hpp
+++ b/src/matrix.hpp
@@ -19,10 +19,10 @@ class matrix {
using size_type = I;
// define storage types
- using vector_type = memory::HostVector<value_type>;
- using view_type = typename vector_type::view_type;
- using index_type = memory::HostVector<size_type>;
- using index_view = typename index_type::view_type;
+ using vector_type = memory::HostVector<value_type>;
+ using vector_view_type = typename vector_type::view_type;
+ using index_type = memory::HostVector<size_type>;
+ using index_view_type = typename index_type::view_type;
matrix() = default;
@@ -81,42 +81,44 @@ class matrix {
}
/// FIXME : make modparser use the correct accessors (l,d,u,rhs) instead of these
- view_type vec_rhs() { return rhs_; }
- view_type vec_d() { return d_; }
- view_type vec_v() { return v_; }
+ vector_view_type vec_rhs() { return rhs(); }
+ vector_view_type vec_d() { return d(); }
+ vector_view_type vec_v() { return v(); }
/// the vector holding the lower part of the matrix
- view_type l()
+ vector_view_type l()
{
return l_;
}
/// the vector holding the diagonal of the matrix
- view_type d()
+ vector_view_type d()
{
return d_;
}
/// the vector holding the upper part of the matrix
- view_type u()
+ vector_view_type u()
{
return u_;
}
/// the vector holding the right hand side of the linear equation system
- view_type rhs()
+ vector_view_type rhs()
{
return rhs_;
}
/// the vector holding the solution (voltage)
- view_type v()
+ vector_view_type v()
{
+ EXPECTS(has_voltage_);
+
return v_;
}
/// the vector holding the parent index
- index_view p()
+ index_view_type p()
{
return parent_index_;
}
@@ -126,7 +128,16 @@ class matrix {
/// and can be accessed via rhs()
void solve()
{
- index_view const& p = parent_index_;
+ /*
+ std::cout << "solving matrix :\n";
+ std::cout << " l " << l_ << "\n";
+ std::cout << " d " << d_ << "\n";
+ std::cout << " u " << u_ << "\n";
+ std::cout << " rhs " << rhs_ << "\n";
+ */
+
+ // FIXME make a const view
+ index_view_type const& p = parent_index_;
auto const ncells = num_cells();
// loop over submatrices
@@ -148,6 +159,17 @@ class matrix {
rhs_[i] /= d_[i];
}
}
+ //std::cout << " v " << rhs_ << "\n";
+ }
+
+ void add_voltage(vector_view_type v_ext)
+ {
+ EXPECTS(v_ext.size()==size());
+
+ std::cout << "============ adding voltage" << std::endl;
+
+ v_ = v_ext;
+ has_voltage_ = true;
}
private:
@@ -163,7 +185,6 @@ class matrix {
d_ = vector_type(n, default_value);
u_ = vector_type(n, default_value);
rhs_ = vector_type(n, default_value);
- v_ = vector_type(n, default_value);
}
/// the parent indice that describe matrix structure
@@ -179,7 +200,9 @@ class matrix {
/// after calling solve, the solution is stored in rhs_
vector_type rhs_;
- vector_type v_; // should this be in the matrix? not so sure...
+ vector_view_type v_;
+
+ bool has_voltage_=false;
};
} // namespace nest
diff --git a/src/mechanism.hpp b/src/mechanism.hpp
index f43abc16..ea69c292 100644
--- a/src/mechanism.hpp
+++ b/src/mechanism.hpp
@@ -6,7 +6,6 @@
#include <string>
#include "indexed_view.hpp"
-#include "matrix.hpp"
#include "parameter_list.hpp"
#include "util.hpp"
#include "ion.hpp"
@@ -32,11 +31,11 @@ public:
using index_view = typename index_type::view_type;
using indexed_view_type = indexed_view<value_type, size_type>;
- using matrix_type = matrix<value_type, size_type>;
using ion_type = ion<value_type, size_type>;
- mechanism(matrix_type *matrix, index_view node_index)
- : matrix_(matrix)
+ mechanism(view_type vec_v, view_type vec_i, index_view node_index)
+ : vec_v_(vec_v)
+ , vec_i_(vec_i)
, node_index_(node_index)
{}
@@ -50,6 +49,16 @@ public:
return node_index_;
}
+ value_type voltage(size_type i) const
+ {
+ return vec_v_[node_index_[i]];
+ }
+
+ value_type current(size_type i) const
+ {
+ return vec_i_[node_index_[i]];
+ }
+
virtual void set_params(value_type t_, value_type dt_) = 0;
virtual std::string name() const = 0;
virtual std::size_t memory() const = 0;
@@ -61,7 +70,8 @@ public:
virtual mechanismKind kind() const = 0;
- matrix_type* matrix_;
+ view_type vec_v_;
+ view_type vec_i_;
index_type node_index_;
};
@@ -71,10 +81,11 @@ using mechanism_ptr = std::unique_ptr<mechanism<T,I>>;
template <typename M>
mechanism_ptr<typename M::value_type, typename M::size_type>
make_mechanism(
- typename M::matrix_type* matrix,
- typename M::index_view node_indices
+ typename M::view_type vec_v,
+ typename M::view_type vec_i,
+ typename M::index_view node_indices
) {
- return util::make_unique<M>(matrix, node_indices);
+ return util::make_unique<M>(vec_v, vec_i, node_indices);
}
} // namespace mechanisms
diff --git a/src/mechanism_interface.hpp b/src/mechanism_interface.hpp
index 33b6e7f2..c59caaad 100644
--- a/src/mechanism_interface.hpp
+++ b/src/mechanism_interface.hpp
@@ -3,7 +3,6 @@
#include <map>
#include <string>
-#include "matrix.hpp"
#include "mechanism.hpp"
#include "parameter_list.hpp"
@@ -23,11 +22,11 @@ struct mechanism_helper {
using index_type = memory::HostVector<I>;
using index_view = typename index_type::view_type;
using mechanism_ptr_type = mechanism_ptr<T, I>;
- using matrix_type = typename mechanism<T,I>::matrix_type;
+ using view_type = typename mechanism<T,I>::view_type;
virtual std::string name() const = 0;
- virtual mechanism_ptr<T,I> new_mechanism(matrix_type*, index_view) const = 0;
+ virtual mechanism_ptr<T,I> new_mechanism(view_type, view_type, index_view) const = 0;
virtual void set_parameters(mechanism_ptr_type&, parameter_list const&) const = 0;
};
diff --git a/src/parameter_list.hpp b/src/parameter_list.hpp
index 1080ac49..0c069887 100644
--- a/src/parameter_list.hpp
+++ b/src/parameter_list.hpp
@@ -146,8 +146,15 @@ namespace mc {
membrane_parameters()
: base("membrane")
{
- base::add_parameter({"r_L", 0.01, {0., 1e9}}); // typically 10 nF/mm^2 == 0.01 F/m2
- base::add_parameter({"c_m", 180.00, {0., 1e9}}); // Ohm.cm
+ // from nrn/src/nrnoc/membdef.h
+ //#define DEF_cm 1. // uF/cm^2
+ //#define DEF_Ra 35.4 // ohm-cm
+ //#define DEF_celsius 6.3 // deg-C
+ //#define DEF_vrest -65. // mV
+ // r_L is called Ra in Neuron
+ //base::add_parameter({"c_m", 10e-6, {0., 1e9}}); // typically 10 nF/mm^2 == 0.01 F/m^2 == 10^-6 F/cm^2
+ base::add_parameter({"c_m", 0.01, {0., 1e9}}); // typically 10 nF/mm^2 == 0.01 F/m^2 == 10^-6 F/cm^2
+ base::add_parameter({"r_L", 180.00, {0., 1e9}}); // equivalent to Ra in Neuron : Ohm.cm
}
};
@@ -169,10 +176,10 @@ namespace mc {
hh_parameters()
: base("hh")
{
- base::add_parameter({"gnabar", 0.12, {0, 1e9}});
- base::add_parameter({"gkbar", 0.036, {0, 1e9}});
- base::add_parameter({"gl", 0.0003,{0, 1e9}});
- base::add_parameter({"el", -54.3});
+ base::add_parameter({"gnabar", 0.12, {0, 1e9}});
+ base::add_parameter({"gkbar", 0.036, {0, 1e9}});
+ base::add_parameter({"gl", 0.0003,{0, 1e9}});
+ base::add_parameter({"el", -54.3});
}
};
diff --git a/tests/test_mechanisms.cpp b/tests/test_mechanisms.cpp
index f4e5a4d2..d61c89d2 100644
--- a/tests/test_mechanisms.cpp
+++ b/tests/test_mechanisms.cpp
@@ -22,13 +22,16 @@ TEST(mechanisms, helpers) {
//0 1 2 3 4 5 6 7 8 9
std::vector<int> parent_index = {0,0,1,2,3,4,0,6,7,8};
memory::HostVector<int> node_indices = std::vector<int>{0,6,7,8,9};
+ auto n = node_indices.size();
- nest::mc::matrix<double, int> matrix(parent_index);
+ //nest::mc::matrix<double, int> matrix(parent_index);
+ memory::HostVector<double> vec_i(n, 0.);
+ memory::HostVector<double> vec_v(n, 0.);
auto& helper = nest::mc::mechanisms::get_mechanism_helper("hh");
- auto mech = helper->new_mechanism(&matrix, node_indices);
+ auto mech = helper->new_mechanism(vec_v, vec_i, node_indices);
EXPECT_EQ(mech->name(), "hh");
EXPECT_EQ(mech->size(), 5u);
- EXPECT_EQ(mech->matrix_, &matrix);
+ //EXPECT_EQ(mech->matrix_, &matrix);
}
diff --git a/tests/test_run.cpp b/tests/test_run.cpp
index 55fb0d0f..0dc0b495 100644
--- a/tests/test_run.cpp
+++ b/tests/test_run.cpp
@@ -3,6 +3,94 @@
#include "../src/cell.hpp"
#include "../src/fvm.hpp"
+// based on hh/Neuron/steps_B.py
+TEST(run, single_compartment)
+{
+ using namespace nest::mc;
+
+ nest::mc::cell cell;
+
+ // setup global state for the mechanisms
+ nest::mc::mechanisms::setup_mechanism_helpers();
+
+ // Soma with diameter 18.8um -> 1.88e-3 cm and HH channel
+ auto soma = cell.add_soma(1.88e-3/2.0);
+ soma->add_mechanism(hh_parameters());
+
+ std::cout << soma->mechanism("membrane");
+
+ // 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 / 1000.; // convert ms to s
+ //auto tfinal = 100./1000.;
+ auto dt = 0.02; // ms
+ auto tfinal = 5; // ms
+ int nt = tfinal/dt;
+ std::vector<double> result;
+ result.push_back(model.voltage()[0]);
+ for(auto i=0; i<nt; ++i) {
+ model.advance(dt);
+ result.push_back(model.voltage()[0]);
+ }
+ std::cout << "took " << nt << " time steps" << std::endl;
+
+ {
+ std::ofstream fid("v.dat");
+ auto t = 0.;
+ for(auto v:result) {
+ fid << t << " " << v << "\n";
+ t += dt;
+ }
+ }
+}
+
+TEST(run, ball_and_stick)
+{
+ using namespace nest::mc;
+
+ nest::mc::cell cell;
+
+ // setup global state for the mechanisms
+ nest::mc::mechanisms::setup_mechanism_helpers();
+
+ // Soma with diameter 12.6157 and HH channel
+ auto soma = cell.add_soma(12.6157/2.0);
+ soma->add_mechanism(hh_parameters());
+
+ // add dendrite with passive channel and 10 compartments
+ auto dendrite = cell.add_cable(0, segmentKind::dendrite, 0.5, 0.5, 200);
+ dendrite->set_compartments(5);
+ dendrite->add_mechanism(pas_parameters());
+
+ // make the lowered finite volume cell
+ fvm::fvm_cell<double, int> model(cell);
+
+ std::cout << "CV areas " << model.cv_areas() << "\n";
+
+ // set initial conditions
+ using memory::all;
+ model.voltage()(all) = -65.;
+
+ // run the simulation
+ auto dt = 0.02; // ms
+ model.advance(dt);
+
+ std::cout << model.voltage() << "\n";
+}
+
TEST(run, cable)
{
using namespace nest::mc;
diff --git a/tests/vis.py b/tests/vis.py
new file mode 100644
index 00000000..751c1fd2
--- /dev/null
+++ b/tests/vis.py
@@ -0,0 +1,13 @@
+from matplotlib import pyplot
+import numpy as np
+
+data = np.loadtxt('v.dat')
+t = data[:,0]
+v = data[:,1]
+
+
+pyplot.plot(t, v)
+pyplot.xlabel('time (ms)')
+pyplot.ylabel('mV')
+pyplot.show()
+
--
GitLab