diff --git a/arbor/CMakeLists.txt b/arbor/CMakeLists.txt
index 923a4c97af3f8f117e3220f21f7c6dc8b40d8275..978d3a30514d3d189d0d51e5f0f944a0a4034fcf 100644
--- a/arbor/CMakeLists.txt
+++ b/arbor/CMakeLists.txt
@@ -10,6 +10,8 @@ set(arbor_sources
communication/dry_run_context.cpp
benchmark_cell_group.cpp
builtin_mechanisms.cpp
+ cable_cell.cpp
+ cable_cell_param.cpp
cell_group_factory.cpp
common_types_io.cpp
execution_context.cpp
@@ -22,7 +24,6 @@ set(arbor_sources
io/locked_ostream.cpp
io/serialize_hex.cpp
lif_cell_group.cpp
- cable_cell.cpp
mc_cell_group.cpp
mechcat.cpp
memory/cuda_wrappers.cpp
diff --git a/arbor/backends/gpu/shared_state.cpp b/arbor/backends/gpu/shared_state.cpp
index 5c29d9fb31dafc5989c376651d0cd83d699bdfdc..e15692beb7962f6f06a534cbc6ad26bbcf03d1bb 100644
--- a/arbor/backends/gpu/shared_state.cpp
+++ b/arbor/backends/gpu/shared_state.cpp
@@ -3,12 +3,12 @@
#include <arbor/constants.hpp>
#include <arbor/fvm_types.hpp>
-#include <arbor/ion_info.hpp>
#include "backends/event.hpp"
#include "backends/gpu/gpu_store_types.hpp"
#include "backends/gpu/shared_state.hpp"
#include "backends/multi_event_stream_state.hpp"
+#include "memory/copy.hpp"
#include "memory/wrappers.hpp"
#include "util/rangeutil.hpp"
@@ -19,14 +19,6 @@ namespace gpu {
// CUDA implementation entry points:
-void init_concentration_impl(
- std::size_t n, fvm_value_type* Xi, fvm_value_type* Xo, const fvm_value_type* weight_Xi,
- const fvm_value_type* weight_Xo, fvm_value_type iconc, fvm_value_type econc);
-
-void nernst_impl(
- std::size_t n, fvm_value_type factor,
- const fvm_value_type* charge, const fvm_value_type* Xi, const fvm_value_type* Xo, fvm_value_type* eX);
-
void update_time_to_impl(
std::size_t n, fvm_value_type* time_to, const fvm_value_type* time,
fvm_value_type dt, fvm_value_type tmax);
@@ -46,6 +38,8 @@ void take_samples_impl(
const multi_event_stream_state<raw_probe_info>& s,
const fvm_value_type* time, fvm_value_type* sample_time, fvm_value_type* sample_value);
+void add_scalar(std::size_t n, fvm_value_type* data, fvm_value_type v);
+
// GPU-side minmax: consider CUDA kernel replacement.
std::pair<fvm_value_type, fvm_value_type> minmax_value_impl(fvm_size_type n, const fvm_value_type* v) {
auto v_copy = memory::on_host(memory::const_device_view<fvm_value_type>(v, n));
@@ -55,10 +49,11 @@ std::pair<fvm_value_type, fvm_value_type> minmax_value_impl(fvm_size_type n, con
// Ion state methods:
ion_state::ion_state(
- ion_info info,
+ int charge,
const std::vector<fvm_index_type>& cv,
- const std::vector<fvm_value_type>& iconc_norm_area,
- const std::vector<fvm_value_type>& econc_norm_area,
+ const std::vector<fvm_value_type>& init_iconc,
+ const std::vector<fvm_value_type>& init_econc,
+ const std::vector<fvm_value_type>& init_erev,
unsigned // alignment/padding ignored.
):
node_index_(make_const_view(cv)),
@@ -66,53 +61,39 @@ ion_state::ion_state(
eX_(cv.size(), NAN),
Xi_(cv.size(), NAN),
Xo_(cv.size(), NAN),
- weight_Xi_(make_const_view(iconc_norm_area)),
- weight_Xo_(make_const_view(econc_norm_area)),
- charge(1u, info.charge),
- default_int_concentration(info.default_int_concentration),
- default_ext_concentration(info.default_ext_concentration)
+ init_Xi_(make_const_view(init_iconc)),
+ init_Xo_(make_const_view(init_econc)),
+ init_eX_(make_const_view(init_erev)),
+ charge(1u, charge)
{
- arb_assert(node_index_.size()==weight_Xi_.size());
- arb_assert(node_index_.size()==weight_Xo_.size());
-}
-
-void ion_state::nernst(fvm_value_type temperature_K) {
- // Nernst equation: reversal potenial eX given by:
- //
- // eX = RT/zF * ln(Xo/Xi)
- //
- // where:
- // R: universal gas constant 8.3144598 J.K-1.mol-1
- // T: temperature in Kelvin
- // z: valency of species (K, Na: +1) (Ca: +2)
- // F: Faraday's constant 96485.33289 C.mol-1
- // Xo/Xi: ratio of out/in concentrations
-
- // 1e3 factor required to scale from V -> mV.
- constexpr fvm_value_type RF = 1e3*constant::gas_constant/constant::faraday;
-
- fvm_value_type factor = RF*temperature_K;
- nernst_impl(Xi_.size(), factor, charge.data(), Xo_.data(), Xi_.data(), eX_.data());
+ arb_assert(node_index_.size()==init_Xi_.size());
+ arb_assert(node_index_.size()==init_Xo_.size());
+ arb_assert(node_index_.size()==init_eX_.size());
}
void ion_state::init_concentration() {
- init_concentration_impl(
- Xi_.size(),
- Xi_.data(), Xo_.data(),
- weight_Xi_.data(), weight_Xo_.data(),
- default_int_concentration, default_ext_concentration);
+ memory::copy(init_Xi_, Xi_);
+ memory::copy(init_Xo_, Xo_);
}
void ion_state::zero_current() {
memory::fill(iX_, 0);
}
+void ion_state::reset() {
+ zero_current();
+ init_concentration();
+ memory::copy(init_eX_, eX_);
+}
+
// Shared state methods:
shared_state::shared_state(
fvm_size_type n_intdom,
const std::vector<fvm_index_type>& cv_to_intdom_vec,
const std::vector<fvm_gap_junction>& gj_vec,
+ const std::vector<fvm_value_type>& init_membrane_potential,
+ const std::vector<fvm_value_type>& temperature_K,
unsigned // alignment parameter ignored.
):
n_intdom(n_intdom),
@@ -127,32 +108,35 @@ shared_state::shared_state(
voltage(n_cv),
current_density(n_cv),
conductivity(n_cv),
- temperature_degC(1),
+ init_voltage(make_const_view(init_membrane_potential)),
+ temperature_degC(make_const_view(temperature_K)),
deliverable_events(n_intdom)
-{}
+{
+ add_scalar(temperature_degC.size(), temperature_degC.data(), -273.15);
+}
void shared_state::add_ion(
const std::string& ion_name,
- ion_info info,
+ int charge,
const std::vector<fvm_index_type>& cv,
- const std::vector<fvm_value_type>& iconc_norm_area,
- const std::vector<fvm_value_type>& econc_norm_area)
+ const std::vector<fvm_value_type>& init_iconc,
+ const std::vector<fvm_value_type>& init_econc,
+ const std::vector<fvm_value_type>& init_erev)
{
ion_data.emplace(std::piecewise_construct,
std::forward_as_tuple(ion_name),
- std::forward_as_tuple(info, cv, iconc_norm_area, econc_norm_area, 1u));
+ std::forward_as_tuple(charge, cv, init_iconc, init_econc, init_erev, 1u));
}
-void shared_state::reset(fvm_value_type initial_voltage, fvm_value_type temperature_K) {
- memory::fill(voltage, initial_voltage);
+void shared_state::reset() {
+ memory::copy(init_voltage, voltage);
memory::fill(current_density, 0);
memory::fill(conductivity, 0);
memory::fill(time, 0);
memory::fill(time_to, 0);
- memory::fill(temperature_degC, temperature_K - 273.15);
for (auto& i: ion_data) {
- i.second.reset(temperature_K);
+ i.second.reset();
}
}
@@ -170,12 +154,6 @@ void shared_state::ions_init_concentration() {
}
}
-void shared_state::ions_nernst_reversal_potential(fvm_value_type temperature_K) {
- for (auto& i: ion_data) {
- i.second.nernst(temperature_K);
- }
-}
-
void shared_state::update_time_to(fvm_value_type dt_step, fvm_value_type tmax) {
update_time_to_impl(n_intdom, time_to.data(), time.data(), dt_step, tmax);
}
@@ -203,21 +181,26 @@ void shared_state::take_samples(const sample_event_stream::state& s, array& samp
// Debug interface
std::ostream& operator<<(std::ostream& o, shared_state& s) {
o << " cv_to_intdom " << s.cv_to_intdom << "\n";
- o << " time " << s.time << "\n";
- o << " time_to " << s.time_to << "\n";
+ o << " time " << s.time << "\n";
+ o << " time_to " << s.time_to << "\n";
o << " dt_intdom " << s.dt_intdom << "\n";
- o << " dt_cv " << s.dt_cv << "\n";
- o << " voltage " << s.voltage << "\n";
- o << " current " << s.current_density << "\n";
+ o << " dt_cv " << s.dt_cv << "\n";
+ o << " voltage " << s.voltage << "\n";
+ o << " init_voltage " << s.init_voltage << "\n";
+ o << " temperature " << s.temperature_degC << "\n";
+ o << " current " << s.current_density << "\n";
o << " conductivity " << s.conductivity << "\n";
for (auto& ki: s.ion_data) {
auto& kn = ki.first;
auto& i = const_cast<ion_state&>(ki.second);
- o << " " << kn << ".current_density " << i.iX_ << "\n";
- o << " " << kn << ".reversal_potential " << i.eX_ << "\n";
- o << " " << kn << ".internal_concentration " << i.Xi_ << "\n";
- o << " " << kn << ".external_concentration " << i.Xo_ << "\n";
- o << " " << kn << ".node_index " << i.node_index_ << "\n";
+ o << " " << kn << "/current_density " << i.iX_ << "\n";
+ o << " " << kn << "/reversal_potential " << i.eX_ << "\n";
+ o << " " << kn << "/internal_concentration " << i.Xi_ << "\n";
+ o << " " << kn << "/external_concentration " << i.Xo_ << "\n";
+ o << " " << kn << "/intconc_initial " << i.init_Xi_ << "\n";
+ o << " " << kn << "/extconc_initial " << i.init_Xo_ << "\n";
+ o << " " << kn << "/revpot_initial " << i.init_eX_ << "\n";
+ o << " " << kn << "/node_index " << i.node_index_ << "\n";
}
return o;
}
diff --git a/arbor/backends/gpu/shared_state.cu b/arbor/backends/gpu/shared_state.cu
index 34912740283b8fa07764c257b0a5b32f9898bc44..78c8a0919147c1789e873fb696c1cb59265cada9 100644
--- a/arbor/backends/gpu/shared_state.cu
+++ b/arbor/backends/gpu/shared_state.cu
@@ -13,27 +13,6 @@ namespace gpu {
namespace kernel {
-template <typename T>
-__global__
-void nernst_impl(unsigned n, T factor, const T* charge, const T* Xo, const T* Xi, T* eX) {
- unsigned i = threadIdx.x+blockIdx.x*blockDim.x;
-
- if (i<n) {
- eX[i] = factor*std::log(Xo[i]/Xi[i])/charge[0];
- }
-}
-
-template <typename T>
-__global__
-void init_concentration_impl(unsigned n, T* Xi, T* Xo, const T* weight_Xi, const T* weight_Xo, T c_int, T c_ext) {
- unsigned i = threadIdx.x+blockIdx.x*blockDim.x;
-
- if (i<n) {
- Xi[i] = c_int*weight_Xi[i];
- Xo[i] = c_ext*weight_Xo[i];
- }
-}
-
template <typename T>
__global__ void update_time_to_impl(unsigned n, T* time_to, const T* time, T dt, T tmax) {
unsigned i = threadIdx.x+blockIdx.x*blockDim.x;
@@ -54,6 +33,15 @@ __global__ void add_gj_current_impl(unsigned n, const T* gj_info, const I* volta
}
}
+// Vector/scalar addition: x[i] += v ∀i
+template <typename T>
+__global__ void add_scalar(unsigned n, T* x, fvm_value_type v) {
+ unsigned i = threadIdx.x+blockIdx.x*blockDim.x;
+ if (i<n) {
+ x[i] += v;
+ }
+}
+
// Vector minus: x = y - z
template <typename T>
__global__ void vec_minus(unsigned n, T* x, const T* y, const T* z) {
@@ -91,26 +79,12 @@ __global__ void take_samples_impl(
using impl::block_count;
-void nernst_impl(
- std::size_t n, fvm_value_type factor,
- const fvm_value_type* charge, const fvm_value_type* Xo, const fvm_value_type* Xi, fvm_value_type* eX)
-{
- if (!n) return;
-
- constexpr int block_dim = 128;
- int nblock = block_count(n, block_dim);
- kernel::nernst_impl<<<nblock, block_dim>>>(n, factor, charge, Xo, Xi, eX);
-}
-
-void init_concentration_impl(
- std::size_t n, fvm_value_type* Xi, fvm_value_type* Xo, const fvm_value_type* weight_Xi,
- const fvm_value_type* weight_Xo, fvm_value_type c_int, fvm_value_type c_ext)
-{
+void add_scalar(std::size_t n, fvm_value_type* data, fvm_value_type v) {
if (!n) return;
constexpr int block_dim = 128;
- int nblock = block_count(n, block_dim);
- kernel::init_concentration_impl<<<nblock, block_dim>>>(n, Xi, Xo, weight_Xi, weight_Xo, c_int, c_ext);
+ const int nblock = block_count(n, block_dim);
+ kernel::add_scalar<<<nblock, block_dim>>>(n, data, v);
}
void update_time_to_impl(
diff --git a/arbor/backends/gpu/shared_state.hpp b/arbor/backends/gpu/shared_state.hpp
index 6aea214e9352a13edeb80622f2bbe45c44265b7e..a2d99456a8bdbc007899f4027384735c0802cb31 100644
--- a/arbor/backends/gpu/shared_state.hpp
+++ b/arbor/backends/gpu/shared_state.hpp
@@ -6,7 +6,6 @@
#include <vector>
#include <arbor/fvm_types.hpp>
-#include <arbor/ion_info.hpp>
#include "backends/gpu/gpu_store_types.hpp"
@@ -31,37 +30,33 @@ struct ion_state {
array eX_; // (mV) reversal potential
array Xi_; // (mM) internal concentration
array Xo_; // (mM) external concentration
- array weight_Xi_; // (1) concentration weight internal
- array weight_Xo_; // (1) concentration weight external
+
+ array init_Xi_; // (mM) area-weighted initial internal concentration
+ array init_Xo_; // (mM) area-weighted initial external concentration
+ array init_eX_; // (mM) initial reversal potential
array charge; // charge of ionic species (global, length 1)
- fvm_value_type default_int_concentration; // (mM) default internal concentration
- fvm_value_type default_ext_concentration; // (mM) default external concentration
ion_state() = default;
ion_state(
- ion_info info,
+ int charge,
const std::vector<fvm_index_type>& cv,
- const std::vector<fvm_value_type>& iconc_norm_area,
- const std::vector<fvm_value_type>& econc_norm_area,
+ const std::vector<fvm_value_type>& init_Xi,
+ const std::vector<fvm_value_type>& init_Xo,
+ const std::vector<fvm_value_type>& init_eX,
unsigned align
);
- // Calculate the reversal potential eX (mV) using Nernst equation
- void nernst(fvm_value_type temperature_K);
-
// Set ion concentrations to weighted proportion of default concentrations.
void init_concentration();
// Set ionic current density to zero.
void zero_current();
- void reset(fvm_value_type temperature_K) {
- zero_current();
- init_concentration();
- nernst(temperature_K);
- }
+ // Zero currents, reset concentrations, and reset reversal potential from
+ // initial values.
+ void reset();
};
struct shared_state {
@@ -78,7 +73,9 @@ struct shared_state {
array voltage; // Maps CV index to membrane voltage [mV].
array current_density; // Maps CV index to current density [A/m²].
array conductivity; // Maps CV index to membrane conductivity [kS/m²].
- array temperature_degC; // Global temperature [°C] (length 1 array).
+
+ array init_voltage; // Maps CV index to initial membrane voltage [mV].
+ array temperature_degC; // Maps CV to local temperature (read only) [°C].
std::unordered_map<std::string, ion_state> ion_data;
@@ -90,22 +87,23 @@ struct shared_state {
fvm_size_type n_intdom,
const std::vector<fvm_index_type>& cv_to_intdom_vec,
const std::vector<fvm_gap_junction>& gj_vec,
+ const std::vector<fvm_value_type>& init_membrane_potential,
+ const std::vector<fvm_value_type>& temperature_K,
unsigned align
);
void add_ion(
const std::string& ion_name,
- ion_info info,
+ int charge,
const std::vector<fvm_index_type>& cv,
- const std::vector<fvm_value_type>& iconc_norm_area,
- const std::vector<fvm_value_type>& econc_norm_area);
+ const std::vector<fvm_value_type>& init_iconc,
+ const std::vector<fvm_value_type>& init_econc,
+ const std::vector<fvm_value_type>& init_erev);
void zero_currents();
void ions_init_concentration();
- void ions_nernst_reversal_potential(fvm_value_type temperature_K);
-
// Set time_to to earliest of time+dt_step and tmax.
void update_time_to(fvm_value_type dt_step, fvm_value_type tmax);
@@ -128,7 +126,7 @@ struct shared_state {
array& sample_time,
array& sample_value);
- void reset(fvm_value_type initial_voltage, fvm_value_type temperature_K);
+ void reset();
};
// For debugging only
diff --git a/arbor/backends/multicore/mechanism.cpp b/arbor/backends/multicore/mechanism.cpp
index 117fdc25719bf66e288fe4d95f7f00a33dacbaed..5c62c1c3df4dba1ef6f05b80b7cbfb9ae083ff0f 100644
--- a/arbor/backends/multicore/mechanism.cpp
+++ b/arbor/backends/multicore/mechanism.cpp
@@ -93,7 +93,7 @@ void mechanism::instantiate(unsigned id, backend::shared_state& shared, const me
vec_i_ = shared.current_density.data();
vec_g_ = shared.conductivity.data();
- temperature_degC_ = &shared.temperature_degC;
+ temperature_degC_ = shared.temperature_degC.data();
auto ion_state_tbl = ion_state_table();
n_ion_ = ion_state_tbl.size();
diff --git a/arbor/backends/multicore/shared_state.cpp b/arbor/backends/multicore/shared_state.cpp
index 73c6605569e57dd08f6a20eaa47bc03052e3dad1..84986e53d14c33cc6590ef71746213c2c2b78bce 100644
--- a/arbor/backends/multicore/shared_state.cpp
+++ b/arbor/backends/multicore/shared_state.cpp
@@ -1,3 +1,4 @@
+#include <algorithm>
#include <cfloat>
#include <cmath>
#include <iostream>
@@ -10,7 +11,6 @@
#include <arbor/common_types.hpp>
#include <arbor/constants.hpp>
#include <arbor/fvm_types.hpp>
-#include <arbor/ion_info.hpp>
#include <arbor/math.hpp>
#include <arbor/simd/simd.hpp>
@@ -47,70 +47,44 @@ using pad = util::padded_allocator<>;
// ion_state methods:
ion_state::ion_state(
- ion_info info,
+ int charge,
const std::vector<fvm_index_type>& cv,
- const std::vector<fvm_value_type>& iconc_norm_area,
- const std::vector<fvm_value_type>& econc_norm_area,
+ const std::vector<fvm_value_type>& init_Xi,
+ const std::vector<fvm_value_type>& init_Xo,
+ const std::vector<fvm_value_type>& init_eX,
unsigned align
):
alignment(min_alignment(align)),
node_index_(cv.begin(), cv.end(), pad(alignment)),
iX_(cv.size(), NAN, pad(alignment)),
- eX_(cv.size(), NAN, pad(alignment)),
+ eX_(init_eX.begin(), init_eX.end(), pad(alignment)),
Xi_(cv.size(), NAN, pad(alignment)),
Xo_(cv.size(), NAN, pad(alignment)),
- weight_Xi_(iconc_norm_area.begin(), iconc_norm_area.end(), pad(alignment)),
- weight_Xo_(econc_norm_area.begin(), econc_norm_area.end(), pad(alignment)),
- charge(1u, info.charge, pad(alignment)),
- default_int_concentration(info.default_int_concentration),
- default_ext_concentration(info.default_ext_concentration)
+ init_Xi_(init_Xi.begin(), init_Xi.end(), pad(alignment)),
+ init_Xo_(init_Xo.begin(), init_Xo.end(), pad(alignment)),
+ init_eX_(init_eX.begin(), init_eX.end(), pad(alignment)),
+ charge(1u, charge, pad(alignment))
{
- arb_assert(node_index_.size()==weight_Xi_.size());
- arb_assert(node_index_.size()==weight_Xo_.size());
-}
-
-void ion_state::nernst(fvm_value_type temperature_K) {
- // Nernst equation: reversal potenial eX given by:
- //
- // eX = RT/zF * ln(Xo/Xi)
- //
- // where:
- // R: universal gas constant 8.3144598 J.K-1.mol-1
- // T: temperature in Kelvin
- // z: valency of species (K, Na: +1) (Ca: +2)
- // F: Faraday's constant 96485.33289 C.mol-1
- // Xo/Xi: ratio of out/in concentrations
-
- // 1e3 factor required to scale from V -> mV.
- constexpr fvm_value_type RF = 1e3*constant::gas_constant/constant::faraday;
-
- simd_value_type factor = RF*temperature_K/charge[0];
- for (std::size_t i=0; i<Xi_.size(); i+=simd_width) {
- simd_value_type xi(Xi_.data()+i);
- simd_value_type xo(Xo_.data()+i);
-
- auto ex = factor*log(xo/xi);
- ex.copy_to(eX_.data()+i);
- }
+ arb_assert(node_index_.size()==init_Xi_.size());
+ arb_assert(node_index_.size()==init_Xo_.size());
+ arb_assert(node_index_.size()==eX_.size());
+ arb_assert(node_index_.size()==init_eX_.size());
}
void ion_state::init_concentration() {
- for (std::size_t i=0u; i<Xi_.size(); i+=simd_width) {
- simd_value_type weight_xi(weight_Xi_.data()+i);
- simd_value_type weight_xo(weight_Xo_.data()+i);
-
- auto xi = default_int_concentration*weight_xi;
- xi.copy_to(Xi_.data()+i);
-
- auto xo = default_ext_concentration*weight_xo;
- xo.copy_to(Xo_.data()+i);
- }
+ std::copy(init_Xi_.begin(), init_Xi_.end(), Xi_.begin());
+ std::copy(init_Xo_.begin(), init_Xo_.end(), Xo_.begin());
}
void ion_state::zero_current() {
util::fill(iX_, 0);
}
+void ion_state::reset() {
+ zero_current();
+ init_concentration();
+ std::copy(init_eX_.begin(), init_eX_.end(), eX_.begin());
+}
// shared_state methods:
@@ -118,6 +92,8 @@ shared_state::shared_state(
fvm_size_type n_intdom,
const std::vector<fvm_index_type>& cv_to_intdom_vec,
const std::vector<fvm_gap_junction>& gj_vec,
+ const std::vector<fvm_value_type>& init_membrane_potential,
+ const std::vector<fvm_value_type>& temperature_K,
unsigned align
):
alignment(min_alignment(align)),
@@ -134,7 +110,8 @@ shared_state::shared_state(
voltage(n_cv, pad(alignment)),
current_density(n_cv, pad(alignment)),
conductivity(n_cv, pad(alignment)),
- temperature_degC(NAN),
+ init_voltage(init_membrane_potential.begin(), init_membrane_potential.end(), pad(alignment)),
+ temperature_degC(n_cv, pad(alignment)),
deliverable_events(n_intdom)
{
// For indices in the padded tail of cv_to_intdom, set index to last valid intdom index.
@@ -146,30 +123,34 @@ shared_state::shared_state(
std::copy(gj_vec.begin(), gj_vec.end(), gap_junctions.begin());
std::fill(gap_junctions.begin()+n_gj, gap_junctions.end(), gj_vec.back());
}
+
+ for (unsigned i = 0; i<n_cv; ++i) {
+ temperature_degC[i] = temperature_K[i] - 273.15;
+ }
}
void shared_state::add_ion(
const std::string& ion_name,
- ion_info info,
+ int charge,
const std::vector<fvm_index_type>& cv,
- const std::vector<fvm_value_type>& iconc_norm_area,
- const std::vector<fvm_value_type>& econc_norm_area)
+ const std::vector<fvm_value_type>& init_iconc,
+ const std::vector<fvm_value_type>& init_econc,
+ const std::vector<fvm_value_type>& init_erev)
{
ion_data.emplace(std::piecewise_construct,
std::forward_as_tuple(ion_name),
- std::forward_as_tuple(info, cv, iconc_norm_area, econc_norm_area, alignment));
+ std::forward_as_tuple(charge, cv, init_iconc, init_econc, init_erev, alignment));
}
-void shared_state::reset(fvm_value_type initial_voltage, fvm_value_type temperature_K) {
- util::fill(voltage, initial_voltage);
+void shared_state::reset() {
+ std::copy(init_voltage.begin(), init_voltage.end(), voltage.begin());
util::fill(current_density, 0);
util::fill(conductivity, 0);
util::fill(time, 0);
util::fill(time_to, 0);
- temperature_degC = temperature_K - 273.15;
for (auto& i: ion_data) {
- i.second.reset(temperature_K);
+ i.second.reset();
}
}
@@ -187,12 +168,6 @@ void shared_state::ions_init_concentration() {
}
}
-void shared_state::ions_nernst_reversal_potential(fvm_value_type temperature_K) {
- for (auto& i: ion_data) {
- i.second.nernst(temperature_K);
- }
-}
-
void shared_state::update_time_to(fvm_value_type dt_step, fvm_value_type tmax) {
for (fvm_size_type i = 0; i<n_intdom; i+=simd_width) {
simd_value_type t(time.data()+i);
@@ -258,23 +233,28 @@ std::ostream& operator<<(std::ostream& out, const shared_state& s) {
using io::csv;
out << "n_intdom " << s.n_intdom << "\n";
- out << "n_cv " << s.n_cv << "\n";
+ out << "n_cv " << s.n_cv << "\n";
out << "cv_to_intdom " << csv(s.cv_to_intdom) << "\n";
- out << "time " << csv(s.time) << "\n";
- out << "time_to " << csv(s.time_to) << "\n";
+ out << "time " << csv(s.time) << "\n";
+ out << "time_to " << csv(s.time_to) << "\n";
out << "dt_intdom " << csv(s.dt_intdom) << "\n";
- out << "dt_cv " << csv(s.dt_cv) << "\n";
- out << "voltage " << csv(s.voltage) << "\n";
- out << "current " << csv(s.current_density) << "\n";
+ out << "dt_cv " << csv(s.dt_cv) << "\n";
+ out << "voltage " << csv(s.voltage) << "\n";
+ out << "init_voltage " << csv(s.init_voltage) << "\n";
+ out << "temperature " << csv(s.temperature_degC) << "\n";
+ out << "current " << csv(s.current_density) << "\n";
out << "conductivity " << csv(s.conductivity) << "\n";
for (const auto& ki: s.ion_data) {
auto& kn = ki.first;
auto& i = const_cast<ion_state&>(ki.second);
- out << kn << ".current_density " << csv(i.iX_) << "\n";
- out << kn << ".reversal_potential " << csv(i.eX_) << "\n";
- out << kn << ".internal_concentration " << csv(i.Xi_) << "\n";
- out << kn << ".external_concentration " << csv(i.Xo_) << "\n";
- out << kn << ".node_index " << csv(i.node_index_) << "\n";
+ out << kn << "/current_density " << csv(i.iX_) << "\n";
+ out << kn << "/reversal_potential " << csv(i.eX_) << "\n";
+ out << kn << "/internal_concentration " << csv(i.Xi_) << "\n";
+ out << kn << "/external_concentration " << csv(i.Xo_) << "\n";
+ out << kn << "/intconc_initial " << csv(i.init_Xi_) << "\n";
+ out << kn << "/extconc_initial " << csv(i.init_Xo_) << "\n";
+ out << kn << "/revpot_initial " << csv(i.init_eX_) << "\n";
+ out << kn << "/node_index " << csv(i.node_index_) << "\n";
}
return out;
diff --git a/arbor/backends/multicore/shared_state.hpp b/arbor/backends/multicore/shared_state.hpp
index 4d97d1cd0d3d3350c4e9ab8e449fca8b91b3eb5a..3a16c03a8ddaec0b8cae780080c1d5b3146eceac 100644
--- a/arbor/backends/multicore/shared_state.hpp
+++ b/arbor/backends/multicore/shared_state.hpp
@@ -10,7 +10,6 @@
#include <arbor/assert.hpp>
#include <arbor/common_types.hpp>
#include <arbor/fvm_types.hpp>
-#include <arbor/ion_info.hpp>
#include <arbor/simd/simd.hpp>
#include "backends/event.hpp"
@@ -46,37 +45,33 @@ struct ion_state {
array eX_; // (mV) reversal potential
array Xi_; // (mM) internal concentration
array Xo_; // (mM) external concentration
- array weight_Xi_; // (1) concentration weight internal
- array weight_Xo_; // (1) concentration weight external
+
+ array init_Xi_; // (mM) area-weighted initial internal concentration
+ array init_Xo_; // (mM) area-weighted initial external concentration
+ array init_eX_; // (mV) initial reversal potential
array charge; // charge of ionic species (global value, length 1)
- fvm_value_type default_int_concentration; // (mM) default internal concentration
- fvm_value_type default_ext_concentration; // (mM) default external concentration
ion_state() = default;
ion_state(
- ion_info info,
+ int charge,
const std::vector<fvm_index_type>& cv,
- const std::vector<fvm_value_type>& iconc_norm_area,
- const std::vector<fvm_value_type>& econc_norm_area,
+ const std::vector<fvm_value_type>& init_Xi,
+ const std::vector<fvm_value_type>& init_Xo,
+ const std::vector<fvm_value_type>& init_eX,
unsigned align
);
- // Calculate the reversal potential eX (mV) using Nernst equation
- void nernst(fvm_value_type temperature_K);
-
// Set ion concentrations to weighted proportion of default concentrations.
void init_concentration();
// Set ionic current density to zero.
void zero_current();
- void reset(fvm_value_type temperature_K) {
- zero_current();
- init_concentration();
- nernst(temperature_K);
- }
+ // Zero currents, reset concentrations, and reset reversal potential from
+ // initial values.
+ void reset();
};
struct shared_state {
@@ -96,7 +91,9 @@ struct shared_state {
array voltage; // Maps CV index to membrane voltage [mV].
array current_density; // Maps CV index to membrane current density contributions [A/m²].
array conductivity; // Maps CV index to membrane conductivity [kS/m²].
- fvm_value_type temperature_degC; // Global temperature [°C].
+
+ array init_voltage; // Maps CV index to initial membrane voltage [mV].
+ array temperature_degC; // Maps CV to local temperature (read only) [°C].
std::unordered_map<std::string, ion_state> ion_data;
@@ -108,15 +105,18 @@ struct shared_state {
fvm_size_type n_intdom,
const std::vector<fvm_index_type>& cv_to_intdom_vec,
const std::vector<fvm_gap_junction>& gj_vec,
+ const std::vector<fvm_value_type>& init_membrane_potential,
+ const std::vector<fvm_value_type>& temperature_K,
unsigned align
);
void add_ion(
const std::string& ion_name,
- ion_info info,
+ int charge,
const std::vector<fvm_index_type>& cv,
- const std::vector<fvm_value_type>& iconc_norm_area,
- const std::vector<fvm_value_type>& econc_norm_area);
+ const std::vector<fvm_value_type>& init_iconc,
+ const std::vector<fvm_value_type>& init_econc,
+ const std::vector<fvm_value_type>& init_erev);
void zero_currents();
@@ -146,7 +146,7 @@ struct shared_state {
array& sample_time,
array& sample_value);
- void reset(fvm_value_type initial_voltage, fvm_value_type temperature_K);
+ void reset();
};
// For debugging only:
diff --git a/arbor/cable_cell.cpp b/arbor/cable_cell.cpp
index 7403dc1e7ab5f7e9aeb3dc066abf3bcc0d87b7be..bea58ad8b7e86b35ba31dfa9908e7e379e307063 100644
--- a/arbor/cable_cell.cpp
+++ b/arbor/cable_cell.cpp
@@ -85,6 +85,12 @@ cable_segment* cable_cell::cable(index_type index) {
return cable? cable: throw cable_cell_error("segment is not a cable segment");
}
+const cable_segment* cable_cell::cable(index_type index) const {
+ assert_valid_segment(index);
+ auto cable = segment(index)->as_cable();
+ return cable? cable: throw cable_cell_error("segment is not a cable segment");
+}
+
std::vector<size_type> cable_cell::compartment_counts() const {
std::vector<size_type> comp_count;
comp_count.reserve(num_segments());
@@ -108,6 +114,55 @@ void cable_cell::add_detector(segment_location loc, double threshold) {
spike_detectors_.push_back({loc, threshold});
}
+// Approximating wildly by ignoring O(x) effects entirely, the attenuation b
+// over a single cable segment with constant resistivity R and membrane
+// capacitance C is given by:
+//
+// b = 2√(Ï€RCf) · Σ 2L/(√dâ‚€ + √dâ‚)
+//
+// where the sum is taken over each piecewise linear segment of length L
+// with diameters dâ‚€ and dâ‚ at each end.
+
+value_type cable_cell::segment_mean_attenuation(
+ value_type frequency, index_type segidx,
+ const cable_cell_parameter_set& global_defaults) const
+{
+ value_type R = default_parameters.axial_resistivity.value_or(global_defaults.axial_resistivity.value());
+ value_type C = default_parameters.membrane_capacitance.value_or(global_defaults.membrane_capacitance.value());
+
+ value_type length_factor = 0; // [1/õm]
+
+ if (segidx==0) {
+ if (const soma_segment* s = soma()) {
+ R = s->parameters.axial_resistivity.value_or(R);
+ C = s->parameters.membrane_capacitance.value_or(C);
+
+ value_type d = 2*s->radius();
+ length_factor = 1/std::sqrt(d);
+ }
+ }
+ else {
+ const cable_segment* s = cable(segidx);
+ const auto& lengths = s->lengths();
+ const auto& radii = s->radii();
+
+ value_type total_length = 0;
+ R = s->parameters.axial_resistivity.value_or(R);
+ C = s->parameters.membrane_capacitance.value_or(C);
+
+ for (std::size_t i = 0; i<lengths.size(); ++i) {
+ length_factor += 2*lengths[i]/(std::sqrt(radii[i])+std::sqrt(radii[i+1]));
+ total_length += lengths[i];
+ }
+ length_factor /= total_length;
+ }
+
+ // R*C is in [s·cm/m²]; need to convert to [s/µm]
+ value_type tau_per_um = R*C*1e-8;
+
+ return 2*std::sqrt(math::pi<double>*tau_per_um*frequency)*length_factor; // [1/µm]
+}
+
// Construct cell from flat morphology specification.
diff --git a/arbor/cable_cell_param.cpp b/arbor/cable_cell_param.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..a60bcfb2a2facd822240752a7cb242d0629dd62b
--- /dev/null
+++ b/arbor/cable_cell_param.cpp
@@ -0,0 +1,67 @@
+#include <cfloat>
+#include <cmath>
+
+#include <arbor/cable_cell_param.hpp>
+
+#include "util/maputil.hpp"
+
+namespace arb {
+
+void check_global_properties(const cable_cell_global_properties& G) {
+ auto& param = G.default_parameters;
+
+ if (!param.init_membrane_potential) {
+ throw cable_cell_error("missing global default parameter value: init_membrane_potential");
+ }
+
+ if (!param.temperature_K) {
+ throw cable_cell_error("missing global default parameter value: temperature");
+ }
+
+ if (!param.axial_resistivity) {
+ throw cable_cell_error("missing global default parameter value: axial_resistivity");
+ }
+
+ if (!param.membrane_capacitance) {
+ throw cable_cell_error("missing global default parameter value: membrane_capacitance");
+ }
+
+ for (const auto& ion: util::keys(G.ion_species)) {
+ if (!param.ion_data.count(ion)) {
+ throw cable_cell_error("missing ion defaults for ion "+ion);
+ }
+ }
+
+ for (const auto& kv: param.ion_data) {
+ auto& ion = kv.first;
+ const cable_cell_ion_data& data = kv.second;
+ if (std::isnan(data.init_int_concentration)) {
+ throw cable_cell_error("missing init_int_concentration for ion "+ion);
+ }
+ if (std::isnan(data.init_ext_concentration)) {
+ throw cable_cell_error("missing init_ext_concentration for ion "+ion);
+ }
+ if (std::isnan(data.init_reversal_potential) && !param.reversal_potential_method.count(ion)) {
+ throw cable_cell_error("missing init_reversal_potential or reversal_potential_method for ion "+ion);
+ }
+ }
+}
+
+cable_cell_local_parameter_set neuron_parameter_defaults = {
+ // ion defaults: internal concentration [mM], external concentration [mM], reversal potential [mV]
+ {
+ {"na", {10.0, 140.0, 115 - 65.}},
+ {"k", {54.4, 2.5, -12 - 65.}},
+ {"ca", {5e-5, 2.0, 12.5*std::log(2.0/5e-5)}}
+ },
+ // initial membrane potential [mV]
+ -65.0,
+ // temperatue [K]
+ 6.3 + 273.15,
+ // axial resistivity [Ω·cm]
+ 35.4,
+ // membrane capacitance [F/m²]
+ 0.01
+};
+
+} // namespace arb
diff --git a/arbor/fvm_layout.cpp b/arbor/fvm_layout.cpp
index dd4d5b51042516ff61153b06ccbe808ffaf4a07d..f035b70b88f8b2290a281a8af1f5a3ce4d08d8c8 100644
--- a/arbor/fvm_layout.cpp
+++ b/arbor/fvm_layout.cpp
@@ -5,6 +5,7 @@
#include <arbor/arbexcept.hpp>
#include <arbor/cable_cell.hpp>
+#include <arbor/util/optional.hpp>
#include "algorithms.hpp"
#include "fvm_compartment.hpp"
@@ -19,7 +20,9 @@
namespace arb {
using util::count_along;
+using util::keys;
using util::make_span;
+using util::optional;
using util::subrange_view;
using util::transform_view;
using util::value_by_key;
@@ -34,6 +37,14 @@ namespace {
}
}
+ template <typename ResizableContainer, typename Index>
+ auto& extend_at(ResizableContainer& c, const Index& i) {
+ if (util::size(c)<=i) {
+ c.resize(i+1);
+ }
+ return c[i];
+ }
+
struct compartment_model {
arb::tree tree;
std::vector<tree::int_type> parent_index;
@@ -59,6 +70,73 @@ namespace {
};
ARB_DEFINE_LEXICOGRAPHIC_ORDERING(cv_param,(a.cv,a.params,a.target),(b.cv,b.params,b.target))
+
+ template <typename V>
+ optional<V> operator|(const optional<V>& a, const optional<V>& b) {
+ return a? a: b;
+ }
+
+ // For each segment given by the provided sorted sequence of segment
+ // indices, call `action` for each CV intersecting the segment, starting
+ // from the most proximal.
+ //
+ // By the ordering of CVs and segments in the discretization, with the
+ // exception of the most proximal CV in a segment, each CV will be visited
+ // once, and the visited CVs will be in increasing order. The most proximal
+ // CV (the 'parent' CV) may be visited multiple times.
+ //
+ // Action is a functional that takes the following arguments:
+ //
+ // size_type cv_index The index into the total (sorted) list of
+ // CVs that constitute the segments.
+ //
+ // index_type cv The CV number (within the discretization).
+ //
+ // value_type cv_area The area of the CV that lies within the
+ // current segment.
+ //
+ // index_type seg_index The index into the provided sequence of
+ // the provided segment_indices.
+ //
+ // index_type seg The segment currently being iterated over.
+
+ template <typename Seq, typename Action>
+ void for_each_cv_in_segments(const fvm_discretization& D, const Seq& segment_indices, const Action& action) {
+ using index_type = fvm_index_type;
+ using size_type = fvm_size_type;
+
+ std::unordered_map<index_type, size_type> parent_cv_indices;
+ size_type cv_index = 0;
+
+ index_type seg_index = 0;
+ for (const auto& seg: segment_indices) {
+ const segment_info& seg_info = D.segments[seg];
+
+ if (seg_info.has_parent()) {
+ index_type cv = seg_info.parent_cv;
+
+ size_type i = parent_cv_indices.insert({cv, cv_index}).first->second;
+ if (i==cv_index) {
+ ++cv_index;
+ }
+
+ action(i, cv, seg_info.parent_cv_area, seg_index, seg);
+ }
+
+ for (index_type cv = seg_info.proximal_cv; cv < seg_info.distal_cv; ++cv) {
+ index_type i = cv_index++;
+ action(i, cv, D.cv_area[cv], seg_index, seg);
+ }
+
+ index_type cv = seg_info.distal_cv;
+ size_type i = cv_index++;
+
+ action(i, seg_info.distal_cv, seg_info.distal_cv_area, seg_index, seg);
+ parent_cv_indices.insert({cv, i});
+ ++seg_index;
+ }
+ }
+
} // namespace
// Cable segment discretization
@@ -123,7 +201,7 @@ namespace {
// = 1/R · hVâ‚Vâ‚‚/(h₂²Vâ‚+h₲Vâ‚‚)
//
-fvm_discretization fvm_discretize(const std::vector<cable_cell>& cells) {
+fvm_discretization fvm_discretize(const std::vector<cable_cell>& cells, const cable_cell_parameter_set& global_defaults) {
using value_type = fvm_value_type;
using index_type = fvm_index_type;
@@ -153,6 +231,8 @@ fvm_discretization fvm_discretize(const std::vector<cable_cell>& cells) {
D.face_conductance.assign(D.ncv, 0.);
D.cv_area.assign(D.ncv, 0.);
D.cv_capacitance.assign(D.ncv, 0.);
+ D.init_membrane_potential.assign(D.ncv, 0.);
+ D.temperature_K.assign(D.ncv, 0.);
D.parent_cv.assign(D.ncv, index_type(-1));
D.cv_to_cell.resize(D.ncv);
for (auto i: make_span(0, D.ncell)) {
@@ -170,6 +250,12 @@ fvm_discretization fvm_discretize(const std::vector<cable_cell>& cells) {
D.parent_cv[k] = cell_graph.parent_index[k-cell_cv_base]+cell_cv_base;
}
+ // Electrical defaults from global defaults, possibly overridden by cell.
+ auto cm_default = c.default_parameters.membrane_capacitance | global_defaults.membrane_capacitance;
+ auto rL_default = c.default_parameters.axial_resistivity | global_defaults.axial_resistivity;
+ auto init_vm_default = c.default_parameters.init_membrane_potential | global_defaults.init_membrane_potential;
+ auto temp_default = c.default_parameters.temperature_K | global_defaults.temperature_K;
+
// Compartment index range for each segment in this cell.
seg_cv_bounds.clear();
auto seg_cv_part = make_partition(
@@ -206,6 +292,9 @@ fvm_discretization fvm_discretize(const std::vector<cable_cell>& cells) {
soma_info.distal_cv_area = soma_area;
D.segments.push_back(soma_info);
+ index_type soma_segment_index = D.segments.size()-1;
+ D.parent_segment.push_back(soma_segment_index);
+
// Other segments must all be cable segments.
for (size_type j = 1; j<nseg; ++j) {
const auto& seg_cv_ival = seg_cv_part[j];
@@ -217,8 +306,12 @@ fvm_discretization fvm_discretize(const std::vector<cable_cell>& cells) {
if (!cable) {
throw arbor_internal_error("fvm_layout: non-root segments of cell must be cable segments");
}
- auto cm = cable->cm; // [F/m²]
- auto rL = cable->rL; // [Ω·cm]
+
+ const auto& params = cable->parameters;
+ auto cm = (params.membrane_capacitance | cm_default).value(); // [F/m²]
+ auto rL = (params.axial_resistivity | rL_default).value(); // [Ω·cm]
+ auto init_vm = (params.init_membrane_potential | init_vm_default).value(); // [mV]
+ auto temp = (params.temperature_K | temp_default).value(); // [mV]
bool soma_parent = c.parent(j)->as_soma() ? true : false; //segment's parent is a soma
@@ -226,7 +319,7 @@ fvm_discretization fvm_discretize(const std::vector<cable_cell>& cells) {
auto lengths = cable->lengths();
// If segment has soma parent, send soma information to div_compartment_integrator
- if(soma_parent) {
+ if (soma_parent) {
radii.insert(radii.begin(), soma->radius());
lengths.insert(lengths.begin(), soma->radius()*2);
}
@@ -242,6 +335,19 @@ fvm_discretization fvm_discretize(const std::vector<cable_cell>& cells) {
seg_info.distal_cv_area = divs(ncv-1).right.area;
D.segments.push_back(seg_info);
+ if (soma_parent) {
+ D.parent_segment.push_back(soma_segment_index);
+ }
+ else {
+ auto opt_index = util::binary_search_index(D.segments, seg_info.parent_cv,
+ [&D](const segment_info& seg_info) { return seg_info.distal_cv; });
+
+ if (!opt_index) {
+ throw arbor_internal_error("fvm_layout: could not find parent segment");
+ }
+ D.parent_segment.push_back(*opt_index);
+ }
+
for (auto i: make_span(seg_cv_ival)) {
const auto& div = divs(i-seg_cv_ival.first);
@@ -261,14 +367,33 @@ fvm_discretization fvm_discretize(const std::vector<cable_cell>& cells) {
auto ar = div.right.area; // [µm²]
D.cv_area[j] += al; // [µm²]
- D.cv_capacitance[j] += al * cm; // [pF]
+ D.cv_capacitance[j] += al*cm; // [pF]
+ D.init_membrane_potential[j] += al*init_vm; // [mV·µm²]
+ D.temperature_K[j] += al*temp; // [K·µm²]
+
D.cv_area[i] += ar; // [µm²]
- D.cv_capacitance[i] += ar * cm; // [pF]
+ D.cv_capacitance[i] += ar*cm; // [pF]
+ D.init_membrane_potential[i] += ar*init_vm; // [mV·µm²]
+ D.temperature_K[i] += ar*temp; // [K·µm²]
}
}
+ auto soma_cm = (soma->parameters.membrane_capacitance | cm_default).value(); // [F/m²]
+ auto soma_init_vm = (soma->parameters.init_membrane_potential | init_vm_default).value(); // [mV]
+ auto soma_temp = (soma->parameters.temperature_K | temp_default).value(); // [mV]
+
D.cv_area[soma_cv] = soma_area; // [µm²]
- D.cv_capacitance[soma_cv] = soma_area*soma->cm; // [pF]
+ D.cv_capacitance[soma_cv] = soma_area*soma_cm; // [pF]
+ D.init_membrane_potential[soma_cv] = soma_area*soma_init_vm; // [mV·µm²]
+ D.temperature_K[soma_cv] = soma_area*soma_temp; // [K·µm²]
+ }
+
+ // Rescale CV init_vm and temperature values to get area-weighted means.
+ for (auto i: make_span(0, D.ncv)) {
+ if (D.cv_area[i]) {
+ D.init_membrane_potential[i] /= D.cv_area[i]; // [mV]
+ D.temperature_K[i] /= D.cv_area[i]; // [mV]
+ }
}
// Number of CVs per cell is exactly number of compartments.
@@ -301,6 +426,7 @@ fvm_mechanism_data fvm_build_mechanism_data(const cable_cell_global_properties&
using string_index_map = std::unordered_map<std::string, size_type>;
const mechanism_catalogue& catalogue = *gprop.catalogue;
+ const cable_cell_parameter_set& gparam = gprop.default_parameters;
fvm_mechanism_data mechdata;
@@ -312,12 +438,24 @@ fvm_mechanism_data fvm_build_mechanism_data(const cable_cell_global_properties&
// 3. Pointer to mechanism metadata from catalogue.
struct density_mech_data {
- std::vector<std::pair<size_type, const mechanism_desc*>> segments; //
+ std::vector<std::pair<size_type, const mechanism_desc*>> segments;
string_set paramset;
std::shared_ptr<mechanism_info> info = nullptr;
};
std::unordered_map<std::string, density_mech_data> density_mech_table;
+ // Temporary table for revpot mechanism info, mapping mechanism name to tuple of:
+ // 1. Sorted map from cell index to mechanism parameter settings.
+ // 2. Set of the names of parameters that are anywhere modified.
+ // 3. Pointer to mechanism metadat from catalogue.
+
+ struct revpot_mech_data {
+ std::map<size_type, const mechanism_desc*> cells;
+ string_set paramset;
+ std::shared_ptr<mechanism_info> info = nullptr;
+ };
+ std::unordered_map<std::string, revpot_mech_data> revpot_mech_table;
+
// Temporary table for point mechanism info, mapping mechanism name to tuple:
// 1. Vector of point info: CV, index into cell group targets, parameter settings.
// 2. Set of the names of parameters that are anywhere modified.
@@ -337,38 +475,34 @@ fvm_mechanism_data fvm_build_mechanism_data(const cable_cell_global_properties&
// Built-in stimulus mechanism data is dealt with especially below.
// Record for each stimulus the CV and clamp data.
+
std::vector<std::pair<size_type, i_clamp>> stimuli;
- // Temporary table for presence of ion channels, mapping ion name to _sorted_
- // collection of segment indices.
+ // Temporary table for presence of ion channels, mapping ion name to a _sorted_
+ // collection of per-segment ion data, viz. a map from segment index to
+ // initial internal and external concentrations and reversal potentials, plus
+ // information about whether there is a mechanism that requires this ion
+ // on this segment, and if that ion writes to internal or external concentration.
+
+ struct ion_segment_data {
+ cable_cell_ion_data ion_data;
+ bool mech_requires = false;
+ bool mech_writes_iconc = false;
+ bool mech_writes_econc = false;
+ };
- std::unordered_map<std::string, std::set<size_type>> ion_segments;
+ std::unordered_map<std::string, std::map<index_type, ion_segment_data>> ion_segments;
- auto update_paramset_and_validate =
- [&catalogue]
- (const mechanism_desc& desc, std::shared_ptr<mechanism_info>& info, string_set& paramset)
- {
- auto& name = desc.name();
- if (!info) {
- info.reset(new mechanism_info(catalogue[name]));
- }
- for (const auto& pv: desc.values()) {
- if (!paramset.count(pv.first)) {
- if (!info->parameters.count(pv.first)) {
- throw no_such_parameter(name, pv.first);
- }
- if (!info->parameters.at(pv.first).valid(pv.second)) {
- throw invalid_parameter_value(name, pv.first, pv.second);
- }
- paramset.insert(pv.first);
- }
- }
- };
+ // Temporary table for presence of mechanisms that read the reversal potential
+ // of an ion channel, mapping ion name and cell index to a _sorted_
+ // collection of segment indices.
+
+ std::unordered_map<std::string, std::unordered_map<size_type, std::set<size_type>>> ion_revpot_segments;
auto verify_ion_usage =
[&gprop](const std::string& mech_name, const mechanism_info* info)
{
- const auto& global_ions = gprop.ion_default;
+ const auto& global_ions = gprop.ion_species;
for (const auto& ion: info->ions) {
const auto& ion_name = ion.first;
@@ -380,7 +514,7 @@ fvm_mechanism_data fvm_build_mechanism_data(const cable_cell_global_properties&
}
if (ion_dep.verify_ion_charge) {
- if (ion_dep.expected_ion_charge!=global_ions.at(ion_name).charge) {
+ if (ion_dep.expected_ion_charge!=global_ions.at(ion_name)) {
throw cable_cell_error(
"mechanism "+mech_name+" uses ion "+ion_name+ " expecting a different valence");
}
@@ -388,6 +522,28 @@ fvm_mechanism_data fvm_build_mechanism_data(const cable_cell_global_properties&
}
};
+ auto update_paramset_and_validate =
+ [&catalogue,&verify_ion_usage]
+ (const mechanism_desc& desc, std::shared_ptr<mechanism_info>& info, string_set& paramset)
+ {
+ auto& name = desc.name();
+ if (!info) {
+ info.reset(new mechanism_info(catalogue[name]));
+ verify_ion_usage(name, info.get());
+ }
+ for (const auto& pv: desc.values()) {
+ if (!paramset.count(pv.first)) {
+ if (!info->parameters.count(pv.first)) {
+ throw no_such_parameter(name, pv.first);
+ }
+ if (!info->parameters.at(pv.first).valid(pv.second)) {
+ throw invalid_parameter_value(name, pv.first, pv.second);
+ }
+ paramset.insert(pv.first);
+ }
+ }
+ };
+
auto cell_segment_part = D.cell_segment_part();
size_type target_id = 0;
@@ -395,16 +551,61 @@ fvm_mechanism_data fvm_build_mechanism_data(const cable_cell_global_properties&
auto& cell = cells[cell_idx];
auto seg_range = cell_segment_part[cell_idx];
+ auto add_ion_segment =
+ [&gparam, &cell, &ion_segments]
+ (const std::string& ion_name, size_type segment_idx, const cable_cell_local_parameter_set& seg_param, const ion_dependency* iondep = nullptr)
+ {
+ const auto& global_ion_data = gparam.ion_data;
+ const auto& cell_ion_data = cell.default_parameters.ion_data;
+ const auto& seg_ion_data = seg_param.ion_data;
+
+ auto& ion_entry = ion_segments[ion_name];
+
+ // New entry?
+ util::optional<cable_cell_ion_data> opt_ion_data;
+ if (!ion_entry.count(segment_idx)) {
+ opt_ion_data = value_by_key(seg_ion_data, ion_name);
+ if (!opt_ion_data) {
+ opt_ion_data = value_by_key(cell_ion_data, ion_name);
+ }
+ if (!opt_ion_data) {
+ opt_ion_data = value_by_key(global_ion_data, ion_name);
+ }
+ if (!opt_ion_data) {
+ throw arbor_internal_error("missing entry for ion "+ion_name+" in cable_cell global defaults");
+ }
+ }
+
+ auto& ion_entry_seg = ion_entry[segment_idx];
+ if (opt_ion_data) {
+ ion_entry_seg.ion_data = *opt_ion_data;
+ }
+ if (iondep) {
+ ion_entry_seg.mech_requires = true;
+ ion_entry_seg.mech_writes_iconc |= iondep->write_concentration_int;
+ ion_entry_seg.mech_writes_econc |= iondep->write_concentration_ext;
+ }
+ };
+
for (auto segment_idx: make_span(seg_range)) {
- for (const mechanism_desc& desc: cell.segments()[segment_idx-seg_range.first]->mechanisms()) {
+ const segment_ptr& seg = cell.segments()[segment_idx-seg_range.first];
+
+ for (const mechanism_desc& desc: seg->mechanisms()) {
const auto& name = desc.name();
density_mech_data& entry = density_mech_table[name];
update_paramset_and_validate(desc, entry.info, entry.paramset);
entry.segments.emplace_back(segment_idx, &desc);
- for (const auto& ion: entry.info->ions) {
- ion_segments[ion.first].insert(segment_idx);
+ for (const auto& ion_entry: entry.info->ions) {
+ const std::string& ion_name = ion_entry.first;
+ const ion_dependency& iondep = ion_entry.second;
+
+ add_ion_segment(ion_name, segment_idx, seg->parameters, &iondep);
+
+ if (ion_entry.second.read_reversal_potential) {
+ ion_revpot_segments[ion_name][cell_idx].insert(segment_idx);
+ }
}
}
}
@@ -418,9 +619,18 @@ fvm_mechanism_data fvm_build_mechanism_data(const cable_cell_global_properties&
update_paramset_and_validate(desc, entry.info, entry.paramset);
entry.points.push_back({cv, target_id++, &desc});
+ const segment_ptr& seg = cell.segments()[cellsyn.location.segment];
size_type segment_idx = D.cell_segment_bounds[cell_idx]+cellsyn.location.segment;
- for (const auto& ion: entry.info->ions) {
- ion_segments[ion.first].insert(segment_idx);
+
+ for (const auto& ion_entry: entry.info->ions) {
+ const std::string& ion_name = ion_entry.first;
+ const ion_dependency& iondep = ion_entry.second;
+
+ add_ion_segment(ion_name, segment_idx, seg->parameters, &iondep);
+
+ if (ion_entry.second.read_reversal_potential) {
+ ion_revpot_segments[ion_name][cell_idx].insert(segment_idx);
+ }
}
}
@@ -428,22 +638,85 @@ fvm_mechanism_data fvm_build_mechanism_data(const cable_cell_global_properties&
size_type cv = D.segment_location_cv(cell_idx, stimulus.location);
stimuli.push_back({cv, stimulus.clamp});
}
- }
- for (auto& entry: density_mech_table) {
- verify_ion_usage(entry.first, entry.second.info.get());
- }
+ // Add segments to ion_segments map which intersect with existing segments, so
+ // that each CV with an ion value is 100% covered.
+
+ for (auto& e: ion_segments) {
+ const auto& name = e.first;
+ const auto& ion_entry = e.second;
+
+ for (auto segment_idx: make_span(seg_range)) {
+ index_type parent_segment_idx = D.parent_segment[segment_idx];
+ const segment_ptr& parent_seg = cell.segments()[parent_segment_idx-seg_range.first];
+
+ if (ion_entry.count(segment_idx)) {
+ add_ion_segment(name, parent_segment_idx, parent_seg->parameters);
+ }
+ }
+
+ for (auto segment_idx: make_span(seg_range)) {
+ const segment_ptr& seg = cell.segments()[segment_idx-seg_range.first];
+ index_type parent_segment_idx = D.parent_segment[segment_idx];
+
+ if (ion_entry.count(parent_segment_idx)) {
+ add_ion_segment(name, segment_idx, seg->parameters);
+ }
+ }
+ }
- for (auto& entry: point_mech_table) {
- verify_ion_usage(entry.first, entry.second.info.get());
+
+ // Maintain a map of reversal potential mechanism to written ions for this
+ // cell, to ensure consistency with assignments from the cell and global
+ // parameters.
+
+ std::unordered_multimap<std::string, std::string> revpot_to_ion;
+ std::unordered_map<std::string, std::string> ion_to_revpot;
+
+ auto add_revpot = [&](const std::string& ion, const mechanism_desc& desc) {
+ const auto& name = desc.name();
+
+ revpot_mech_data& entry = revpot_mech_table[name];
+ update_paramset_and_validate(desc, entry.info, entry.paramset);
+ ion_to_revpot[ion] = desc.name();
+ for (auto dep: entry.info->ions) {
+ if (dep.second.write_reversal_potential) {
+ revpot_to_ion.insert({desc.name(), dep.first});
+ }
+ }
+ entry.cells.insert({cell_idx, &desc});
+ };
+
+ const auto& cellrevpot = cell.default_parameters.reversal_potential_method;
+ for (const auto& revpot: cellrevpot) {
+ add_revpot(revpot.first, revpot.second);
+ }
+
+ const auto& globalrevpot = gparam.reversal_potential_method;
+ for (const auto& revpot: globalrevpot) {
+ if (!cellrevpot.count(revpot.first)) {
+ add_revpot(revpot.first, revpot.second);
+ }
+ }
+
+ // Ensure that if a revpot mechanism writes to multiple ions, that
+ // that mechanism is associated with all of them.
+ for (auto& entry: revpot_to_ion) {
+ auto declared_revpot = value_by_key(ion_to_revpot, entry.second);
+ if (!declared_revpot || declared_revpot.value()!=entry.first) {
+ throw cable_cell_error("mechanism "+entry.first+" writes reversal potential for "
+ +entry.second+", which is not configured to use "+entry.first);
+ }
+ }
}
+
// II. Build ion and mechanism configs.
// Shared temporary lookup info across mechanism instances, set by build_param_data.
- string_index_map param_index;
- std::vector<std::string> param_name;
- std::vector<value_type> param_default;
+ string_index_map param_index; // maps parameter name to parameter index
+ std::vector<std::string> param_name; // indexed by parameter index
+ std::vector<value_type> param_default; // indexed by parameter index
auto build_param_data =
[¶m_name, ¶m_index, ¶m_default](const string_set& paramset, const mechanism_info* info)
@@ -463,33 +736,109 @@ fvm_mechanism_data fvm_build_mechanism_data(const cable_cell_global_properties&
// IIa. Ion channel CVs.
- for (auto& ionseg: ion_segments) {
- auto& ion = mechdata.ions[ionseg.first];
+ for (const auto& ionseg: ion_segments) {
+ auto& seg_ion_map = ion_segments[ionseg.first];
+ fvm_ion_config& ion_config = mechdata.ions[ionseg.first];
+ auto& seg_ion_data = ionseg.second;
+
+ for_each_cv_in_segments(D, keys(seg_ion_data),
+ [&ion_config, &seg_ion_map](size_type cv_index, index_type cv, value_type area, index_type seg_index, index_type seg) {
+ if (seg_ion_map[seg].mech_requires) {
+ if (!util::binary_search_index(ion_config.cv, cv)) {
+ ion_config.cv.push_back(cv);
+ }
+ }
+ });
- for (size_type segment: ionseg.second) {
- const segment_info& seg_info = D.segments[segment];
+ ion_config.init_iconc.resize(ion_config.cv.size());
+ ion_config.init_econc.resize(ion_config.cv.size());
+ ion_config.init_revpot.resize(ion_config.cv.size());
- if (seg_info.has_parent()) {
- index_type cv = seg_info.parent_cv;
- optional<std::size_t> parent_idx = util::binary_search_index(ion.cv, cv);
- if (!parent_idx) {
- ion.cv.push_back(cv);
- ion.iconc_norm_area.push_back(D.cv_area[cv]);
- ion.econc_norm_area.push_back(D.cv_area[cv]);
+ for_each_cv_in_segments(D, keys(seg_ion_data),
+ [&ion_config, &seg_ion_map, &D](size_type cv_index, index_type cv, value_type area, index_type seg_index, index_type seg) {
+ auto opt_i = util::binary_search_index(ion_config.cv, cv);
+ if (!opt_i) return;
+
+ std::size_t i = *opt_i;
+ auto& seg_ion_entry = seg_ion_map[seg];
+
+ value_type weight = area/D.cv_area[cv];
+ if (!seg_ion_entry.mech_writes_iconc) {
+ ion_config.init_iconc[i] += weight*seg_ion_entry.ion_data.init_int_concentration;
+ }
+
+ if (!seg_ion_entry.mech_writes_econc) {
+ ion_config.init_econc[i] += weight*seg_ion_entry.ion_data.init_ext_concentration;
+ }
+
+ // Reversal potentials are not area weighted, and are overridden at the
+ // per-cell level by any supplied revpot mechanisms.
+ ion_config.init_revpot[i] = seg_ion_entry.ion_data.init_reversal_potential;
+ });
+ }
+
+ // IIb. Reversal potential mechanism CVs and parameters.
+
+ for (const auto& entry: revpot_mech_table) {
+ const std::string& name = entry.first;
+ const mechanism_info& info = *entry.second.info;
+
+ fvm_mechanism_config& config = mechdata.mechanisms[name];
+ config.kind = mechanismKind::revpot;
+
+ auto nparam = build_param_data(entry.second.paramset, &info);
+ config.param_values.resize(nparam);
+ for (auto pidx: make_span(nparam)) {
+ config.param_values[pidx].first = param_name[pidx];
+ }
+
+ for (auto& cell_entry: entry.second.cells) {
+ auto cell_idx = cell_entry.first;
+ std::vector<index_type> segment_indices;
+
+ for (auto& mech_ion_dep_entry: info.ions) {
+ auto& ion_name = mech_ion_dep_entry.first;
+ auto& ion_dep = mech_ion_dep_entry.second;
+
+ if (!ion_dep.write_reversal_potential) continue;
+
+ const auto& segments = value_by_key(ion_revpot_segments[ion_name], cell_idx);
+ if (!segments) continue;
+
+ for (auto seg_idx: segments.value()) {
+ segment_indices.push_back(seg_idx);
}
}
- for (auto cv: make_span(seg_info.cv_range())) {
- ion.cv.push_back(cv);
- ion.iconc_norm_area.push_back(D.cv_area[cv]);
- ion.econc_norm_area.push_back(D.cv_area[cv]);
+ const mechanism_desc& desc = *cell_entry.second;
+ std::vector<value_type> pval = param_default;
+
+ for (auto pidx: make_span(nparam)) {
+ if (auto opt_v = value_by_key(desc.values(), param_name[pidx])) {
+ pval[pidx] = opt_v.value();
+ }
}
+
+ size_type config_offset = config.cv.size();
+ for_each_cv_in_segments(D, segment_indices,
+ [&](size_type cv_index, index_type cv, auto, auto, auto) {
+ extend_at(config.cv, config_offset+cv_index) = cv;
+
+ for (auto pidx: make_span(nparam)) {
+ extend_at(config.param_values[pidx].second, config_offset+cv_index) = pval[pidx];
+ }
+ });
}
}
- // IIb. Density mechanism CVs, parameters and ionic default concentration contributions.
+ // Remove any reversal potential mechanisms that ultimately have no extent.
+ for (auto i = mechdata.mechanisms.begin(); i!=mechdata.mechanisms.end(); ) {
+ i = i->second.cv.empty()? mechdata.mechanisms.erase(i): std::next(i);
+ }
+
+ // IIc. Density mechanism CVs, parameters and ionic default concentration contributions.
- // Ameliorate area sum rounding areas by clamping normalized area contributions to [0, 1]
+ // Ameliorate area sum rounding errors by clamping normalized area contributions to [0, 1]
// and rounding values within an epsilon of 0 or 1 to that value.
auto trim = [](value_type& v) {
constexpr value_type eps = std::numeric_limits<value_type>::epsilon()*4;
@@ -503,7 +852,7 @@ fvm_mechanism_data fvm_build_mechanism_data(const cable_cell_global_properties&
auto nparam = build_param_data(entry.second.paramset, entry.second.info.get());
- // In order to properly account for partially overriden paramaters in CVs
+ // In order to properly account for partially overriden parameters in CVs
// that are shared between segments, we need to track not only the area-weighted
// sum of parameter values, but also the total area for each CV for each parameter
// that has been overriden — the remaining area demands a contribution from the
@@ -512,65 +861,23 @@ fvm_mechanism_data fvm_build_mechanism_data(const cable_cell_global_properties&
std::vector<std::vector<value_type>> param_value(nparam);
std::vector<std::vector<value_type>> param_area_contrib(nparam);
- // (gcc 6.x bug fails to deduce const in lambda capture reference initialization)
- const auto& info = *entry.second.info;
- auto accumulate_mech_data =
- [
- &info,
- &ion_configs = mechdata.ions,
- ¶m_index, ¶m_value, ¶m_area_contrib, &config
- ]
- (size_type index, index_type cv, value_type area, const mechanism_desc& desc)
- {
- for (auto& kv: desc.values()) {
- int pidx = param_index.at(kv.first);
- value_type v = kv.second;
+ for_each_cv_in_segments(D, keys(entry.second.segments),
+ [&](size_type cv_index, index_type cv, value_type area, index_type seg_index, index_type seg)
+ {
+ const mechanism_desc& mech_desc = *entry.second.segments[seg_index].second;
- extend_to(param_area_contrib[pidx], index);
- param_area_contrib[pidx][index] += area;
+ extend_at(config.cv, cv_index) = cv;
- extend_to(param_value[pidx], index);
- param_value[pidx][index] += area*v;
- }
+ for (auto& kv: mech_desc.values()) {
+ int pidx = param_index.at(kv.first);
+ value_type v = kv.second;
- for (auto& ion: info.ions) {
- fvm_ion_config& ion_config = ion_configs[ion.first];
- size_type index = util::binary_search_index(ion_config.cv, cv).value();
- if (ion.second.write_concentration_int) {
- ion_config.iconc_norm_area[index] -= area;
- }
- if (ion.second.write_concentration_ext) {
- ion_config.econc_norm_area[index] -= area;
+ extend_at(param_area_contrib[pidx], cv_index) += area;
+ extend_at(param_value[pidx], cv_index) += area*v;
}
- }
- extend_to(config.norm_area, index);
- config.norm_area[index] += area;
- };
-
- for (auto& seg_entry: entry.second.segments) {
- const segment_info& seg_info = D.segments[seg_entry.first];
- const mechanism_desc& mech_desc = *seg_entry.second;
-
- if (seg_info.has_parent()) {
- index_type cv = seg_info.parent_cv;
- optional<std::size_t> parent_idx = util::binary_search_index(config.cv, cv);
- if (!parent_idx) {
- parent_idx = config.cv.size();
- config.cv.push_back(cv);
- }
-
- accumulate_mech_data(*parent_idx, cv, seg_info.parent_cv_area, mech_desc);
- }
-
- for (auto cv: make_span(seg_info.cv_range())) {
- size_type idx = config.cv.size();
- config.cv.push_back(cv);
-
- value_type area = cv==seg_info.distal_cv? seg_info.distal_cv_area: D.cv_area[cv];
- accumulate_mech_data(idx, cv, area, mech_desc);
- }
- }
+ extend_at(config.norm_area, cv_index) += area;
+ });
// Complete parameter values with default values.
@@ -599,20 +906,6 @@ fvm_mechanism_data fvm_build_mechanism_data(const cable_cell_global_properties&
}
}
- // Normalize ion norm_area entries.
-
- for (auto& entry: mechdata.ions) {
- auto& ion_config = entry.second;
- for (auto i: count_along(ion_config.cv)) {
- auto cv_area = D.cv_area[ion_config.cv[i]];
- ion_config.iconc_norm_area[i] /= cv_area;
- trim(ion_config.iconc_norm_area[i]);
-
- ion_config.econc_norm_area[i] /= cv_area;
- trim(ion_config.econc_norm_area[i]);
- }
- }
-
// II.3 Point mechanism CVs, targets, parameters and stimuli.
for (const auto& entry: point_mech_table) {
@@ -638,7 +931,7 @@ fvm_mechanism_data fvm_build_mechanism_data(const cable_cell_global_properties&
// Generate config.param_values: contains parameters of group of synapses that can be coalesced into one instance
// Generate multiplicity: contains number of synapses in each coalesced group of synapses
// Generate target: contains the synapse target number
- if(catalogue[name].linear && gprop.coalesce_synapses) {
+ if (catalogue[name].linear && gprop.coalesce_synapses) {
// cv_param_vec used to lexicographically sort the cv, parameters and target, which are stored in that order
std::vector<cv_param> cv_param_vec(cv_order.size());
diff --git a/arbor/fvm_layout.hpp b/arbor/fvm_layout.hpp
index a87814ef9396dcdc901edc9f6503cde05f895495..c0172e06e21ba10a15cba1c117a6bcc1d6298b70 100644
--- a/arbor/fvm_layout.hpp
+++ b/arbor/fvm_layout.hpp
@@ -66,8 +66,14 @@ struct fvm_discretization {
std::vector<value_type> face_conductance; // [µS]
std::vector<value_type> cv_area; // [µm²]
std::vector<value_type> cv_capacitance; // [pF]
+ std::vector<value_type> init_membrane_potential; // [mV]
+ std::vector<value_type> temperature_K; // [K]
std::vector<segment_info> segments;
+
+ // If segment has no parent segment, parent_segment[j] = j.
+ std::vector<index_type> parent_segment;
+
std::vector<size_type> cell_segment_bounds; // Partitions segment indices by cell.
std::vector<index_type> cell_cv_bounds; // Partitions CV indices by cell.
@@ -88,7 +94,7 @@ struct fvm_discretization {
}
};
-fvm_discretization fvm_discretize(const std::vector<cable_cell>& cells);
+fvm_discretization fvm_discretize(const std::vector<cable_cell>& cells, const cable_cell_parameter_set& params);
// Post-discretization data for point and density mechanism instantiation.
@@ -103,7 +109,7 @@ struct fvm_mechanism_config {
// duplicates for point mechanisms.
std::vector<index_type> cv;
- // Coalesced synapse multiplier
+ // Coalesced synapse multiplier (point mechanisms only).
std::vector<index_type> multiplicity;
// Normalized area contribution in corresponding CV (density mechanisms only).
@@ -126,8 +132,11 @@ struct fvm_ion_config {
std::vector<index_type> cv;
// Normalized area contribution of default concentration contribution in corresponding CV.
- std::vector<value_type> iconc_norm_area;
- std::vector<value_type> econc_norm_area;
+ std::vector<value_type> init_iconc;
+ std::vector<value_type> init_econc;
+
+ // Ion-specific (initial) reversal potential per CV.
+ std::vector<value_type> init_revpot;
};
struct fvm_mechanism_data {
diff --git a/arbor/fvm_lowered_cell_impl.hpp b/arbor/fvm_lowered_cell_impl.hpp
index b1b13406ca835ca02b6fe67a9964927d102979e3..a466077b19611dcc2132e511a8c49c1ab7f30290 100644
--- a/arbor/fvm_lowered_cell_impl.hpp
+++ b/arbor/fvm_lowered_cell_impl.hpp
@@ -18,6 +18,7 @@
#include <arbor/assert.hpp>
#include <arbor/common_types.hpp>
+#include <arbor/cable_cell_param.hpp>
#include <arbor/recipe.hpp>
#include "builtin_mechanisms.hpp"
@@ -100,9 +101,8 @@ private:
threshold_watcher threshold_watcher_;
value_type tmin_ = 0;
- value_type initial_voltage_ = NAN;
- value_type temperature_ = NAN;
- std::vector<mechanism_ptr> mechanisms_;
+ std::vector<mechanism_ptr> mechanisms_; // excludes reversal potential calculators.
+ std::vector<mechanism_ptr> revpot_mechanisms_;
// Non-physical voltage check threshold, 0 => no check.
value_type check_voltage_mV = 0;
@@ -151,9 +151,13 @@ void fvm_lowered_cell_impl<Backend>::assert_tmin() {
template <typename Backend>
void fvm_lowered_cell_impl<Backend>::reset() {
- state_->reset(initial_voltage_, temperature_);
+ state_->reset();
set_tmin(0);
+ for (auto& m: revpot_mechanisms_) {
+ m->initialize();
+ }
+
for (auto& m: mechanisms_) {
m->initialize();
}
@@ -198,6 +202,15 @@ fvm_integration_result fvm_lowered_cell_impl<Backend>::integrate(
// complete fvm state into shared state object.
while (remaining_steps) {
+ // Update any required reversal potentials based on ionic concs.
+
+ PE(advance_update_revpot)
+ for (auto& m: revpot_mechanisms_) {
+ m->nrn_current();
+ }
+ PL();
+
+
// Deliver events and accumulate mechanism current contributions.
PE(advance_integrate_events);
@@ -299,7 +312,6 @@ void fvm_lowered_cell_impl<B>::update_ion_state() {
for (auto& m: mechanisms_) {
m->write_ions();
}
- state_->ions_nernst_reversal_potential(temperature_);
}
template <typename B>
@@ -355,9 +367,11 @@ void fvm_lowered_cell_impl<B>::initialize(
throw bad_global_property(cell_kind::cable);
}
+ // Assert that all global default parameters have been set.
+ // (Throws cable_cell_error on failure.)
+ check_global_properties(global_props);
+
const mechanism_catalogue* catalogue = global_props.catalogue;
- initial_voltage_ = global_props.init_membrane_potential_mV;
- temperature_ = global_props.temperature_K;
// Mechanism instantiator helper.
auto mech_instance = [&catalogue](const std::string& name) {
@@ -373,7 +387,7 @@ void fvm_lowered_cell_impl<B>::initialize(
// Discretize cells, build matrix.
- fvm_discretization D = fvm_discretize(cells);
+ fvm_discretization D = fvm_discretize(cells, global_props.default_parameters);
std::vector<index_type> cv_to_intdom(D.ncv);
std::transform(D.cv_to_cell.begin(), D.cv_to_cell.end(), cv_to_intdom.begin(),
@@ -387,7 +401,7 @@ void fvm_lowered_cell_impl<B>::initialize(
fvm_mechanism_data mech_data = fvm_build_mechanism_data(global_props, cells, D);
- // Discritize and build gap junction info
+ // Discretize and build gap junction info.
auto gj_vector = fvm_gap_junctions(cells, gids, rec, D);
@@ -398,15 +412,17 @@ void fvm_lowered_cell_impl<B>::initialize(
util::transform_view(keys(mech_data.mechanisms),
[&](const std::string& name) { return mech_instance(name).mech->data_alignment(); }));
- state_ = std::make_unique<shared_state>(num_intdoms, cv_to_intdom, gj_vector, data_alignment? data_alignment: 1u);
+ state_ = std::make_unique<shared_state>(
+ num_intdoms, cv_to_intdom, gj_vector, D.init_membrane_potential, D.temperature_K,
+ data_alignment? data_alignment: 1u);
// Instantiate mechanisms and ions.
for (auto& i: mech_data.ions) {
const std::string& ion_name = i.first;
- if (auto ion = value_by_key(global_props.ion_default, ion_name)) {
- state_->add_ion(ion_name, ion.value(), i.second.cv, i.second.iconc_norm_area, i.second.econc_norm_area);
+ if (auto charge = value_by_key(global_props.ion_species, ion_name)) {
+ state_->add_ion(ion_name, *charge, i.second.cv, i.second.init_iconc, i.second.init_econc, i.second.init_revpot);
}
else {
throw cable_cell_error("unrecognized ion '"+ion_name+"' in mechanism");
@@ -415,10 +431,10 @@ void fvm_lowered_cell_impl<B>::initialize(
target_handles.resize(mech_data.ntarget);
+ unsigned mech_id = 0;
for (auto& m: mech_data.mechanisms) {
auto& name = m.first;
auto& config = m.second;
- unsigned mech_id = mechanisms_.size();
mechanism_layout layout;
layout.cv = config.cv;
@@ -432,7 +448,8 @@ void fvm_lowered_cell_impl<B>::initialize(
// contribution in the CV, and F is the scaling factor required
// to convert from the mechanism current contribution units to A/m².
- if (config.kind==mechanismKind::point) {
+ switch (config.kind) {
+ case mechanismKind::point:
// Point mechanism contributions are in [nA]; CV area A in [µm^2].
// F = 1/A * [nA/µm²] / [A/m²] = 1000/A.
@@ -452,22 +469,33 @@ void fvm_lowered_cell_impl<B>::initialize(
};
}
}
- }
- else {
+ break;
+ case mechanismKind::density:
// Current density contributions from mechanism are already in [A/m²].
for (auto i: count_along(layout.cv)) {
layout.weight[i] = config.norm_area[i];
}
+ break;
+ case mechanismKind::revpot:
+ // Mechanisms that set reversal potential should not be contributing
+ // to any currents, so leave weights as zero.
+ break;
}
auto minst = mech_instance(name);
- minst.mech->instantiate(mech_id, *state_, minst.overrides, layout);
+ minst.mech->instantiate(mech_id++, *state_, minst.overrides, layout);
for (auto& pv: config.param_values) {
minst.mech->set_parameter(pv.first, pv.second);
}
- mechanisms_.push_back(mechanism_ptr(minst.mech.release()));
+
+ if (config.kind==mechanismKind::revpot) {
+ revpot_mechanisms_.push_back(mechanism_ptr(minst.mech.release()));
+ }
+ else {
+ mechanisms_.push_back(mechanism_ptr(minst.mech.release()));
+ }
}
// Collect detectors, probe handles.
diff --git a/arbor/include/arbor/cable_cell.hpp b/arbor/include/arbor/cable_cell.hpp
index b91a6e11f2434c10535c620a72620c6d375dcbee..f790627a60bf5a515fe7c12dd9a5847c22e1775a 100644
--- a/arbor/include/arbor/cable_cell.hpp
+++ b/arbor/include/arbor/cable_cell.hpp
@@ -5,23 +5,15 @@
#include <vector>
#include <arbor/arbexcept.hpp>
+#include <arbor/cable_cell_param.hpp>
#include <arbor/common_types.hpp>
#include <arbor/constants.hpp>
-#include <arbor/ion_info.hpp>
#include <arbor/mechcat.hpp>
#include <arbor/morphology.hpp>
#include <arbor/segment.hpp>
namespace arb {
-// Specialize arbor exception for errors in cell building.
-
-struct cable_cell_error: arbor_exception {
- cable_cell_error(const std::string& what):
- arbor_exception("cable_cell: "+what)
- {}
-};
-
// Location specification for point processes.
struct segment_location {
@@ -64,34 +56,6 @@ struct cell_probe_address {
probe_kind kind;
};
-// Global parameter type for cell descriptions.
-
-struct cable_cell_global_properties {
- const mechanism_catalogue* catalogue = &global_default_catalogue();
-
- // If >0, check membrane voltage magnitude is less than limit
- // during integration.
- double membrane_voltage_limit_mV = 0;
-
- // TODO: consider making some/all of the following parameters
- // cell or even segment-local.
- //
- // Consider also a model-level dictionary of default values that
- // can be used to initialize per-cell-kind info?
- //
- // Defaults below chosen to match NEURON.
-
- std::unordered_map<std::string, ion_info> ion_default = {
- {"ca", { 2, 5e-5, 2. }},
- {"na", { 1, 10., 140.}},
- {"k", { 1, 54.4, 2.5 }}
- };
-
- double temperature_K = constant::hh_squid_temp; // [K]
- double init_membrane_potential_mV = -65; // [mV]
- bool coalesce_synapses = true;
-};
-
/// high-level abstract representation of a cell and its segments
class cable_cell {
public:
@@ -117,11 +81,14 @@ public:
double threshold;
};
+ cable_cell_parameter_set default_parameters;
+
/// Default constructor
cable_cell();
/// Copy constructor
cable_cell(const cable_cell& other):
+ default_parameters(other.default_parameters),
parents_(other.parents_),
stimuli_(other.stimuli_),
synapses_(other.synapses_),
@@ -176,6 +143,7 @@ public:
/// will throw an cable_cell_error exception if
/// the cable index is not valid
cable_segment* cable(index_type index);
+ const cable_segment* cable(index_type index) const;
/// the total number of compartments over all segments
size_type num_compartments() const;
@@ -251,6 +219,20 @@ public:
return parents_;
}
+ // Approximate per-segment mean attenuation b(f) at given frequency f,
+ // ignoring membrane resistance [1/µm].
+ value_type segment_mean_attenuation(value_type frequency, index_type segidx,
+ const cable_cell_parameter_set& global_defaults) const;
+
+ // Estimate of length constant λ(f) = 1/2 · 1/b(f), following
+ // Hines and Carnevale (2001), "NEURON: A Tool for Neuroscientists",
+ // Neuroscientist 7, pp. 123-135.
+ value_type segment_length_constant(value_type frequency, index_type segidx,
+ const cable_cell_parameter_set& global_defaults) const
+ {
+ return 0.5/segment_mean_attenuation(frequency, segidx, global_defaults);
+ }
+
private:
void assert_valid_segment(index_type) const;
diff --git a/arbor/include/arbor/cable_cell_param.hpp b/arbor/include/arbor/cable_cell_param.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..2c147c0ce25667eb91ee1b0ee12e285e4cbfc754
--- /dev/null
+++ b/arbor/include/arbor/cable_cell_param.hpp
@@ -0,0 +1,157 @@
+#pragma once
+
+#include <arbor/arbexcept.hpp>
+#include <arbor/mechcat.hpp>
+#include <arbor/util/optional.hpp>
+
+#include <unordered_map>
+#include <string>
+
+namespace arb {
+
+// Specialized arbor exception for errors in cell building.
+
+struct cable_cell_error: arbor_exception {
+ cable_cell_error(const std::string& what):
+ arbor_exception("cable_cell: "+what) {}
+};
+
+// Mechanism description, viz. mechanism name and
+// (non-global) parameter settings. Used to assign
+// density and point mechanisms to segments and
+// reversal potential computations to cells.
+
+struct mechanism_desc {
+ struct field_proxy {
+ mechanism_desc* m;
+ std::string key;
+
+ field_proxy& operator=(double v) {
+ m->set(key, v);
+ return *this;
+ }
+
+ operator double() const {
+ return m->get(key);
+ }
+ };
+
+ // implicit
+ mechanism_desc(std::string name): name_(std::move(name)) {}
+ mechanism_desc(const char* name): name_(name) {}
+
+ mechanism_desc() = default;
+ mechanism_desc(const mechanism_desc&) = default;
+ mechanism_desc(mechanism_desc&&) = default;
+
+ mechanism_desc& operator=(const mechanism_desc&) = default;
+ mechanism_desc& operator=(mechanism_desc&&) = default;
+
+ mechanism_desc& set(const std::string& key, double value) {
+ param_[key] = value;
+ return *this;
+ }
+
+ double operator[](const std::string& key) const {
+ return get(key);
+ }
+
+ field_proxy operator[](const std::string& key) {
+ return {this, key};
+ }
+
+ double get(const std::string& key) const {
+ auto i = param_.find(key);
+ if (i==param_.end()) {
+ throw std::out_of_range("no field "+key+" set");
+ }
+ return i->second;
+ }
+
+ const std::unordered_map<std::string, double>& values() const {
+ return param_;
+ }
+
+ const std::string& name() const { return name_; }
+
+private:
+ std::string name_;
+ std::unordered_map<std::string, double> param_;
+};
+
+// Cable cell ion and electrical defaults.
+//
+// Parameters can be overridden with `cable_cell_local_parameter_set`
+// on unbranched segments within a cell; per-cell and global defaults
+// use `cable_cell_parameter_set`, which extends the parameter set
+// to supply per-cell or global ion reversal potential calculation
+// mechanisms.
+
+struct cable_cell_ion_data {
+ double init_int_concentration = NAN;
+ double init_ext_concentration = NAN;
+ double init_reversal_potential = NAN;
+};
+
+struct cable_cell_local_parameter_set {
+ std::unordered_map<std::string, cable_cell_ion_data> ion_data;
+ util::optional<double> init_membrane_potential; // [mV]
+ util::optional<double> temperature_K; // [K]
+ util::optional<double> axial_resistivity; // [Ω·cm]
+ util::optional<double> membrane_capacitance; // [F/m²]
+};
+
+struct cable_cell_parameter_set: public cable_cell_local_parameter_set {
+ std::unordered_map<std::string, mechanism_desc> reversal_potential_method;
+
+ // We'll need something like this until C++17, for sane initialization syntax.
+ cable_cell_parameter_set() = default;
+ cable_cell_parameter_set(cable_cell_local_parameter_set p, std::unordered_map<std::string, mechanism_desc> m = {}):
+ cable_cell_local_parameter_set(std::move(p)), reversal_potential_method(std::move(m))
+ {}
+};
+
+extern cable_cell_local_parameter_set neuron_parameter_defaults;
+
+// Global cable cell data.
+
+struct cable_cell_global_properties {
+ const mechanism_catalogue* catalogue = &global_default_catalogue();
+
+ // If >0, check membrane voltage magnitude is less than limit
+ // during integration.
+ double membrane_voltage_limit_mV = 0;
+
+ // True => combine linear synapses for performance.
+ bool coalesce_synapses = true;
+
+ // Available ion species, together with charge.
+ std::unordered_map<std::string, int> ion_species = {
+ {"na", 1},
+ {"k", 1},
+ {"ca", 2}
+ };
+
+ cable_cell_parameter_set default_parameters;
+
+ // Convenience methods for adding a new ion together with default ion values.
+ void add_ion(const std::string& ion_name, int charge, double init_iconc, double init_econc, double init_revpot) {
+ ion_species[ion_name] = charge;
+
+ auto &ion_data = default_parameters.ion_data[ion_name];
+ ion_data.init_int_concentration = init_iconc;
+ ion_data.init_ext_concentration = init_econc;
+ ion_data.init_reversal_potential = init_revpot;
+ }
+
+ void add_ion(const std::string& ion_name, int charge, double init_iconc, double init_econc, mechanism_desc revpot_mechanism) {
+ add_ion(ion_name, charge, init_iconc, init_econc, 0);
+ default_parameters.reversal_potential_method[ion_name] = std::move(revpot_mechanism);
+ }
+};
+
+// Throw cable_cell_error if any default parameters are left unspecified,
+// or if the supplied ion data is incomplete.
+void check_global_properties(const cable_cell_global_properties&);
+
+} // namespace arb
diff --git a/arbor/include/arbor/ion_info.hpp b/arbor/include/arbor/ion_info.hpp
deleted file mode 100644
index a3e97263b091bc592cfdc7decc375e6fd4a24cf5..0000000000000000000000000000000000000000
--- a/arbor/include/arbor/ion_info.hpp
+++ /dev/null
@@ -1,11 +0,0 @@
-#pragma once
-
-namespace arb {
-
-struct ion_info {
- int charge; // charge of ionic species
- double default_int_concentration; // (mM) default internal concentration
- double default_ext_concentration; // (mM) default external concentration
-};
-
-} // namespace arb
diff --git a/arbor/include/arbor/mechanism.hpp b/arbor/include/arbor/mechanism.hpp
index b6ed27f974ea9eff7400549abd211aa58f0dcf22..37e497f865e7965c8591c7d94f6f188e697f9ffb 100644
--- a/arbor/include/arbor/mechanism.hpp
+++ b/arbor/include/arbor/mechanism.hpp
@@ -9,7 +9,7 @@
namespace arb {
-enum class mechanismKind {point, density};
+enum class mechanismKind {point, density, revpot};
class mechanism;
using mechanism_ptr = std::unique_ptr<mechanism>;
diff --git a/arbor/include/arbor/segment.hpp b/arbor/include/arbor/segment.hpp
index 1f408850f1b0757a6e30fa49d22c8343cf0f48b2..38f1849c65e0b8a8c81ba722ca5d7744f0446c38 100644
--- a/arbor/include/arbor/segment.hpp
+++ b/arbor/include/arbor/segment.hpp
@@ -10,6 +10,7 @@
#include <vector>
#include <arbor/assert.hpp>
+#include <arbor/cable_cell_param.hpp>
#include <arbor/common_types.hpp>
#include <arbor/math.hpp>
#include <arbor/morphology.hpp>
@@ -19,59 +20,6 @@
namespace arb {
-// Mechanism information attached to a segment.
-
-struct mechanism_desc {
- struct field_proxy {
- mechanism_desc* m;
- std::string key;
-
- field_proxy& operator=(double v) {
- m->set(key, v);
- return *this;
- }
-
- operator double() const {
- return m->get(key);
- }
- };
-
- // implicit
- mechanism_desc(std::string name): name_(std::move(name)) {}
- mechanism_desc(const char* name): name_(name) {}
-
- mechanism_desc& set(const std::string& key, double value) {
- param_[key] = value;
- return *this;
- }
-
- double operator[](const std::string& key) const {
- return get(key);
- }
-
- field_proxy operator[](const std::string& key) {
- return {this, key};
- }
-
- double get(const std::string& key) const {
- auto i = param_.find(key);
- if (i==param_.end()) {
- throw std::out_of_range("no field "+key+" set");
- }
- return i->second;
- }
-
- const std::unordered_map<std::string, double>& values() const {
- return param_;
- }
-
- const std::string& name() const { return name_; }
-
-private:
- std::string name_;
- std::unordered_map<std::string, double> param_;
-};
-
// forward declarations of segment specializations
class soma_segment;
class cable_segment;
@@ -164,9 +112,7 @@ public:
return mechanisms_;
}
- // common electrical properties
- value_type rL = 100.0; // resistivity [Ohm.cm]
- value_type cm = 0.01; // capacitance [F/m^2] : 10 nF/mm^2 = 0.01 F/m^2
+ cable_cell_local_parameter_set parameters; // override cell and global defaults
protected:
segment(section_kind kind): kind_(kind) {}
@@ -297,25 +243,6 @@ public:
return std::accumulate(lengths_.begin(), lengths_.end(), value_type{});
}
- value_type length_constant(value_type freq_Hz) const override {
- // Following Hine and Carnevale (2001), "NEURON: A Tool for Neuroscientists",
- // Neuroscientist 7, pp. 123-135.
- //
- // λ(f) = approx. sqrt(diameter/(pi*f*rL*cm))/2.
- //
- // Pick smallest non-zero diameter in the segment.
-
- value_type r_min = 0;
- for (auto r: radii_) {
- if (r>0 && (r_min==0 || r<r_min)) r_min = r;
- }
- value_type d_min = r_min*2e-6; // [m]
- value_type rc = rL*0.01*cm; // [s/m]
- value_type lambda = std::sqrt(d_min/(math::pi<double>*freq_Hz*rc))/2.; // [m]
-
- return lambda*1e6; // [µm]
- }
-
bool has_locations() const
{
return locations_.size() > 0;
diff --git a/arbor/util/maputil.hpp b/arbor/util/maputil.hpp
index afb2f3ad254b323bd5f68022a8f5b4ade883c60a..44956499488ef3d4e5828a7583fa958472638f08 100644
--- a/arbor/util/maputil.hpp
+++ b/arbor/util/maputil.hpp
@@ -121,6 +121,19 @@ optional<typename sequence_traits<C>::difference_type> binary_search_index(const
return it!=strict.end() && key==*it? just(std::distance(strict.begin(), it)): nullopt;
}
+// As binary_search_index above, but compare key against the proj(x) for elements
+// x in the sequence. The image of proj applied to the sequence must be monotonically
+// increasing.
+
+template <typename C, typename Key, typename Proj>
+optional<typename sequence_traits<C>::difference_type> binary_search_index(const C& c, const Key& key, const Proj& proj) {
+ auto strict = strict_view(c);
+ auto projected = transform_view(strict, proj);
+ auto it = std::lower_bound(projected.begin(), projected.end(), key);
+ return it!=strict.end() && key==*it? just(std::distance(projected.begin(), it)): nullopt;
+}
+
+
// Key equality helper for NUL-terminated strings.
struct cstr_equal {
diff --git a/doc/cpp_cable_cell.rst b/doc/cpp_cable_cell.rst
new file mode 100644
index 0000000000000000000000000000000000000000..4cf1dce867ae1eeae7a954ce6afde27f96db56a6
--- /dev/null
+++ b/doc/cpp_cable_cell.rst
@@ -0,0 +1,260 @@
+.. _cppcablecell:
+
+Cable cells
+===============
+
+.. Warning::
+ The interface for building and modifying cable cell objects
+ will be thoroughly revised in the near future. The documentation
+ here is primarily a place holder.
+
+Cable cells, which use the :cpp:enum:`cell_kind` :cpp:expr:`cable`,
+represent morphologically-detailed neurons as 1-d trees, with
+electrical and biophysical properties mapped onto those trees.
+
+A single cell is represented by an object of type :cpp:type:`cable_cell`.
+Properties shared by all cable cells, as returned by the recipe
+:cpp:expr:`get_global_properties` method, are described by an
+object of type :cpp:type:`cable_cell_global_properties`.
+
+
+The :cpp:type:`cable_cell` object
+---------------------------------
+
+Cable cells are built up from a series of :cpp:type:`segment`
+objects, which themselves describe an unbranched component of the
+cell morphology. These segments are added via the methods:
+
+.. cpp:function:: soma_segment* cable_cell::add_soma(double radius)
+
+ Add the soma to the cable cell with the given radius. There
+ can be only one per cell.
+
+ The soma segment has index 0, and must be added before any
+ cable segments.
+
+.. cpp:function:: cable_segment* cable_cell::add_cable(cell_lid_type index, Args&&... args)
+
+ Add a unbranched section of the cell morphology, with its proximal
+ end attached to the segment given by :cpp:expr:`index`. The
+ following arguments are forwarded to the :cpp:type:`cable_segment`
+ constructor.
+
+Segment indices are exactly the order in which they have been added
+to a cell, counting from zero (for the soma). Both :cpp:type:`soma_segment`
+and :cpp:type:`cable_segment` are derived from the abstract base
+class :cpp:type:`segment`.
+
+.. todo::
+
+ Describe cable_segment constructor arguments, unless we get to the
+ replace cell building/morphology implementation first.
+
+Each segment will inherit the electrical properties of the cell, unless
+otherwise overriden (see below).
+
+
+Cell dynamics
+-------------
+
+Each segment in a cell may have attached to it one or more density *mechanisms*,
+which describe biophysical processes. These are processes
+that are distributed in space, but whose behaviour is defined purely
+by the state of the cell and the process at any given point.
+
+Cells may also have *point* mechanisms, which are added directly to the
+:cpp:type:`cable_cell` object.
+
+A third type of mechanism, which describes ionic reversal potential
+behaviour, can be specified for cells or the whole model via cell parameter
+settings, described below.
+
+Mechanisms are described by a :cpp:type:`mechanism_desc` object. These specify the
+name of the mechanism (used to find the mechanism in the mechanism catalogue)
+and parameter values for the mechanism that apply within a segment.
+A :cpp:type:`mechanism_desc` is effectively a wrapper around a name and
+a dictionary of parameter/value settings.
+
+Mechanism descriptions can be constructed implicitly from the
+mechanism name, and mechanism parameter values then set with the
+:cpp:expr:`set` method. Relevant :cpp:type:`mechanism_desc` methods:
+
+.. cpp:function:: mechanism_desc::mechanism_desc(std::string name)
+
+ Construct a mechanism description for the mechanism named `name`.
+
+.. cpp:function:: mechanism_desc& mechanism_desc::set(const std::string& key, double value)
+
+ Sets the parameter associated with :cpp:expr:`key` in the description.
+ Returns a reference to the mechanism description, so that calls to
+ :cpp:expr:`set` can be chained in a single expression.
+
+
+Density mechanisms are associated with a cable cell object with:
+
+.. cpp:function:: void segment::add_mechanism(mechanism_desc mech)
+
+Point mechanisms, which are associated with connection end points on a
+cable cell, are attached to a cell with:
+
+.. cpp:function:: void cable_cell::add_synapse(segment_location loc, mechanism_desc mech)
+
+where :cpp:type:`segment_location` is a simple struct holding a segment index
+and a relative position (from 0, proximal, to 1, distal) along that segment:
+
+.. cpp:function:: segment_location::segment_location(cell_lid_type index, double position)
+
+
+Electrical properities and ion values
+-------------------------------------
+
+On each cell segment, electrical and ion properties can be specified by the
+:cpp:expr:`parameters` field, of type :cpp:type:`cable_cell_local_parameter_set`.
+
+The :cpp:type:`cable_cell_local_parameter_set` has the following members,
+where an empty optional value or missing map key indicates that the corresponding
+value should be taken from the cell or global parameter set.
+
+.. cpp:class:: cable_cell_local_parameter_set
+
+ .. cpp:member:: std::unordered_map<std::string, cable_cell_ion_data> ion_data
+
+ The keys of this map are names of ions, whose parameters will be locally overriden.
+ The struct :cpp:type:`cable_cell_ion_data` has three fields:
+ :cpp:type:`init_int_concentration`, :cpp:type:`init_ext_concentration`, and
+ :cpp:type:`init_reversal_potential`.
+
+ Internal and external concentrations are given in millimolars, i.e. mol/m³.
+ Reversal potential is given in millivolts.
+
+ .. cpp:member:: util::optional<double> init_membrane_potential
+
+ Initial membrane potential in millivolts.
+
+ .. cpp:member:: util::optional<double> temperature_K
+
+ Local temperature in Kelvin.
+
+ .. cpp:member:: util::optional<double> axial_resistivity
+
+ Local resistivity of the intracellular medium, in ohm-centimetres.
+
+ .. cpp:member:: util::optional<double> membrane_capacitance
+
+ Local areal capacitance of the cell membrane, in Farads per square metre.
+
+
+Default parameters for a cell are given by the :cpp:expr:`default_parameters`
+field in the :cpp:type:`cable_cell` object. This is a value of type :cpp:type:`cable_cell_parameter_set`,
+which extends :cpp:type:`cable_cell_local_parameter_set` by adding an additional
+field describing reversal potential computation:
+
+.. cpp:class:: cable_cell_parameter_set: public cable_cell_local_parameter_set
+
+ .. cpp:member:: std::unordered_map<std::string, mechanism_desc> reversal_potential_method
+
+ Maps the name of an ion to a 'reversal potential' mechanism that describes
+ how it should be computed. When no mechanism is provided for an ionic reversal
+ potential, the reversal potential will be kept at its initial value.
+
+Default parameters for all cells are supplied in the :cpp:type:`cable_cell_global_properties`
+struct.
+
+Global properties
+-----------------
+
+.. cpp:class:: cable_cell_global_properties
+
+ .. cpp:member:: const mechanism_catalogue* catalogue
+
+ All mechanism names refer to mechanism instances in this mechanism catalogue.
+ By default, this is set to point to `global_default_catalogue()`, the catalogue
+ that contains all mechanisms bundled with Arbor.
+
+ .. cpp:member:: double membrane_voltage_limit_mV
+
+ If non-zero, check to see if the membrane voltage ever exceeds this value
+ in magnitude during the course of a simulation. If so, throw an exception
+ and abort the simulation.
+
+ .. cpp:member:: bool coalesce_synapses
+
+ When synapse dynamics are sufficiently simple, the states of synapses within
+ the same discretized element can be combined for better performance. This
+ is true by default.
+
+ .. cpp:member:: std::unordered_map<std::string, int> ion_species
+
+ Every ion species used by cable cells in the simulation must have an entry in
+ this map, which takes an ion name to its charge, expressed as a multiple of
+ the elementary charge. By default, it is set to include sodium "na" with
+ charge 1, calcium "ca" with charge 2, and potassium "k" with charge 1.
+
+ .. cpp:member:: cable_cell_parameter_set default_parameters
+
+ The default electrical and physical properties associated with each cable
+ cell, unless overridden locally. In the global properties, *every
+ optional field must be given a value*, and every ion must have its default
+ values set in :cpp:expr:`default_parameters.ion_data`.
+
+ .. cpp:function:: add_ion(const std::string& ion_name, int charge, double init_iconc, double init_econc, double init_revpot)
+
+ Convenience function for adding a new ion to the global :cpp:expr:`ion_species`
+ table, and setting up its default values in the `ion_data` table.
+
+ .. cpp:function:: add_ion(const std::string& ion_name, int charge, double init_iconc, double init_econc, mechanism_desc revpot_mechanism)
+
+ As above, but set the initial reversal potential to zero, and use the given mechanism
+ for reversal potential calculation.
+
+
+For convenience, :cpp:expr:`neuron_parameter_defaults` is a predefined :cpp:type:`cable_cell_local_parameter_set`
+value that holds values that correspond to NEURON defaults. To use these values,
+assign them to the :cpp:expr:`default_parameters` field of the global properties
+object returned in the recipe.
+
+
+Reversal potential dynamics
+---------------------------
+
+If no reversal potential mechanism is specified for an ion species, the initial
+reversal potential values are maintained for the course of a simulation. Otherwise,
+a provided mechanism does the work, but it subject to some strict restrictions.
+A reversal potential mechanism described in NMODL:
+
+* May not maintain any STATE variables.
+* Can only write to the "eX" value associated with an ion.
+* Can not given as a POINT mechanism.
+
+Essentially, reversal potential mechanisms must be pure functions of cellular
+and ionic state.
+
+If a reversal potential mechanism writes to multiple ions, then if the mechanism
+is given for one of the ions in the global or per-cell parameters, it must be
+given for all of them.
+
+Arbor's default catalogue includes a "nernst" reversal potential, which is
+parameterized over a single ion, and so can be assigned to e.g. calcium in
+the global parameters via
+
+.. code::
+
+ cable_cell_global_properties gprop;
+ // ...
+ gprop.default_parameters.reversal_potential_method["ca"] = "nernst/ca";
+
+
+This mechanism has global scalar parameters for the gas constant *R* and
+Faraday constant *F*, corresponding to the exact values given by the 2019
+redifinition of the SI base units. These values can be changed in a derived
+mechanism in order to use, for example, older values of these physical
+constants.
+
+.. code::
+
+ mechanism_catalogue mycat(global_default_catalogue());
+ mycat.derive("nernst1998", "nernst", {{"R", 8.314472}, {"F", 96485.3415}});
+
+ gprop.catalogue = &mycat;
+ gprop.default_parameters.reversal_potential_method["ca"] = "nernst1998/ca";
+
diff --git a/doc/index.rst b/doc/index.rst
index d027f7c7d082acf76d2569cfe79d34d0834a8cb9..b87dc5b1aa19b83fb6e97ee72d0396125c88de82 100644
--- a/doc/index.rst
+++ b/doc/index.rst
@@ -104,6 +104,7 @@ Alternative citation formats for the paper can be `downloaded here <https://ieee
cpp_recipe
cpp_domdec
cpp_simulation
+ cpp_cable_cell
.. toctree::
:caption: Developers:
diff --git a/mechanisms/CMakeLists.txt b/mechanisms/CMakeLists.txt
index 3597672e3af91f650dbaaa468a40317c735fe023..8ce3653b6dc45a91e5538cd5f030fb98d6c1091b 100644
--- a/mechanisms/CMakeLists.txt
+++ b/mechanisms/CMakeLists.txt
@@ -1,7 +1,7 @@
include(BuildModules.cmake)
# The list of library-provided mechanisms used to populate the default catalogue:
-set(mechanisms pas hh expsyn exp2syn test_kin1 test_kinlva test_ca nax kdrmt kamt)
+set(mechanisms pas hh expsyn exp2syn nernst test_kin1 test_kinlva test_ca nax kdrmt kamt)
set(mod_srcdir "${CMAKE_CURRENT_SOURCE_DIR}/mod")
diff --git a/mechanisms/mod/nernst.mod b/mechanisms/mod/nernst.mod
new file mode 100644
index 0000000000000000000000000000000000000000..aca2a7170ff23933f143804e72d02491863b5ca1
--- /dev/null
+++ b/mechanisms/mod/nernst.mod
@@ -0,0 +1,33 @@
+: Compute reversal potential for generic ion using Nernst equation.
+:
+: TODO: support global derived parameters via e.g. GLOBAL + ASSIGNED,
+: so that we don't need to store a copy of coeff at every CV.
+
+NEURON {
+ SUFFIX nernst
+ USEION x READ xi, xo WRITE ex VALENCE zx
+ GLOBAL R, F
+ RANGE coeff
+}
+
+PARAMETER {
+ R = 8.31446261815324 (joule/kelvin/mole) : gas constant
+ F = 96485.3321233100184 (coulomb/mole) : Faraday's constant
+ celsius : temperature in °C (set externally)
+}
+
+ASSIGNED {
+ coeff
+}
+
+INITIAL {
+ coeff = R*(celsius+273.15)/(zx*F)*1000
+}
+
+STATE {
+}
+
+BREAKPOINT {
+ ex = coeff*log(xo/xi)
+}
+
diff --git a/modcc/blocks.hpp b/modcc/blocks.hpp
index 5d9d1cf85fcbb52a31966927d8834f4b42b1e8e4..a2c098aa5e81cbc0dd6df3ec5bb64243a5126c27 100644
--- a/modcc/blocks.hpp
+++ b/modcc/blocks.hpp
@@ -34,6 +34,9 @@ struct IonDep {
bool uses_concentration_ext() const {
return has_variable(name+"o");
};
+ bool writes_current() const {
+ return writes_variable("i"+name);
+ };
bool writes_concentration_int() const {
return writes_variable(name+"i");
};
diff --git a/modcc/identifier.hpp b/modcc/identifier.hpp
index 8e62534a999740c71d76a280c90c25cd702ca329..d1c900c978606bc7f45434df4f6cfa6f3c66a011 100644
--- a/modcc/identifier.hpp
+++ b/modcc/identifier.hpp
@@ -5,7 +5,7 @@
#include <stdexcept>
enum class moduleKind {
- point, density
+ point, density, revpot
};
/// indicate how a variable is accessed
@@ -117,7 +117,7 @@ inline sourceKind ion_source(const std::string& ion, const std::string& var, mod
else if (var=="i"+ion) return mkind==moduleKind::point? sourceKind::ion_current: sourceKind::ion_current_density;
else if (var=="e"+ion) return sourceKind::ion_revpot;
else if (var==ion+"i") return sourceKind::ion_iconc;
- else if (var==ion+"e") return sourceKind::ion_econc;
+ else if (var==ion+"o") return sourceKind::ion_econc;
else return sourceKind::no_source;
}
diff --git a/modcc/module.cpp b/modcc/module.cpp
index 0e1ba7674b124c47fe6dd696359bc509735b7613..86a3ebec8777c7e8e706f86e046bb23563269143 100644
--- a/modcc/module.cpp
+++ b/modcc/module.cpp
@@ -447,6 +447,12 @@ bool Module::semantic() {
}
linear_ = linear;
+ // Are we writing an ionic reversal potential? If so, change the moduleKind to
+ // `revpot` and assert that the mechanism is 'pure': it has no state variables;
+ // it contributes to no currents, ionic or otherwise; it isn't a point mechanism;
+ // and it writes to only one ionic reversal potential.
+ check_revpot_mechanism();
+
// Perform semantic analysis and inlining on function and procedure bodies
// in order to inline calls inside the newly crated API methods.
semantic_func_proc();
@@ -765,3 +771,43 @@ int Module::semantic_func_proc() {
return errors;
}
+void Module::check_revpot_mechanism() {
+ int n_write_revpot = 0;
+ for (auto& iondep: neuron_block_.ions) {
+ if (iondep.writes_rev_potential()) ++n_write_revpot;
+ }
+
+ if (n_write_revpot==0) return;
+
+ bool pure = true;
+
+ // Are we marked as a point mechanism?
+ if (kind()==moduleKind::point) {
+ error("Cannot write reversal potential from a point mechanism.");
+ return;
+ }
+
+ // Do we write any other ionic variables or a nonspecific current?
+ for (auto& iondep: neuron_block_.ions) {
+ pure &= !iondep.writes_concentration_int();
+ pure &= !iondep.writes_concentration_ext();
+ pure &= !iondep.writes_current();
+ }
+ pure &= !neuron_block_.has_nonspecific_current();
+
+ if (!pure) {
+ error("Cannot write reversal potential and also to other currents or ionic state.");
+ return;
+ }
+
+ // Do we have any state variables?
+ pure &= state_block_.state_variables.size()==0;
+ if (!pure) {
+ error("Cannot write reversal potential and also maintain state variables.");
+ return;
+ }
+
+ kind_ = moduleKind::revpot;
+}
+
+
diff --git a/modcc/module.hpp b/modcc/module.hpp
index 4cf5e6ac5a1a0b8727a0c09042e5d1c9f8e631ff..09e29d407594c1f3ef2201ca364083f065e185c7 100644
--- a/modcc/module.hpp
+++ b/modcc/module.hpp
@@ -139,6 +139,9 @@ private:
return s == symbols_.end() ? false : s->second->kind() == kind;
}
+ // Check requirements for reversal potential setters.
+ void check_revpot_mechanism();
+
// Perform semantic analysis on functions and procedures.
// Returns the number of errors that were encountered.
int semantic_func_proc();
diff --git a/python/cells.cpp b/python/cells.cpp
index 26f86c22a4650f2d1774deb96fddb95d2a3ac698..ed332548bada6508f605b220fb15fc3bdeb5e1d0 100644
--- a/python/cells.cpp
+++ b/python/cells.cpp
@@ -86,9 +86,10 @@ double interp(const std::array<T,2>& r, unsigned i, unsigned n) {
arb::cable_cell make_cable_cell(arb::cell_gid_type gid, const cell_parameters& params) {
arb::cable_cell cell;
+ cell.default_parameters.axial_resistivity = 100;
+
// Add soma.
auto soma = cell.add_soma(12.6157/2.0); // For area of 500 μm².
- soma->rL = 100;
soma->add_mechanism("hh");
std::vector<std::vector<unsigned>> levels;
@@ -117,7 +118,6 @@ arb::cable_cell make_cable_cell(arb::cell_gid_type gid, const cell_parameters& p
auto dend = cell.add_cable(sec, arb::section_kind::dendrite, dend_radius, dend_radius, l);
dend->set_compartments(nc);
dend->add_mechanism("pas");
- dend->rL = 100;
}
}
}
diff --git a/python/recipe.hpp b/python/recipe.hpp
index f8ccaef3448481c78ff38efc453491ccdb21506b..3ee71adbb2004e4c374436643869cd5f97581416 100644
--- a/python/recipe.hpp
+++ b/python/recipe.hpp
@@ -7,6 +7,7 @@
#include <pybind11/stl.h>
#include <arbor/event_generator.hpp>
+#include <arbor/cable_cell_param.hpp>
#include <arbor/recipe.hpp>
namespace pyarb {
@@ -150,6 +151,11 @@ public:
// TODO: wrap
arb::util::any get_global_properties(arb::cell_kind kind) const override {
+ if (kind==arb::cell_kind::cable) {
+ arb::cable_cell_global_properties gprop;
+ gprop.default_parameters = arb::neuron_parameter_defaults;
+ return gprop;
+ }
return arb::util::any{};
}
};
diff --git a/test/common_cells.hpp b/test/common_cells.hpp
index ae15233e463645caec0ae67c0f74af6919dad4c4..0694bf6daf874a97904bde52815fb77b5dd52eab 100644
--- a/test/common_cells.hpp
+++ b/test/common_cells.hpp
@@ -249,11 +249,12 @@ inline cable_cell make_cell_ball_and_3stick(bool with_stim = true) {
inline cable_cell make_cell_simple_cable(bool with_stim = true) {
cable_cell c;
+ c.default_parameters.axial_resistivity = 100;
+ c.default_parameters.membrane_capacitance = 0.01;
+
c.add_soma(0);
c.add_cable(0, section_kind::dendrite, 0.5, 0.5, 1000);
- double r_L = 100;
- double c_m = 0.01;
double gbar = 0.000025;
double I = 0.1;
@@ -262,8 +263,6 @@ inline cable_cell make_cell_simple_cable(bool with_stim = true) {
for (auto& seg: c.segments()) {
seg->add_mechanism(pas);
- seg->rL = r_L;
- seg->cm = c_m;
if (seg->is_dendrite()) {
seg->set_compartments(4);
diff --git a/test/simple_recipes.hpp b/test/simple_recipes.hpp
index e6da789a09a8b0bf347d4fbb1c1a825ef13c637a..331148055ce581197a6e83bd6bfa575fe50f8aa6 100644
--- a/test/simple_recipes.hpp
+++ b/test/simple_recipes.hpp
@@ -7,6 +7,7 @@
#include <arbor/event_generator.hpp>
#include <arbor/cable_cell.hpp>
+#include <arbor/cable_cell_param.hpp>
#include <arbor/recipe.hpp>
namespace arb {
@@ -20,6 +21,7 @@ public:
catalogue_(global_default_catalogue())
{
cell_gprop_.catalogue = &catalogue_;
+ cell_gprop_.default_parameters = neuron_parameter_defaults;
}
cell_size_type num_probes(cell_gid_type i) const override {
@@ -50,8 +52,12 @@ public:
return catalogue_;
}
- void add_ion(const char* name, int charge, double iconc, double econc) {
- cell_gprop_.ion_default[name] = {charge, iconc, econc};
+ void add_ion(const std::string& ion_name, int charge, double init_iconc, double init_econc, double init_revpot) {
+ cell_gprop_.add_ion(ion_name, charge, init_iconc, init_econc, init_revpot);
+ }
+
+ void nernst_ion(const std::string& ion_name) {
+ cell_gprop_.default_parameters.reversal_potential_method[ion_name] = "nernst/"+ion_name;
}
protected:
diff --git a/test/unit/CMakeLists.txt b/test/unit/CMakeLists.txt
index 36e66d6f508be26d9a1e35c0ed5b9cd40a968398..45d8b30b01036e3690f5f8b69d476d68d9f8783f 100644
--- a/test/unit/CMakeLists.txt
+++ b/test/unit/CMakeLists.txt
@@ -7,6 +7,9 @@ set(test_mechanisms
linear_ca_conc
test_cl_valence
test_ca_read_valence
+ read_eX
+ write_multiple_eX
+ write_eX
)
include(${PROJECT_SOURCE_DIR}/mechanisms/BuildModules.cmake)
diff --git a/test/unit/mod/read_eX.mod b/test/unit/mod/read_eX.mod
new file mode 100644
index 0000000000000000000000000000000000000000..072334bff4f2d58220b282605640d9f5b911a605
--- /dev/null
+++ b/test/unit/mod/read_eX.mod
@@ -0,0 +1,23 @@
+: Test mechanism for consuming ionic reversal potential
+
+NEURON {
+ SUFFIX read_eX
+ USEION x READ ex
+}
+
+PARAMETER {}
+
+ASSIGNED {}
+
+STATE {
+ record_ex
+}
+
+INITIAL {
+ record_ex = ex
+}
+
+BREAKPOINT {
+ record_ex = ex
+}
+
diff --git a/test/unit/mod/write_eX.mod b/test/unit/mod/write_eX.mod
new file mode 100644
index 0000000000000000000000000000000000000000..62302eb980ddb9074fbc679d2f123a002beceb88
--- /dev/null
+++ b/test/unit/mod/write_eX.mod
@@ -0,0 +1,22 @@
+: Reversal potential writer that writes to a single ion
+: as a simple function of internal concentration.
+
+NEURON {
+ SUFFIX write_eX
+ USEION x READ xi WRITE ex
+}
+
+PARAMETER {}
+
+ASSIGNED {}
+
+STATE {}
+
+INITIAL {
+ ex = 100 + xi
+}
+
+BREAKPOINT {
+ ex = 100 + xi
+}
+
diff --git a/test/unit/mod/write_multiple_eX.mod b/test/unit/mod/write_multiple_eX.mod
new file mode 100644
index 0000000000000000000000000000000000000000..a48c571a3910e759386c1e31b0b460c7b7601843
--- /dev/null
+++ b/test/unit/mod/write_multiple_eX.mod
@@ -0,0 +1,25 @@
+: Reversal potential writer that writes to two different ions
+: as a simple function of internal concentration.
+
+NEURON {
+ SUFFIX write_multiple_eX
+ USEION x READ xi WRITE ex
+ USEION y READ yi WRITE ey
+}
+
+PARAMETER {}
+
+ASSIGNED {}
+
+STATE {}
+
+INITIAL {
+ ex = 100 + xi
+ ey = 200 + yi
+}
+
+BREAKPOINT {
+ ex = 100 + xi
+ ey = 200 + yi
+}
+
diff --git a/test/unit/test_fvm_layout.cpp b/test/unit/test_fvm_layout.cpp
index ae394dd5f0158f6c7f607bca015da74f44c12cb3..3bd981cd580c884e1dd9fcbaef2af08f4785fce1 100644
--- a/test/unit/test_fvm_layout.cpp
+++ b/test/unit/test_fvm_layout.cpp
@@ -94,23 +94,24 @@ namespace {
cable_cell c2;
segment* s;
+ c2.default_parameters.axial_resistivity = 90;
+
s = c2.add_soma(14./2);
s->add_mechanism("hh");
s = c2.add_cable(0, section_kind::dendrite, 1.0/2, 1.0/2, 200);
- s->cm = 0.017;
+ s->parameters.membrane_capacitance = 0.017;
s = c2.add_cable(1, section_kind::dendrite, 0.8/2, 0.8/2, 300);
- s->cm = 0.013;
+ s->parameters.membrane_capacitance = 0.013;
s = c2.add_cable(1, section_kind::dendrite, 0.7/2, 0.7/2, 180);
- s->cm = 0.018;
+ s->parameters.membrane_capacitance = 0.018;
c2.add_stimulus({2,1}, {5., 80., 0.45});
c2.add_stimulus({3,1}, {40., 10.,-0.2});
for (auto& seg: c2.segments()) {
- seg->rL = 90.;
if (seg->is_dendrite()) {
seg->add_mechanism("pas");
seg->set_compartments(4);
@@ -134,7 +135,7 @@ TEST(fvm_layout, topology) {
std::vector<cable_cell> cells = two_cell_system();
check_two_cell_system(cells);
- fvm_discretization D = fvm_discretize(cells);
+ fvm_discretization D = fvm_discretize(cells, neuron_parameter_defaults);
// Expected CV layouts for cells, segment indices in paren.
//
@@ -215,7 +216,7 @@ TEST(fvm_layout, area) {
std::vector<cable_cell> cells = two_cell_system();
check_two_cell_system(cells);
- fvm_discretization D = fvm_discretize(cells);
+ fvm_discretization D = fvm_discretize(cells, neuron_parameter_defaults);
// Note: stick models have constant diameter segments.
// Refer to comment above for CV vs. segment layout.
@@ -256,14 +257,13 @@ TEST(fvm_layout, area) {
// capacitance from segments 3, 4 and 5 (cell 1 segments
// 1, 2 and 3 respectively).
- double cm1 = cells[1].segment(1)->cm;
- double cm2 = cells[1].segment(2)->cm;
- double cm3 = cells[1].segment(3)->cm;
+ double cm1 = 0.017, cm2 = 0.013, cm3 = 0.018;
double c = A[3]/(2*n)*cm1+A[4]/(2*n)*cm2+A[5]/(2*n)*cm3;
EXPECT_FLOAT_EQ(c, D.cv_capacitance[11]);
- double cm0 = cells[1].soma()->cm;
+ //double cm0 = cells[1].soma()->cm;
+ double cm0 = neuron_parameter_defaults.membrane_capacitance.value();
c = A[2]*cm0;
EXPECT_FLOAT_EQ(c, D.cv_capacitance[6]);
@@ -276,8 +276,9 @@ TEST(fvm_layout, area) {
double a = volume(cable)/cable->length();
EXPECT_FLOAT_EQ(math::pi<double>*0.8*0.8/4, a);
+ double rL = 90;
double h = cable->length()/4;
- double g = a/h/cable->rL; // [µm·S/cm]
+ double g = a/h/rL; // [µm·S/cm]
g *= 100; // [µS]
EXPECT_FLOAT_EQ(g, D.face_conductance[13]);
@@ -294,7 +295,9 @@ TEST(fvm_layout, mech_index) {
cells[1].add_synapse({3, 0.4}, "expsyn");
cable_cell_global_properties gprop;
- fvm_discretization D = fvm_discretize(cells);
+ gprop.default_parameters = neuron_parameter_defaults;
+
+ fvm_discretization D = fvm_discretize(cells, gprop.default_parameters);
fvm_mechanism_data M = fvm_build_mechanism_data(gprop, cells, D);
auto& hh_config = M.mechanisms.at("hh");
@@ -353,9 +356,11 @@ TEST(fvm_layout, coalescing_synapses) {
};
cable_cell_global_properties gprop_no_coalesce;
+ gprop_no_coalesce.default_parameters = neuron_parameter_defaults;
gprop_no_coalesce.coalesce_synapses = false;
cable_cell_global_properties gprop_coalesce;
+ gprop_coalesce.default_parameters = neuron_parameter_defaults;
gprop_coalesce.coalesce_synapses = true;
{
@@ -367,7 +372,7 @@ TEST(fvm_layout, coalescing_synapses) {
cell.add_synapse({1, 0.7}, "expsyn");
cell.add_synapse({1, 0.9}, "expsyn");
- fvm_discretization D = fvm_discretize({cell});
+ fvm_discretization D = fvm_discretize({cell}, neuron_parameter_defaults);
fvm_mechanism_data M = fvm_build_mechanism_data(gprop_coalesce, {cell}, D);
auto &expsyn_config = M.mechanisms.at("expsyn");
@@ -383,7 +388,7 @@ TEST(fvm_layout, coalescing_synapses) {
cell.add_synapse({1, 0.7}, "expsyn");
cell.add_synapse({1, 0.9}, "exp2syn");
- fvm_discretization D = fvm_discretize({cell});
+ fvm_discretization D = fvm_discretize({cell}, neuron_parameter_defaults);
fvm_mechanism_data M = fvm_build_mechanism_data(gprop_coalesce, {cell}, D);
auto &expsyn_config = M.mechanisms.at("expsyn");
@@ -402,7 +407,7 @@ TEST(fvm_layout, coalescing_synapses) {
cell.add_synapse({1, 0.7}, "expsyn");
cell.add_synapse({1, 0.9}, "expsyn");
- fvm_discretization D = fvm_discretize({cell});
+ fvm_discretization D = fvm_discretize({cell}, neuron_parameter_defaults);
fvm_mechanism_data M = fvm_build_mechanism_data(gprop_no_coalesce, {cell}, D);
auto &expsyn_config = M.mechanisms.at("expsyn");
@@ -418,7 +423,7 @@ TEST(fvm_layout, coalescing_synapses) {
cell.add_synapse({1, 0.7}, "expsyn");
cell.add_synapse({1, 0.9}, "exp2syn");
- fvm_discretization D = fvm_discretize({cell});
+ fvm_discretization D = fvm_discretize({cell}, neuron_parameter_defaults);
fvm_mechanism_data M = fvm_build_mechanism_data(gprop_no_coalesce, {cell}, D);
auto &expsyn_config = M.mechanisms.at("expsyn");
@@ -438,7 +443,7 @@ TEST(fvm_layout, coalescing_synapses) {
cell.add_synapse({1, 0.7}, "expsyn");
cell.add_synapse({1, 0.7}, "expsyn");
- fvm_discretization D = fvm_discretize({cell});
+ fvm_discretization D = fvm_discretize({cell}, neuron_parameter_defaults);
fvm_mechanism_data M = fvm_build_mechanism_data(gprop_coalesce, {cell}, D);
auto &expsyn_config = M.mechanisms.at("expsyn");
@@ -454,7 +459,7 @@ TEST(fvm_layout, coalescing_synapses) {
cell.add_synapse({1, 0.3}, syn_desc("expsyn", 0.1, 0.2));
cell.add_synapse({1, 0.7}, syn_desc("expsyn", 0.1, 0.2));
- fvm_discretization D = fvm_discretize({cell});
+ fvm_discretization D = fvm_discretize({cell}, neuron_parameter_defaults);
fvm_mechanism_data M = fvm_build_mechanism_data(gprop_coalesce, {cell}, D);
auto &expsyn_config = M.mechanisms.at("expsyn");
@@ -476,7 +481,7 @@ TEST(fvm_layout, coalescing_synapses) {
cell.add_synapse({1, 0.3}, syn_desc("expsyn", 0, 2));
cell.add_synapse({1, 0.3}, syn_desc("expsyn", 1, 2));
- fvm_discretization D = fvm_discretize({cell});
+ fvm_discretization D = fvm_discretize({cell}, neuron_parameter_defaults);
fvm_mechanism_data M = fvm_build_mechanism_data(gprop_coalesce, {cell}, D);
auto &expsyn_config = M.mechanisms.at("expsyn");
@@ -501,7 +506,7 @@ TEST(fvm_layout, coalescing_synapses) {
cell.add_synapse({1, 0.7}, syn_desc_2("exp2syn", 2, 1));
cell.add_synapse({1, 0.7}, syn_desc_2("exp2syn", 2, 2));
- fvm_discretization D = fvm_discretize({cell});
+ fvm_discretization D = fvm_discretize({cell}, neuron_parameter_defaults);
fvm_mechanism_data M = fvm_build_mechanism_data(gprop_coalesce, {cell}, D);
auto &expsyn_config = M.mechanisms.at("expsyn");
@@ -548,7 +553,9 @@ TEST(fvm_layout, synapse_targets) {
cells[1].add_synapse({3, 0.7}, syn_desc("exp2syn", 6));
cable_cell_global_properties gprop;
- fvm_discretization D = fvm_discretize(cells);
+ gprop.default_parameters = neuron_parameter_defaults;
+
+ fvm_discretization D = fvm_discretize(cells, gprop.default_parameters);
fvm_mechanism_data M = fvm_build_mechanism_data(gprop, cells, D);
ASSERT_EQ(1u, M.mechanisms.count("expsyn"));
@@ -585,9 +592,6 @@ TEST(fvm_layout, synapse_targets) {
}
-// TODO: migrate tests for proportional parameter setting.
-
-
namespace {
double wm_impl(double wa, double xa) {
return wa? xa/wa: 0;
@@ -708,7 +712,9 @@ TEST(fvm_layout, density_norm_area) {
expected_gl[10] = seg3_gl;
cable_cell_global_properties gprop;
- fvm_discretization D = fvm_discretize(cells);
+ gprop.default_parameters = neuron_parameter_defaults;
+
+ fvm_discretization D = fvm_discretize(cells, gprop.default_parameters);
fvm_mechanism_data M = fvm_build_mechanism_data(gprop, cells, D);
// Check CV area assumptions.
@@ -747,19 +753,26 @@ TEST(fvm_layout, valence_verify) {
auto soma = c.add_soma(6);
soma->add_mechanism("test_cl_valence");
- fvm_discretization D = fvm_discretize(cells);
- mechanism_catalogue testcat = make_unit_test_catalogue();
cable_cell_global_properties gprop;
- gprop.catalogue = &testcat;
+ gprop.default_parameters = neuron_parameter_defaults;
- EXPECT_THROW(fvm_build_mechanism_data(gprop, cells, D), cable_cell_error);
+ fvm_discretization D = fvm_discretize(cells, neuron_parameter_defaults);
+
+ mechanism_catalogue testcat = make_unit_test_catalogue();
+ gprop.catalogue = &testcat;
- gprop.ion_default["cl"] = { -2, 1., 1. };
+ // Missing the 'cl' ion:
EXPECT_THROW(fvm_build_mechanism_data(gprop, cells, D), cable_cell_error);
- gprop.ion_default["cl"] = { -1, 1., 1. };
+ // Adding ion, should be fine now:
+ gprop.default_parameters.ion_data["cl"] = { 1., 1., 0. };
+ gprop.ion_species["cl"] = -1;
EXPECT_NO_THROW(fvm_build_mechanism_data(gprop, cells, D));
+
+ // 'cl' ion has wrong charge:
+ gprop.ion_species["cl"] = -2;
+ EXPECT_THROW(fvm_build_mechanism_data(gprop, cells, D), cable_cell_error);
}
TEST(fvm_layout, ion_weights) {
@@ -809,13 +822,24 @@ TEST(fvm_layout, ion_weights) {
{0}, {0, 2, 3}, {2, 3, 4}, {0, 1, 2, 3, 4}, {2, 4}
};
- fvec expected_iconc_norm_area[] = {
+ fvec expected_init_iconc[] = {
{0.}, {0., 1./2, 0.}, {1./4, 0., 0.}, {0., 0., 0., 0., 0.}, {3./4, 0.}
};
cable_cell_global_properties gprop;
+ gprop.default_parameters = neuron_parameter_defaults;
+
+ fvm_value_type cai = gprop.default_parameters.ion_data["ca"].init_int_concentration;
+ fvm_value_type cao = gprop.default_parameters.ion_data["ca"].init_ext_concentration;
+
+ for (auto& v: expected_init_iconc) {
+ for (auto& iconc: v) {
+ iconc *= cai;
+ }
+ }
for (auto run: count_along(mech_segs)) {
+ SCOPED_TRACE("run "+std::to_string(run));
std::vector<cable_cell> cells(1);
cable_cell& c = cells[0];
construct_cell(c);
@@ -824,7 +848,7 @@ TEST(fvm_layout, ion_weights) {
c.segments()[i]->add_mechanism("test_ca");
}
- fvm_discretization D = fvm_discretize(cells);
+ fvm_discretization D = fvm_discretize(cells, gprop.default_parameters);
fvm_mechanism_data M = fvm_build_mechanism_data(gprop, cells, D);
ASSERT_EQ(1u, M.ions.count("ca"s));
@@ -832,8 +856,84 @@ TEST(fvm_layout, ion_weights) {
EXPECT_EQ(expected_ion_cv[run], ca.cv);
- EXPECT_TRUE(testing::seq_almost_eq<fvm_value_type>(expected_iconc_norm_area[run], ca.iconc_norm_area));
+ EXPECT_TRUE(testing::seq_almost_eq<fvm_value_type>(expected_init_iconc[run], ca.init_iconc));
- EXPECT_TRUE(util::all_of(ca.econc_norm_area, [](fvm_value_type v) { return v==1.; }));
+ EXPECT_TRUE(util::all_of(ca.init_econc, [cao](fvm_value_type v) { return v==cao; }));
}
}
+
+TEST(fvm_layout, revpot) {
+ // Create two cells with three ions 'a', 'b' and 'c'.
+ // Configure a reversal potential mechanism that writes to 'a' and
+ // another that writes to 'b' and 'c'.
+ //
+ // Confirm:
+ // * Inconsistencies between revpot mech assignments are caught at discretization.
+ // * Reversal potential mechanisms are only extended where there exists another
+ // mechanism that reads them.
+
+ auto construct_cell = [](cable_cell& c) {
+ c.add_soma(5);
+
+ c.add_cable(0, section_kind::dendrite, 0.5, 0.5, 100);
+ c.add_cable(1, section_kind::dendrite, 0.5, 0.5, 200);
+ c.add_cable(1, section_kind::dendrite, 0.5, 0.5, 100);
+
+ for (auto& s: c.segments()) s->set_compartments(1);
+
+ // Read ea everywhere, ec only on soma.
+ for (auto& s: c.segments()) s->add_mechanism("read_eX/a");
+ c.soma()->add_mechanism("read_eX/c");
+ };
+
+ mechanism_catalogue testcat = make_unit_test_catalogue();
+
+ std::vector<cable_cell> cells(2);
+ for (auto& c: cells) construct_cell(c);
+
+ cable_cell_global_properties gprop;
+ gprop.default_parameters = neuron_parameter_defaults;
+ gprop.catalogue = &testcat;
+
+ gprop.ion_species = {{"a", 1}, {"b", 2}, {"c", 3}};
+ gprop.add_ion("a", 1, 10., 0, 0);
+ gprop.add_ion("b", 2, 30., 0, 0);
+ gprop.add_ion("c", 3, 50., 0, 0);
+
+ gprop.default_parameters.reversal_potential_method["a"] = "write_eX/a";
+ mechanism_desc write_eb_ec = "write_multiple_eX/x=b,y=c";
+
+ {
+ // need to specify ion "c" as well.
+ auto test_gprop = gprop;
+ test_gprop.default_parameters.reversal_potential_method["b"] = write_eb_ec;
+
+ fvm_discretization D = fvm_discretize(cells, test_gprop.default_parameters);
+ EXPECT_THROW(fvm_build_mechanism_data(test_gprop, cells, D), cable_cell_error);
+ }
+
+ {
+ // conflict with ion "c" on second cell.
+ auto test_gprop = gprop;
+ test_gprop.default_parameters.reversal_potential_method["b"] = write_eb_ec;
+ test_gprop.default_parameters.reversal_potential_method["c"] = write_eb_ec;
+ cells[1].default_parameters.reversal_potential_method["c"] = "write_eX/c";
+
+ fvm_discretization D = fvm_discretize(cells, test_gprop.default_parameters);
+ EXPECT_THROW(fvm_build_mechanism_data(test_gprop, cells, D), cable_cell_error);
+ }
+
+ auto& cell1_prop = cells[1].default_parameters;
+ cell1_prop.reversal_potential_method.clear();
+ cell1_prop.reversal_potential_method["b"] = write_eb_ec;
+ cell1_prop.reversal_potential_method["c"] = write_eb_ec;
+
+ fvm_discretization D = fvm_discretize(cells, gprop.default_parameters);
+ fvm_mechanism_data M = fvm_build_mechanism_data(gprop, cells, D);
+
+ // Only CV which needs write_multiple_eX/x=b,y=c is the soma (first CV)
+ // of the second cell.
+ auto soma1_index = D.cell_cv_bounds[1];
+ ASSERT_EQ(1u, M.mechanisms.count(write_eb_ec.name()));
+ EXPECT_EQ((std::vector<fvm_index_type>(1, soma1_index)), M.mechanisms.at(write_eb_ec.name()).cv);
+}
diff --git a/test/unit/test_fvm_lowered.cpp b/test/unit/test_fvm_lowered.cpp
index 7bd1b18cea57dddbb1b5d188178ff9bf44244551..42726b6b15511743a58fdec44518521cc4f23fc9 100644
--- a/test/unit/test_fvm_lowered.cpp
+++ b/test/unit/test_fvm_lowered.cpp
@@ -81,7 +81,6 @@ ACCESS_BIND(\
private_ion_index_table_ptr,\
&arb::multicore::mechanism::ion_index_table)
-
using namespace arb;
class gap_recipe_0: public recipe {
@@ -336,7 +335,10 @@ TEST(fvm_lowered, stimulus) {
// as during the stimulus windows.
std::vector<fvm_index_type> cell_to_intdom(cells.size(), 0);
- fvm_discretization D = fvm_discretize(cells);
+ cable_cell_global_properties gprop;
+ gprop.default_parameters = neuron_parameter_defaults;
+
+ fvm_discretization D = fvm_discretize(cells, gprop.default_parameters);
const auto& A = D.cv_area;
std::vector<target_handle> targets;
@@ -537,8 +539,8 @@ TEST(fvm_lowered, read_valence) {
rec.catalogue() = make_unit_test_catalogue();
rec.catalogue().derive("na_read_valence", "test_ca_read_valence", {}, {{"ca", "na"}});
- rec.catalogue().derive("cr_read_valence", "na_read_valence", {}, {{"na", "cr"}});
- rec.add_ion("cr", 7, 0, 0);
+ rec.catalogue().derive("cr_read_valence", "na_read_valence", {}, {{"na", "mn"}});
+ rec.add_ion("mn", 7, 0, 0, 0);
fvm_cell fvcell(context);
fvcell.initialize({0}, rec, cell_to_intdom, targets, probe_map);
@@ -649,12 +651,12 @@ TEST(fvm_lowered, weighted_write_ion) {
// - Soma (segment 0) plus three dendrites (1, 2, 3) meeting at a branch point.
// - Dendritic segments are given 1 compartments each.
//
- // /
- // d2
- // /
+ // /
+ // d2
+ // /
// s0-d1
- // \.
- // d3
+ // \.
+ // d3
//
// The CV corresponding to the branch point should comprise the terminal
// 1/2 of segment 1 and the initial 1/2 of segments 2 and 3.
@@ -666,7 +668,7 @@ TEST(fvm_lowered, weighted_write_ion) {
// dend 3: 100 µm long, 1 µm diameter cylinder
//
// The radius of the soma is chosen such that the surface area of soma is
- // the same as a 100µm dendrite, which makes it easier to describe the
+ // the same as a 100 µm dendrite, which makes it easier to describe the
// expected weights.
execution_context context;
@@ -689,26 +691,27 @@ TEST(fvm_lowered, weighted_write_ion) {
// Ca ion writer test_ca on CV 2 and 4 via segment 3:
c.segments()[3] ->add_mechanism("test_ca");
+ cable1d_recipe rec(c);
+ rec.add_ion("ca", 2, con_int, con_ext, 0.0);
+
std::vector<target_handle> targets;
std::vector<fvm_index_type> cell_to_intdom;
probe_association_map<probe_handle> probe_map;
fvm_cell fvcell(context);
- fvcell.initialize({0}, cable1d_recipe(c), cell_to_intdom, targets, probe_map);
+ fvcell.initialize({0}, rec, cell_to_intdom, targets, probe_map);
auto& state = *(fvcell.*private_state_ptr).get();
auto& ion = state.ion_data.at("ca"s);
- ion.default_int_concentration = con_int;
- ion.default_ext_concentration = con_ext;
ion.init_concentration();
std::vector<unsigned> ion_nodes = util::assign_from(ion.node_index_);
std::vector<unsigned> expected_ion_nodes = {2, 3, 4};
EXPECT_EQ(expected_ion_nodes, ion_nodes);
- std::vector<double> ion_iconc_weights = util::assign_from(ion.weight_Xi_);
- std::vector<double> expected_ion_iconc_weights = {0.75, 1., 0};
- EXPECT_EQ(expected_ion_iconc_weights, ion_iconc_weights);
+ std::vector<double> ion_init_iconc = util::assign_from(ion.init_Xi_);
+ std::vector<double> expected_init_iconc = {0.75*con_int, 1.*con_int, 0};
+ EXPECT_EQ(expected_init_iconc, ion_init_iconc);
auto test_ca = dynamic_cast<multicore::mechanism*>(find_mechanism(fvcell, "test_ca"));
@@ -721,13 +724,15 @@ TEST(fvm_lowered, weighted_write_ion) {
auto& test_ca_ca_index = *opt_ca_index_ptr.value();
double cai_contrib[3] = {200., 0., 300.};
+ double test_ca_weight[3] = {0.25, 0., 1.};
+
for (int i = 0; i<2; ++i) {
test_ca_cai[i] = cai_contrib[test_ca_ca_index[i]];
}
std::vector<double> expected_iconc(3);
for (int i = 0; i<3; ++i) {
- expected_iconc[i] = math::lerp(cai_contrib[i], con_int, ion_iconc_weights[i]);
+ expected_iconc[i] = test_ca_weight[i]*cai_contrib[i] + ion_init_iconc[i];
}
ion.init_concentration();
@@ -776,7 +781,7 @@ TEST(fvm_lowered, gj_coords_simple) {
d.add_gap_junction({1, 1});
cells.push_back(std::move(d));
- fvm_discretization D = fvm_discretize(cells);
+ fvm_discretization D = fvm_discretize(cells, neuron_parameter_defaults);
std::vector<cell_gid_type> gids = {0, 1};
auto GJ = fvcell.fvm_gap_junctions(cells, gids, rec, D);
@@ -886,7 +891,7 @@ TEST(fvm_lowered, gj_coords_complex) {
std::vector<cell_gid_type> gids = {0, 1, 2};
fvcell.fvm_intdom(rec, gids, cell_to_intdom);
- fvm_discretization D = fvm_discretize(cells);
+ fvm_discretization D = fvm_discretize(cells, neuron_parameter_defaults);
auto GJ = fvcell.fvm_gap_junctions(cells, gids, rec, D);
EXPECT_EQ(10u, GJ.size());
@@ -971,8 +976,8 @@ TEST(fvm_lowered, cell_group_gj) {
auto num_dom0 = fvcell.fvm_intdom(rec, gids_cg0, cell_to_intdom0);
auto num_dom1 = fvcell.fvm_intdom(rec, gids_cg1, cell_to_intdom1);
- fvm_discretization D0 = fvm_discretize(cell_group0);
- fvm_discretization D1 = fvm_discretize(cell_group1);
+ fvm_discretization D0 = fvm_discretize(cell_group0, neuron_parameter_defaults);
+ fvm_discretization D1 = fvm_discretize(cell_group1, neuron_parameter_defaults);
auto GJ0 = fvcell.fvm_gap_junctions(cell_group0, gids_cg0, rec, D0);
auto GJ1 = fvcell.fvm_gap_junctions(cell_group1, gids_cg1, rec, D1);
diff --git a/test/unit/test_maputil.cpp b/test/unit/test_maputil.cpp
index 96ebd6d0724f554fd637598464df46130235ce9d..4cb6ed06d547de60005bbf56d41715ea2d748330 100644
--- a/test/unit/test_maputil.cpp
+++ b/test/unit/test_maputil.cpp
@@ -188,4 +188,18 @@ TEST(maputil, binary_search_index) {
ASSERT_TRUE(opti);
EXPECT_EQ(x, ns[*opti]);
}
+
+ // With projection:
+ const char* ss[] = {"bit", "short", "small", "longer", "very long", "also long", "sesquipedalian"};
+ auto proj = [](const char* s) -> int { return std::strlen(s); };
+
+ for (int x: {1, 4, 7, 10, 20}) {
+ EXPECT_FALSE(binary_search_index(ss, x, proj));
+ }
+
+ for (int x: {3, 5, 6, 9, 14}) {
+ auto opti = binary_search_index(ss, x, proj);
+ ASSERT_TRUE(opti);
+ EXPECT_EQ(x, (int)std::strlen(ss[*opti]));
+ }
}
diff --git a/test/unit/test_mc_cell_group.cpp b/test/unit/test_mc_cell_group.cpp
index ed9cc54345e8ed05a7878877152ebd0268894c36..5a453dd6d6056979e14c311e97222f05c26ca330 100644
--- a/test/unit/test_mc_cell_group.cpp
+++ b/test/unit/test_mc_cell_group.cpp
@@ -42,7 +42,12 @@ TEST(mc_cell_group, get_kind) {
}
TEST(mc_cell_group, test) {
- mc_cell_group group{{0}, cable1d_recipe(make_cell()), lowered_cell()};
+ auto rec = cable1d_recipe(make_cell());
+ rec.nernst_ion("na");
+ rec.nernst_ion("ca");
+ rec.nernst_ion("k");
+
+ mc_cell_group group{{0}, rec, lowered_cell()};
group.advance(epoch(0, 50), 0.01, {});
// Model is expected to generate 4 spikes as a result of the
@@ -65,7 +70,12 @@ TEST(mc_cell_group, sources) {
}
std::vector<cell_gid_type> gids = {3u, 4u, 10u, 16u, 17u, 18u};
- mc_cell_group group{gids, cable1d_recipe(cells), lowered_cell()};
+ auto rec = cable1d_recipe(cells);
+ rec.nernst_ion("na");
+ rec.nernst_ion("ca");
+ rec.nernst_ion("k");
+
+ mc_cell_group group{gids, rec, lowered_cell()};
// Expect group sources to be lexicographically sorted by source id
// with gids in cell group's range and indices starting from zero.
diff --git a/test/unit/test_mech_temperature.cpp b/test/unit/test_mech_temperature.cpp
index 0fae8d43deea58d324a5b0807f958545175f6be1..1a616c76e6699bb4769ca224924a1ad79119acda 100644
--- a/test/unit/test_mech_temperature.cpp
+++ b/test/unit/test_mech_temperature.cpp
@@ -28,8 +28,14 @@ void run_celsius_test() {
auto instance = cat.instance<backend>("celsius_test");
auto& celsius_test = instance.mech;
+ double temperature_K = 300.;
+ double temperature_C = temperature_K-273.15;
+
+ std::vector<fvm_value_type> temp(ncv, temperature_K);
+ std::vector<fvm_value_type> vinit(ncv, -65);
+
auto shared_state = std::make_unique<typename backend::shared_state>(
- ncell, cv_to_intdom, gj, celsius_test->data_alignment());
+ ncell, cv_to_intdom, gj, vinit, temp, celsius_test->data_alignment());
mechanism_layout layout;
mechanism_overrides overrides;
@@ -40,11 +46,7 @@ void run_celsius_test() {
}
celsius_test->instantiate(0, *shared_state, overrides, layout);
-
- double temperature_K = 300.;
- double temperature_C = temperature_K-273.15;
-
- shared_state->reset(-65., temperature_K);
+ shared_state->reset();
// expect 0 value in state 'c' after init:
@@ -59,19 +61,6 @@ void run_celsius_test() {
expected_c_values.assign(ncv, temperature_C);
EXPECT_EQ(expected_c_values, mechanism_field(celsius_test.get(), "c"));
-
- // reset with new temperature and repeat test:
-
- temperature_K = 290.;
- temperature_C = temperature_K-273.15;
-
- shared_state->reset(-65., temperature_K);
- celsius_test->initialize();
-
- celsius_test->nrn_state();
- expected_c_values.assign(ncv, temperature_C);
-
- EXPECT_EQ(expected_c_values, mechanism_field(celsius_test.get(), "c"));
}
TEST(mech_temperature, celsius) {
diff --git a/test/unit/test_synapses.cpp b/test/unit/test_synapses.cpp
index 0efc50f57f0838573a34e0f40230749326f62dcd..56700e8352fb9311bb6f8a740ba03ab678cda8c3 100644
--- a/test/unit/test_synapses.cpp
+++ b/test/unit/test_synapses.cpp
@@ -72,13 +72,15 @@ auto unique_cast(std::unique_ptr<B> p) {
TEST(synapses, syn_basic_state) {
using util::fill;
- using value_type = multicore::backend::value_type;
- using index_type = multicore::backend::index_type;
+ using value_type = fvm_value_type;
+ using index_type = fvm_index_type;
int num_syn = 4;
int num_comp = 4;
int num_intdom = 1;
+ value_type temp_K = *neuron_parameter_defaults.temperature_K;
+
auto expsyn = unique_cast<multicore::mechanism>(global_default_catalogue().instance<backend>("expsyn").mech);
ASSERT_TRUE(expsyn);
@@ -87,9 +89,15 @@ TEST(synapses, syn_basic_state) {
std::vector<fvm_gap_junction> gj = {};
auto align = std::max(expsyn->data_alignment(), exp2syn->data_alignment());
- shared_state state(num_intdom, std::vector<index_type>(num_comp, 0), gj, align);
- state.reset(-65., constant::hh_squid_temp);
+ shared_state state(num_intdom,
+ std::vector<index_type>(num_comp, 0),
+ {},
+ std::vector<value_type>(num_comp, -65),
+ std::vector<value_type>(num_comp, temp_K),
+ align);
+
+ state.reset();
fill(state.current_density, 1.0);
fill(state.time_to, 0.1);
state.set_dt();
diff --git a/test/unit/unit_test_catalogue.cpp b/test/unit/unit_test_catalogue.cpp
index b12af8da2580092de82048956c9c750f4eed6266..fc4987662ab6ead860b7af4044fac0072851613d 100644
--- a/test/unit/unit_test_catalogue.cpp
+++ b/test/unit/unit_test_catalogue.cpp
@@ -13,6 +13,9 @@
#include "mechanisms/linear_ca_conc.hpp"
#include "mechanisms/test_cl_valence.hpp"
#include "mechanisms/test_ca_read_valence.hpp"
+#include "mechanisms/read_eX.hpp"
+#include "mechanisms/write_multiple_eX.hpp"
+#include "mechanisms/write_eX.hpp"
#include "../gtest.h"
@@ -38,6 +41,9 @@ mechanism_catalogue make_unit_test_catalogue() {
ADD_MECH(cat, linear_ca_conc)
ADD_MECH(cat, test_cl_valence)
ADD_MECH(cat, test_ca_read_valence)
+ ADD_MECH(cat, read_eX)
+ ADD_MECH(cat, write_multiple_eX)
+ ADD_MECH(cat, write_eX)
return cat;
}