diff --git a/arbor/fvm_layout.cpp b/arbor/fvm_layout.cpp
index dc6364325399b3f7a604a5d95a0317e4334bc241..5fbfafbe20e42f52e83b4ebc87adaa713d12bd43 100644
--- a/arbor/fvm_layout.cpp
+++ b/arbor/fvm_layout.cpp
@@ -592,7 +592,12 @@ fvm_mechanism_data fvm_build_mechanism_data(const cable_cell_global_properties&
assign(param_names, util::keys(info.parameters));
sort(param_names);
- for (auto& p: param_names) {
+ std::size_t n_param = param_names.size();
+ param_dflt.reserve(n_param);
+ config.param_values.reserve(n_param);
+
+ for (std::size_t i = 0; i<n_param; ++i) {
+ const auto& p = param_names[i];
config.param_values.emplace_back(p, std::vector<value_type>{});
param_dflt.push_back(info.parameters.at(p).default_value);
}
@@ -600,23 +605,21 @@ fvm_mechanism_data fvm_build_mechanism_data(const cable_cell_global_properties&
mcable_map<double> support;
std::vector<mcable_map<double>> param_maps;
- {
- std::unordered_map<std::string, mcable_map<double>> keyed_param_maps;
- for (auto& on_cable: entry.second) {
- verify_mechanism(info, on_cable.second);
- mcable cable = on_cable.first;
-
- support.insert(cable, 1.);
- for (auto param_assign: on_cable.second.values()) {
- keyed_param_maps[param_assign.first].insert(cable, param_assign.second);
- }
- }
- for (auto& p: param_names) {
- param_maps.push_back(std::move(keyed_param_maps[p]));
+ param_maps.resize(n_param);
+
+ for (auto& on_cable: entry.second) {
+ verify_mechanism(info, on_cable.second);
+ mcable cable = on_cable.first;
+ const auto& set_params = on_cable.second.values();
+
+ support.insert(cable, 1.);
+ for (std::size_t i = 0; i<n_param; ++i) {
+ double value = value_by_key(set_params, param_names[i]).value_or(param_dflt[i]);
+ param_maps[i].insert(cable, value);
}
}
- std::vector<double> param_on_cv(config.param_values.size());
+ std::vector<double> param_on_cv(n_param);
for (auto cv: D.geometry.cell_cvs(cell_idx)) {
double area = 0;
@@ -627,17 +630,18 @@ fvm_mechanism_data fvm_build_mechanism_data(const cable_cell_global_properties&
if (!area_on_cable) continue;
area += area_on_cable;
- for (auto i: count_along(param_on_cv)) {
- param_on_cv[i] += embedding.integrate_area(c.branch, pw_over_cable(param_maps[i], c, param_dflt[i]));
+ for (std::size_t i = 0; i<n_param; ++i) {
+ param_on_cv[i] += embedding.integrate_area(c.branch, pw_over_cable(param_maps[i], c, 0.));
}
}
if (area>0) {
- double oo_cv_area = 1./D.cv_area[cv];
config.cv.push_back(cv);
- config.norm_area.push_back(area*oo_cv_area);
+ config.norm_area.push_back(area/D.cv_area[cv]);
+
+ double oo_area = 1./area;
for (auto i: count_along(param_on_cv)) {
- config.param_values[i].second.push_back(param_on_cv[i]*oo_cv_area);
+ config.param_values[i].second.push_back(param_on_cv[i]*oo_area);
}
}
}
diff --git a/test/unit/test_fvm_layout.cpp b/test/unit/test_fvm_layout.cpp
index 2b6fda52e5be9b475dd139522ed3f0d5b186f6f3..bbdc25fd98ff962bb85476bc3a7b68b36e7f8463 100644
--- a/test/unit/test_fvm_layout.cpp
+++ b/test/unit/test_fvm_layout.cpp
@@ -36,27 +36,27 @@ namespace {
// Bulk resistivity: 90 Ω·cm
// capacitance:
// soma: 0.01 F/m² [default]
- // segment 1: 0.017 F/m²
- // segment 2: 0.013 F/m²
- // segment 3: 0.018 F/m²
+ // branch 1: 0.017 F/m²
+ // branch 2: 0.013 F/m²
+ // branch 3: 0.018 F/m²
//
// Soma diameter: 14 µm
// Some mechanisms: HH (default params)
//
- // Segment 1 diameter: 1 µm
- // Segment 1 length: 200 µm
+ // Branch 1 diameter: 1 µm
+ // Branch 1 length: 200 µm
//
- // Segment 2 diameter: 0.8 µm
- // Segment 2 length: 300 µm
+ // Branch 2 diameter: 0.8 µm
+ // Branch 2 length: 300 µm
//
- // Segment 3 diameter: 0.7 µm
- // Segment 3 length: 180 µm
+ // Branch 3 diameter: 0.7 µm
+ // Branch 3 length: 180 µm
//
// Dendrite mechanisms: passive (default params).
- // Stimulus at end of segment 2, amplitude 0.45.
- // Stimulus at end of segment 3, amplitude -0.2.
+ // Stimulus at end of branch 2, amplitude 0.45.
+ // Stimulus at end of branch 3, amplitude -0.2.
//
- // All dendrite segments with 4 compartments.
+ // All dendrite branches with 4 compartments.
soma_cell_builder builder(7.);
auto b1 = builder.add_branch(0, 200, 0.5, 0.5, 4, "dend");
@@ -117,13 +117,13 @@ TEST(fvm_layout, mech_index) {
EXPECT_EQ(mechanismKind::density, hh_config.kind);
EXPECT_EQ(ivec({0,6}), hh_config.cv);
- // Three expsyn synapses, two 0.4 along segment 1, and one 0.4 along segment 5.
+ // Three expsyn synapses, two 0.4 along branch 1, and one 0.4 along branch 5.
// These two synapses can be coalesced into 1 synapse
- // 0.4 along => second (non-parent) CV for segment.
+ // 0.4 along => second (non-parent) CV for branch.
EXPECT_EQ(ivec({3, 17}), expsyn_config.cv);
- // One exp2syn synapse, 0.4 along segment 4.
+ // One exp2syn synapse, 0.4 along branch 4.
EXPECT_EQ(ivec({13}), exp2syn_config.cv);
@@ -459,28 +459,26 @@ namespace {
TEST(fvm_layout, density_norm_area) {
// Test area-weighted linear combination of density mechanism parameters.
- // Create a cell with 4 segments:
+ // Create a cell with 4 branches:
// - Soma (branch 0) plus three dendrites (1, 2, 3) meeting at a branch point.
- // - HH mechanism on all segments.
+ // - HH mechanism on all branches.
// - Discretize with 3 CVs per non-soma branch, centred on forks.
//
// The CV corresponding to the branch point should comprise the terminal
- // 1/6 of segment 1 and the initial 1/6 of segments 2 and 3.
+ // 1/6 of branch 1 and the initial 1/6 of branches 2 and 3.
//
// The HH mechanism current density parameters ('gnabar', 'gkbar' and 'gl') are set
- // differently for each segment:
+ // differently for each branch:
//
// soma: all default values (gnabar = 0.12, gkbar = .036, gl = .0003)
- // segment 1: gl = .0002
- // segment 2: gkbar = .05
- // segment 3: gkbar = .07, gl = .0004
+ // branch 1: gl = .0002
+ // branch 2: gkbar = .05
+ // branch 3: gkbar = .07, gl = .0004
//
// Geometry:
- // segment 1: 100 µm long, 1 µm diameter cylinder.
- // segment 2: 200 µm long, diameter linear taper from 1 µm to 0.2 µm.
- // segment 3: 150 µm long, 0.8 µm diameter cylinder.
- //
- // Use divided compartment view on segments to compute area contributions.
+ // branch 1: 100 µm long, 1 µm diameter cylinder.
+ // branch 2: 200 µm long, diameter linear taper from 1 µm to 0.2 µm.
+ // branch 3: 150 µm long, 0.8 µm diameter cylinder.
soma_cell_builder builder(12.6157/2.0);
@@ -521,31 +519,31 @@ TEST(fvm_layout, density_norm_area) {
std::vector<double> expected_gkbar(ncv, dflt_gkbar);
std::vector<double> expected_gl(ncv, dflt_gl);
- // Last 1/6 of segment 1
+ // Last 1/6 of branch 1
double seg1_area_right = cells[0].embedding().integrate_area(mcable{1, 5./6., 1.});
- // First 1/6 of segment 2
+ // First 1/6 of branch 2
double seg2_area_left = cells[0].embedding().integrate_area(mcable{2, 0., 1./6.});
- // First 1/6 of segment 3
+ // First 1/6 of branch 3
double seg3_area_left = cells[0].embedding().integrate_area(mcable{3, 0., 1./6.});
// CV 0: soma
- // CV1: left of segment 1
+ // CV1: left of branch 1
expected_gl[0] = dflt_gl;
expected_gl[1] = seg1_gl;
expected_gl[2] = seg1_gl;
expected_gl[3] = seg1_gl;
- // CV 4: mix of right of segment 1 and left of segments 2 and 3.
+ // CV 4: mix of right of branch 1 and left of branches 2 and 3.
expected_gkbar[4] = wmean(seg1_area_right, dflt_gkbar, seg2_area_left, seg2_gkbar, seg3_area_left, seg3_gkbar);
expected_gl[4] = wmean(seg1_area_right, seg1_gl, seg2_area_left, dflt_gl, seg3_area_left, seg3_gl);
- // CV 5-7: just segment 2
+ // CV 5-7: just branch 2
expected_gkbar[5] = seg2_gkbar;
expected_gkbar[6] = seg2_gkbar;
expected_gkbar[7] = seg2_gkbar;
- // CV 8-10: just segment 3
+ // CV 8-10: just branch 3
expected_gkbar[8] = seg3_gkbar;
expected_gkbar[9] = seg3_gkbar;
expected_gkbar[10] = seg3_gkbar;
@@ -572,6 +570,91 @@ TEST(fvm_layout, density_norm_area) {
EXPECT_TRUE(testing::seq_almost_eq<double>(expected_gl, gl));
}
+TEST(fvm_layout, density_norm_area_partial) {
+ // Test area-weighted linear combination of density mechanism parameters,
+ // when mechanism covers only part of CV.
+
+ // Create a cell with 2 unbranched cables:
+ // - Soma (branch 0) plus one constant-diameter dendrite.
+ // - HH mechanism on part of the dendrite.
+ // - Discretize with 1 CV per branch.
+ //
+ // The HH mechanism is applied to the first 30% and last 60% of the dendrite:
+ //
+ // first 30%: all default values (gnabar = 0.12, gkbar = .036, gl = .0003)
+ // last 60%: gl = .0002, gkbar = .05
+ //
+ // Geometry:
+ // dendrite: 200 µm long, diameter linear taper from 1 µm to 0.2 µm.
+
+ soma_cell_builder builder(12.6157/2.0);
+
+ // p len r1 r2 ncomp tag
+ builder.add_branch(0, 200, 0.5, 0.1, 1, "dend");
+
+ double dflt_gnabar = .12;
+ double dflt_gkbar = .036;
+ double dflt_gl = 0.0003;
+
+ double end_gl = .0002;
+ double end_gkbar = .05;
+
+ auto hh_begin = mechanism_desc("hh");
+
+ auto hh_end = mechanism_desc("hh");
+ hh_end["gl"] = end_gl;
+ hh_end["gkbar"] = end_gkbar;
+
+ auto cell = builder.make_cell();
+ cell.default_parameters.discretization = cv_policy_fixed_per_branch(1);
+
+ cell.paint(mcable{1, 0., 0.3}, hh_begin);
+ cell.paint(mcable{1, 0.4, 1.}, hh_end);
+
+ std::vector<cable_cell> cells{std::move(cell)};
+
+ // First 30% of branch 1.
+ double b1_area_begin = cells[0].embedding().integrate_area(mcable{1, 0., 0.3});
+ // Last 60% of branch 1.
+ double b1_area_end = cells[0].embedding().integrate_area(mcable{1, 0.4, 1.});
+ // All of branch 1.
+ double b1_area = cells[0].embedding().integrate_area(mcable{1, 0., 1.});
+
+ double expected_norm_area = (b1_area_begin+b1_area_end)/b1_area;
+ double expected_gnabar = dflt_gnabar;
+ double expected_gkbar = (dflt_gkbar*b1_area_begin + end_gkbar*b1_area_end)/(b1_area_begin + b1_area_end);
+ double expected_gl = (dflt_gl*b1_area_begin + end_gl*b1_area_end)/(b1_area_begin + b1_area_end);
+
+ cable_cell_global_properties gprop;
+ gprop.default_parameters = neuron_parameter_defaults;
+
+ fvm_cv_discretization D = fvm_cv_discretize(cells, gprop.default_parameters);
+ fvm_mechanism_data M = fvm_build_mechanism_data(gprop, cells, D);
+
+ // Grab the HH parameters from the mechanism.
+
+ EXPECT_EQ(1u, M.mechanisms.size());
+ ASSERT_EQ(1u, M.mechanisms.count("hh"));
+
+ auto& norm_area = M.mechanisms.at("hh").norm_area;
+ ASSERT_EQ(1u, norm_area.size());
+ EXPECT_DOUBLE_EQ(expected_norm_area, norm_area[0]);
+
+ auto& hh_params = M.mechanisms.at("hh").param_values;
+
+ auto& gkbar = value_by_key(hh_params, "gkbar"s).value();
+ auto& gnabar = value_by_key(hh_params, "gnabar"s).value();
+ auto& gl = value_by_key(hh_params, "gl"s).value();
+
+ ASSERT_EQ(1u, gkbar.size());
+ ASSERT_EQ(1u, gnabar.size());
+ ASSERT_EQ(1u, gl.size());
+
+ EXPECT_DOUBLE_EQ(expected_gkbar, gkbar[0]);
+ EXPECT_DOUBLE_EQ(expected_gnabar, gnabar[0]);
+ EXPECT_DOUBLE_EQ(expected_gl, gl[0]);
+}
+
TEST(fvm_layout, valence_verify) {
auto cell = soma_cell_builder(6).make_cell();
cell.paint("soma", "test_cl_valence");
@@ -599,9 +682,9 @@ TEST(fvm_layout, valence_verify) {
}
TEST(fvm_layout, ion_weights) {
- // Create a cell with 4 segments:
- // - Soma (segment 0) plus three dendrites (1, 2, 3) meeting at a branch point.
- // - Dendritic segments are given 1 compartments each.
+ // Create a cell with 4 branches:
+ // - Soma (branch 0) plus three dendrites (1, 2, 3) meeting at a branch point.
+ // - Dendritic branches are given 1 compartments each.
//
// /
// d2
@@ -611,7 +694,7 @@ TEST(fvm_layout, ion_weights) {
// 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.
+ // 1/2 of branch 1 and the initial 1/2 of branches 2 and 3.
//
// Geometry:
// soma 0: radius 5 µm