diff --git a/example/CMakeLists.txt b/example/CMakeLists.txt
index 0c33de290175c5ebca1f2b0a5b7e04f268679bea..7a665e5c3b9fe66f3e7fea1282ae18ce7e430456 100644
--- a/example/CMakeLists.txt
+++ b/example/CMakeLists.txt
@@ -10,3 +10,4 @@ add_subdirectory(ring)
add_subdirectory(gap_junctions)
add_subdirectory(single)
add_subdirectory(probe-demo)
+add_subdirectory(lfp)
diff --git a/example/lfp/CMakeLists.txt b/example/lfp/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..a3046c48121e81509f63aeb205eef364b559a3f9
--- /dev/null
+++ b/example/lfp/CMakeLists.txt
@@ -0,0 +1,4 @@
+add_executable(lfp EXCLUDE_FROM_ALL lfp.cpp)
+add_dependencies(examples lfp)
+target_link_libraries(lfp PRIVATE arbor ext-tinyopt)
+file(COPY plot-lfp.py DESTINATION "${CMAKE_RUNTIME_OUTPUT_DIRECTORY}")
diff --git a/example/lfp/README.md b/example/lfp/README.md
new file mode 100644
index 0000000000000000000000000000000000000000..a185549c527f88c5a78d1c206266f27cafd4008b
--- /dev/null
+++ b/example/lfp/README.md
@@ -0,0 +1,30 @@
+# Local field potential demo
+
+How might one use Arbor to compute the local field potential near a cell?
+
+This example provide a simple demonstration, simulating one cell and computing
+the LFP from the total membrane current.
+
+The code attempts to provide an Arbor version of the supplied NEURON LFP example
+`neuron_lfp_example.py`. The plot from the NEURON code is included as `example_nrn_EP.png`.
+
+## How to run
+
+The example builds an executable `lfp` that performs a 100 ms simulation of
+a single ball-and-stick neuron, measuring the local field potential at two
+electrode sites.
+
+Running `lfp` generates a JSON file with the simulation output written to stdout.
+The included `plot-lfp.py` script will parse this output and generate a plot.
+
+Run `lfp` and display the output in a window:
+```
+lfp | plot-lfp.py
+```
+
+Run `lfp` and save the results, then generate a plot image:
+```
+lfp > out.json
+plot-lfp.py -o out.png out.json
+```
+
diff --git a/example/lfp/example_nrn_EP.png b/example/lfp/example_nrn_EP.png
new file mode 100644
index 0000000000000000000000000000000000000000..d3b3fa042ddaf518eae6271953484ce33e378c1d
Binary files /dev/null and b/example/lfp/example_nrn_EP.png differ
diff --git a/example/lfp/lfp.cpp b/example/lfp/lfp.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..268fd286f283d8238d67d3e873e2a8fa4b60688f
--- /dev/null
+++ b/example/lfp/lfp.cpp
@@ -0,0 +1,303 @@
+#include <cassert>
+#include <vector>
+#include <iostream>
+
+#include <arbor/load_balance.hpp>
+#include <arbor/cable_cell.hpp>
+#include <arbor/morph/morphology.hpp>
+#include <arbor/morph/place_pwlin.hpp>
+#include <arbor/morph/region.hpp>
+#include <arbor/simple_sampler.hpp>
+#include <arbor/simulation.hpp>
+#include <arbor/sampling.hpp>
+#include <arbor/util/any.hpp>
+#include <arbor/util/any_ptr.hpp>
+
+using arb::util::any;
+using arb::util::any_cast;
+using arb::util::any_ptr;
+using arb::cell_gid_type;
+
+// Recipe represents one cable cell with one synapse, together with probes for total trans-membrane current, membrane voltage,
+// ionic current density, and synaptic conductance. A sequence of spikes are presented to the one synapse on the cell.
+
+struct lfp_demo_recipe: public arb::recipe {
+ explicit lfp_demo_recipe(arb::event_generator events): events_(std::move(events))
+ {
+ make_cell(); // initializes cell_ and synapse_location_.
+ }
+
+ arb::cell_size_type num_cells() const override { return 1; }
+ arb::cell_size_type num_targets(cell_gid_type) const override { return 1; }
+
+ std::vector<arb::probe_info> get_probes(cell_gid_type) const override {
+ // Four probes:
+ // 0. Total membrane current across cell.
+ // 1. Voltage at synapse location.
+ // 2. Total ionic current density at synapse location.
+ // 3. Expsyn synapse conductance value.
+ return {
+ arb::cable_probe_total_current_cell{},
+ arb::cable_probe_membrane_voltage{synapse_location_},
+ arb::cable_probe_total_ion_current_density{synapse_location_},
+ arb::cable_probe_point_state{0, "expsyn", "g"}};
+ }
+
+ arb::cell_kind get_cell_kind(cell_gid_type) const override {
+ return arb::cell_kind::cable;
+ }
+
+ arb::util::unique_any get_cell_description(cell_gid_type) const override {
+ return cell_;
+ }
+
+ virtual std::vector<arb::event_generator> event_generators(cell_gid_type) const {
+ return {events_};
+ }
+
+ any get_global_properties(arb::cell_kind) const {
+ arb::cable_cell_global_properties gprop;
+ gprop.default_parameters = arb::neuron_parameter_defaults;
+ return gprop;
+ }
+
+private:
+ arb::cable_cell cell_;
+ arb::locset synapse_location_;
+ arb::event_generator events_;
+
+ void make_cell() {
+ using namespace arb;
+
+ // Set up morphology as two branches:
+ // * soma, length 20 μm radius 10 μm, with SWC tag 1.
+ // * apical dendrite, length 490 μm, radius 1 μm, with SWC tag 4.
+ sample_tree tree;
+ tree.append({{0, 0, +10, 10}, 1}); // (root point)
+ tree.append({{0, 0, -10, 10}, 1});
+ tree.append(0, {{0, 0, 10, 1}, 4}); // attach to root point.
+ tree.append({{0, 0, 500, 1}, 4});
+
+ cell_ = cable_cell(tree);
+
+ // Use NEURON defaults for reversal potentials, ion concentrations etc., but override ra, cm.
+ cell_.default_parameters.axial_resistivity = 100; // [Ω·cm]
+ cell_.default_parameters.membrane_capacitance = 0.01; // [F/m²]
+
+ // Twenty CVs per branch, except for the soma.
+ cell_.default_parameters.discretization = cv_policy_fixed_per_branch(20, cv_policy_flag::single_root_cv);
+
+ // Add pas and hh mechanisms:
+ cell_.paint(reg::tagged(1), "hh"); // (default parameters)
+ cell_.paint(reg::tagged(4), mechanism_desc("pas").set("e", -70));
+
+ // Add exponential synapse at centre of soma (0.5 along branch 0).
+ synapse_location_ = mlocation{0, 0.5};
+ cell_.place(synapse_location_, mechanism_desc("expsyn").set("e", 0).set("tau", 2));
+ }
+};
+
+struct position { double x, y, z; };
+
+struct lfp_sampler {
+ lfp_sampler(const arb::place_pwlin& p, std::vector<position> electrodes, double sigma):
+ placement(p), electrodes(std::move(electrodes)), sigma(sigma) {}
+
+ // Compute response coefficients for each electrode, given a set of cable-like sources.
+ void initialize(const arb::mcable_list& cables) {
+ const unsigned n_electrode = electrodes.size();
+ response.assign(n_electrode, std::vector<double>(cables.size()));
+
+ std::vector<arb::mpoint> midpoints;
+ std::transform(cables.begin(), cables.end(), std::back_inserter(midpoints),
+ [this](const auto& c) { return placement.at({c.branch, 0.5*(c.prox_pos+c.dist_pos)}); });
+
+ const double coef = 1/(4*M_PI*sigma); // [Ω·m]
+ for (unsigned i = 0; i<n_electrode; ++i) {
+ const position& e = electrodes[i];
+
+ std::transform(midpoints.begin(), midpoints.end(), response[i].begin(),
+ [coef, &e](auto p) {
+ p.x -= e.x;
+ p.y -= e.y;
+ p.z -= e.z;
+ double r = std::sqrt(p.x*p.x+p.y*p.y+p.z*p.z); // [μm]
+ return coef/r; // [MΩ]
+ });
+ }
+ }
+
+ void reset() {
+ response.clear();
+ lfp_time.clear();
+ lfp_voltage.clear();
+ }
+
+ bool is_initialized() const {
+ return !response.empty();
+ }
+
+ // On receipt of a sequence of cell-wide current samples, apply response matrix and save results to lfp_voltage.
+ arb::sampler_function callback() {
+ return [this](arb::probe_metadata pm, std::size_t n, const arb::sample_record* samples) {
+ auto cables_ptr = any_cast<const arb::mcable_list*>(pm.meta);
+ assert(cables_ptr);
+
+ if (!is_initialized()) {
+ // The first time we get metadata, build the response matrix.
+ initialize(*cables_ptr);
+ }
+
+ std::vector<double> currents;
+ lfp_voltage.resize(response.size());
+
+ for (std::size_t i = 0; i<n; ++i) {
+ lfp_time.push_back(samples[i].time);
+
+ auto data_ptr = any_cast<const arb::cable_sample_range*>(samples[i].data);
+ assert(data_ptr);
+
+ for (unsigned j = 0; j<response.size(); ++j) {
+ lfp_voltage[j].push_back(std::inner_product(data_ptr->first, data_ptr->second, response[j].begin(), 0.));
+ }
+ }
+ };
+ }
+
+ std::vector<double> lfp_time;
+ std::vector<std::vector<double>> lfp_voltage; // [mV] (one vector per electrode)
+
+private:
+ const arb::place_pwlin placement; // Represents cell morphology in space.
+ const std::vector<position> electrodes; // [μm]
+ const double sigma; // [S/m]
+ std::vector<std::vector<double>> response; // [MΩ]
+};
+
+// JSON output helpers:
+
+template <typename T, typename F>
+struct as_json_array_wrap {
+ const T& data;
+ F fn;
+ as_json_array_wrap(const T& data, const F& fn): data(data), fn(fn) {}
+
+ friend std::ostream& operator<<(std::ostream& out, const as_json_array_wrap& a) {
+ out << '[';
+ bool first = true;
+ for (auto& x: a.data) out << (!first? ", ": (first=false, "")) << a.fn(x);
+ return out << ']';
+ }
+};
+
+struct {
+ template <typename F>
+ auto operator()(const F& fn) const {
+ return [&fn](const auto& data) { return as_json_array_wrap<decltype(data), F>(data, fn); };
+ }
+
+ auto operator()() const {
+ return this->operator()([](const auto& x) { return x; });
+ }
+} as_json_array;
+
+// Run simulation.
+
+int main(int argc, char** argv) {
+ auto context = arb::make_context();
+
+ // Weight 0.005 μS, onset at t = 0 ms, mean frequency 0.1 kHz.
+ auto events = arb::poisson_generator({0, 0}, .005, 0., 0.1, std::minstd_rand{});
+ lfp_demo_recipe R(events);
+
+ const double t_stop = 100; // [ms]
+ const double sample_dt = 0.1; // [ms]
+ const double dt = 0.1; // [ms]
+
+ arb::simulation sim(R, arb::partition_load_balance(R, context), context);
+
+ std::vector<position> electrodes = {
+ {30, 0, 0},
+ {30, 0, 100}
+ };
+
+ auto sample_schedule = arb::regular_schedule(sample_dt);
+
+ arb::morphology cell_morphology = any_cast<arb::cable_cell>(R.get_cell_description(0)).morphology();
+ arb::place_pwlin placed_cell(cell_morphology);
+ lfp_sampler lfp(placed_cell, electrodes, 3.0);
+ sim.add_sampler(arb::one_probe({0, 0}), sample_schedule, lfp.callback(), arb::sampling_policy::exact);
+
+ arb::trace_vector<double, arb::mlocation> membrane_voltage;
+ sim.add_sampler(arb::one_probe({0, 1}), sample_schedule, make_simple_sampler(membrane_voltage), arb::sampling_policy::exact);
+
+ arb::trace_vector<double> ionic_current_density;
+ sim.add_sampler(arb::one_probe({0, 2}), sample_schedule, make_simple_sampler(ionic_current_density), arb::sampling_policy::exact);
+
+ arb::trace_vector<double> synapse_g;
+ sim.add_sampler(arb::one_probe({0, 3}), sample_schedule, make_simple_sampler(synapse_g), arb::sampling_policy::exact);
+
+ sim.run(t_stop, dt);
+
+ // Output results in JSON format suitable for plotting by plot-lfp.py script.
+
+ auto get_t = [](const auto& x) { return x.t; };
+ auto get_v = [](const auto& x) { return x.v; };
+ auto scale = [](double s) { return [s](const auto& x) { return x*s; }; };
+ auto to_xz = [](const auto& p) { return std::array<double, 2>{p.x, p.z}; };
+
+ // Compute synaptic current from synapse conductance and membrane potential.
+ std::vector<double> syn_i;
+ assert(synapse_g.get(0).size()==membrane_voltage.get(0).size());
+ std::transform(synapse_g.get(0).begin(), synapse_g.get(0).end(), membrane_voltage.get(0).begin(), std::back_inserter(syn_i),
+ [](arb::trace_entry<double> g, arb::trace_entry<double> v) {
+ assert(g.t==v.t);
+ return g.v*v.v;
+ });
+
+ // Collect points from 2-d morphology in vectors of arrays (x, z, radius), one per branch.
+ // (This process will be simplified with improvements to the place_pwlin API.)
+ std::vector<std::vector<std::array<double, 3>>> samples;
+ for (unsigned branch = 0; branch<cell_morphology.num_branches(); ++branch) {
+ samples.push_back({});
+ auto branch_range = cell_morphology.branch_indexes(branch);
+ for (auto i_ptr = branch_range.first; i_ptr!=branch_range.second; ++i_ptr) {
+ arb::msample s = cell_morphology.samples()[*i_ptr];
+ samples.back().push_back(std::array<double, 3>{s.loc.x, s.loc.z, s.loc.radius});
+ }
+ }
+
+ auto probe_xz = to_xz(placed_cell.at(membrane_voltage.get(0).meta));
+ std::vector<std::array<double, 2>> electrodes_xz;
+ std::transform(electrodes.begin(), electrodes.end(), std::back_inserter(electrodes_xz), to_xz);
+
+ std::cout <<
+ "{\n"
+ "\"morphology\": {\n"
+ "\"unit\": \"μm\",\n"
+ "\"samples\": " << as_json_array(as_json_array(as_json_array()))(samples) << ",\n"
+ "\"probe\": " << as_json_array()(probe_xz) << ",\n"
+ "\"electrodes\": " << as_json_array(as_json_array())(electrodes_xz) << "\n"
+ "},\n"
+ "\"extracellular potential\": {\n"
+ "\"unit\": \"μV\",\n"
+ "\"time\": " << as_json_array()(lfp.lfp_time) << ",\n"
+ "\"values\": " << as_json_array(as_json_array(scale(1e3)))(lfp.lfp_voltage) << "\n"
+ "},\n"
+ "\"synaptic current\": {\n"
+ "\"unit\": \"nA\",\n"
+ "\"time\": " << as_json_array(get_t)(synapse_g.get(0)) << ",\n"
+ "\"value\": " << as_json_array()(syn_i) << "\n"
+ "},\n"
+ "\"membrane potential\": {\n"
+ "\"unit\": \"mV\",\n"
+ "\"time\": " << as_json_array(get_t)(membrane_voltage.get(0)) << ",\n"
+ "\"value\": " << as_json_array(get_v)(membrane_voltage.get(0)) << "\n"
+ "},\n"
+ "\"ionic current density\": {\n"
+ "\"unit\": \"A/m²\",\n"
+ "\"time\": " << as_json_array(get_t)(ionic_current_density.get(0)) << ",\n"
+ "\"value\": " << as_json_array(get_v)(ionic_current_density.get(0)) << "\n"
+ "}\n"
+ "}\n";
+}
diff --git a/example/lfp/neuron_lfp_example.py b/example/lfp/neuron_lfp_example.py
new file mode 100755
index 0000000000000000000000000000000000000000..1d1023ac427c3afc1688311e783f0e2388a71de1
--- /dev/null
+++ b/example/lfp/neuron_lfp_example.py
@@ -0,0 +1,296 @@
+#!/usr/env/bin python
+# -*- coding: utf-8 -*-
+# Author: Torbjørn Ness <torbjorn.ness@nmbu.no>
+'''
+NEURON and Python - Creating a multi-compartment model with synaptic input
+with randomized activation times
+'''
+# Import modules for plotting and NEURON itself
+import matplotlib.pyplot as plt
+import neuron
+import numpy as np
+
+
+class Cell:
+ """
+ Cell class that handles interactions with NEURON. It finds the center position
+ of each cellular compartment (cell.xmid, cell.ymid, cell.zmid), and the transmembrane currents
+ cell.imem
+ """
+
+ def __init__(self):
+ cvode = neuron.h.CVode()
+ cvode.use_fast_imem(1)
+ self.tstop = 100. # simulation duration in ms
+ self.v_init = -65 # membrane voltage(s) at t = 0
+
+ neuron.h.dt = 0.1
+
+ self.make_cell()
+ neuron.h.define_shape()
+ self.seclist = []
+ counter = 0
+ for sec in neuron.h.allsec():
+ self.seclist.append(sec)
+ for seg in sec:
+ counter += 1
+ self.totnsegs = counter # Total number of compartments in cell model
+ self.collect_geometry()
+
+
+ self.insert_synapse(self.seclist[0])
+
+ self.initiate_recorders()
+
+ def make_cell(self):
+ neuron.h("""
+ create soma[1]
+ create apic[1]
+ objref all
+ all = new SectionList()
+ soma[0] {pt3dclear()
+ pt3dadd(0, 0, -10, 20)
+ pt3dadd(0, 0, 10, 20)}
+
+ apic[0] {pt3dclear()
+ pt3dadd(0, 0, 10, 2)
+ pt3dadd(0, 0, 500, 2)}
+
+ connect apic[0](0), soma[0](1)
+
+ apic[0] {nseg = 20}
+ forall {
+ Ra = 100.
+ cm = 1.
+ all.append()
+ }
+ apic[0] { insert pas }
+ soma[0] { insert hh }
+
+ """)
+
+ def initiate_recorders(self):
+ self.imem = [] # Record membrane currents
+ self.vmem = [] # Record membrane potentials
+
+ for sec in neuron.h.allsec():
+ for seg in sec:
+ v_ = neuron.h.Vector()
+ v_.record(seg._ref_v, neuron.h.dt)
+ self.vmem.append(v_)
+
+ i_ = neuron.h.Vector()
+ i_.record(seg._ref_i_membrane_, neuron.h.dt)
+ self.imem.append(i_)
+
+ self.tvec = neuron.h.Vector()
+ self.tvec.record(neuron.h._ref_t)
+
+ if hasattr(self, "syn"):
+ self.syn_i = neuron.h.Vector() # Record synaptic current
+ self.syn_i.record(self.syn._ref_i, neuron.h.dt)
+
+ def insert_synapse(self, syn_sec):
+ """
+ Function to insert a single synapse into cell model, as section syn_sec
+ """
+ print(syn_sec.diam)
+ syn = neuron.h.ExpSyn(0.5, sec=syn_sec)
+ syn.e = 0. # reversal potential of synapse conductance in mV
+ syn.tau = 2. # time constant of synapse conductance in ms
+
+ ns = neuron.h.NetStim(0.5) # spike time generator object (~presynaptic)
+ ns.noise = 1. # Fractional randomness (intervals from exp dist)
+ ns.start = 0. # approximate time of first spike
+ ns.number = 1000 # number of spikes
+ ns.interval = 10. # average interspike interval
+ nc = neuron.h.NetCon(ns, syn) # Connect generator to synapse
+ nc.weight[0] = .005 # Set synapse weight
+
+ # Everything must be stored or NEURON will forget they ever existed
+ self.ns = ns
+ self.nc = nc
+ self.syn = syn
+
+ def collect_geometry(self):
+ """
+ Function to get positions, diameters etc of each segment in NEURON
+ """
+
+ areavec = np.zeros(self.totnsegs)
+ diamvec = np.zeros(self.totnsegs)
+ lengthvec = np.zeros(self.totnsegs)
+
+ xstartvec = np.zeros(self.totnsegs)
+ xendvec = np.zeros(self.totnsegs)
+ ystartvec = np.zeros(self.totnsegs)
+ yendvec = np.zeros(self.totnsegs)
+ zstartvec = np.zeros(self.totnsegs)
+ zendvec = np.zeros(self.totnsegs)
+
+ counter = 0
+
+ #loop over all segments
+ for sec in neuron.h.allsec():
+ n3d = int(neuron.h.n3d())
+ nseg = sec.nseg
+ gsen2 = 1./2/nseg
+ if n3d > 0:
+ #create interpolation objects for the xyz pt3d info:
+ L = np.zeros(n3d)
+ x = np.zeros(n3d)
+ y = np.zeros(n3d)
+ z = np.zeros(n3d)
+ for i in range(n3d):
+ L[i] = neuron.h.arc3d(i)
+ x[i] = neuron.h.x3d(i)
+ y[i] = neuron.h.y3d(i)
+ z[i] = neuron.h.z3d(i)
+
+ #normalize as seg.x [0, 1]
+ L /= sec.L
+
+ #temporary store position of segment midpoints
+ segx = np.zeros(nseg)
+ for i, seg in enumerate(sec):
+ segx[i] = seg.x
+
+ #can't be >0 which may happen due to NEURON->Python float transfer:
+ segx0 = (segx - gsen2).round(decimals=6)
+ segx1 = (segx + gsen2).round(decimals=6)
+
+ #fill vectors with interpolated coordinates of start and end points
+ xstartvec[counter:counter+nseg] = np.interp(segx0, L, x)
+ xendvec[counter:counter+nseg] = np.interp(segx1, L, x)
+
+ ystartvec[counter:counter+nseg] = np.interp(segx0, L, y)
+ yendvec[counter:counter+nseg] = np.interp(segx1, L, y)
+
+ zstartvec[counter:counter+nseg] = np.interp(segx0, L, z)
+ zendvec[counter:counter+nseg] = np.interp(segx1, L, z)
+
+ #fill in values area, diam, length
+ for i, seg in enumerate(sec):
+ areavec[counter] = neuron.h.area(seg.x)
+ diamvec[counter] = seg.diam
+ lengthvec[counter] = sec.L/nseg
+ counter += 1
+
+ # starting position of each compartment (segment)
+ self.xstart = xstartvec
+ self.ystart = ystartvec
+ self.zstart = zstartvec
+
+ # ending position of each compartment (segment)
+ self.xend = xendvec
+ self.yend = yendvec
+ self.zend = zendvec
+
+ # Calculates the center position of each compartment (segment)
+ self.xmid = 0.5 * (self.xstart + self.xend)
+ self.ymid = 0.5 * (self.ystart + self.yend)
+ self.zmid = 0.5 * (self.zstart + self.zend)
+ self.area = areavec
+ self.diam = diamvec
+
+ def simulate(self):
+
+ neuron.h.finitialize(self.v_init)
+ neuron.h.fcurrent()
+
+ while neuron.h.t < self.tstop:
+ neuron.h.fadvance()
+
+ self.vmem = np.array(self.vmem)
+ self.imem = np.array(self.imem)
+ self.syn_i = np.array(self.syn_i)
+ self.tvec = np.array(self.tvec)[:self.vmem.shape[1]]
+
+
+class ExtElectrode:
+ def __init__(self, elec_x, elec_y, elec_z):
+ """
+
+ :param elec_x, elec_y, elec_z : x,y,z-positions (um) of each electrode. Must
+ be numpy arrays of equal length
+
+ """
+ self.sigma = 0.3 # Extracellular conductivity (S/m)
+ self.elec_x = elec_x
+ self.elec_y = elec_y
+ self.elec_z = elec_z
+ self.num_elecs = len(self.elec_x)
+
+ # Give electrodes different colors for plotting purposes:
+ self.elec_clr = lambda idx: plt.cm.viridis(idx / self.num_elecs)
+
+ def calc_extracellular_potential(self, cell):
+ self.calc_mapping(cell)
+ self.extracellular_potential = np.dot(electrode.mapping, cell.imem)
+
+ def calc_mapping(self, cell):
+ """
+ Calculates 'mapping' of size (number of electrodes) * (number of cell compartments)
+ Extracellular potential can then be calculated as
+ :param cell: class containing x,y,z-positions (um) of each
+ compartment (segment) center as cell.xmid, cell.ymid, cell.zmid
+
+ """
+ self.mapping = np.zeros((self.num_elecs, cell.totnsegs))
+ for e_idx in range(self.num_elecs):
+ r2 = ((cell.xmid - self.elec_x[e_idx])**2 +
+ (cell.ymid - self.elec_y[e_idx])**2 +
+ (cell.zmid - self.elec_z[e_idx])**2)
+
+ self.mapping[e_idx] = 1 / (4 * np.pi * self.sigma * np.sqrt(r2))
+
+
+def plot_results(cell, electrode):
+ ################################################################################
+ # Plot simulated output
+ ################################################################################
+ fig = plt.figure(figsize=(9, 5))
+ fig.subplots_adjust(wspace=0.5, hspace=0.9)
+ ax_morph = fig.add_subplot(131, aspect=1, xlim=[-150, 150], ylim=[-100, 600],
+ title="morphology", xlabel="x ($\mu$m)", ylabel="y ($\mu$m)")
+ ax_syn = fig.add_subplot(332, ylabel='nA', title="synaptic current", xlabel='time (ms)')
+ ax_vmem = fig.add_subplot(335, ylabel='mV', xlabel='time (ms)', title="membrane potential")
+ ax_imem = fig.add_subplot(338, ylabel='nA', xlabel='time (ms)', title="membrane current")
+ ax_ep = fig.add_subplot(133, ylabel='$\mu$V', xlabel="time (ms)", title="Extracellular potential")
+
+ plot_comp_idx = 0
+ plot_comp_clr = 'r'
+
+ for idx in range(cell.totnsegs):
+ ax_morph.plot([cell.xstart[idx], cell.xend[idx]],
+ [cell.zstart[idx], cell.zend[idx]], lw=cell.diam[idx] / 2, c='k')
+ ax_morph.plot(cell.xmid[plot_comp_idx], cell.zmid[plot_comp_idx], marker='*', c=plot_comp_clr)
+
+ ax_syn.plot(cell.tvec, cell.syn_i, c='k', lw=2)
+ ax_vmem.plot(cell.tvec, cell.vmem[0, :], c=plot_comp_clr, lw=2)
+ ax_imem.plot(cell.tvec, cell.imem[0, :], c=plot_comp_clr, lw=2)
+
+ for e_idx in range(electrode.num_elecs):
+ e_clr = electrode.elec_clr(e_idx)
+ sig = 1000 * electrode.extracellular_potential[e_idx] # convert to uV
+ ax_ep.plot(cell.tvec, sig, c=e_clr)
+ ax_morph.plot(electrode.elec_x[e_idx], electrode.elec_z[e_idx], marker='o', c=e_clr)
+
+ fig.savefig('example_nrn_EP.png')
+
+ plt.close(fig)
+
+if __name__ == '__main__':
+ cell = Cell()
+ cell.simulate()
+
+ num_elecs = 2
+ elec_x = 30 * np.ones(num_elecs)
+ elec_y = np.zeros(num_elecs)
+ elec_z = np.linspace(0, 100, num_elecs)
+
+ electrode = ExtElectrode(elec_x, elec_y, elec_z)
+ electrode.calc_extracellular_potential(cell)
+
+ plot_results(cell, electrode)
+
diff --git a/example/lfp/plot-lfp.py b/example/lfp/plot-lfp.py
new file mode 100755
index 0000000000000000000000000000000000000000..14ddde3eacd33da10de67a8879498c2d21eb6001
--- /dev/null
+++ b/example/lfp/plot-lfp.py
@@ -0,0 +1,71 @@
+#!/usr/bin/env python
+
+import argparse
+import matplotlib.pyplot as plt
+import json
+import math
+import sys
+
+# Read JSON output of lfp example from stdin and plot.
+#
+# The JSON data for timeseries is structured as:
+# <string>: { "unit": <string>, "time": [ <number>... ], "value": [ <number>... ] }
+# or:
+# <string>: { "unit": <string>, "time": [ <number>... ], "values": [[ <number>... ] ...] }
+#
+# 2-d morphology data is represented as samples (x, z, r) where r is the radius, one array
+# per branch. Extra point data (probe location, electrode sites) as pairs (x, z):
+# "morphology": { "unit": <string>, "samples": [[[<number> <number> <number>] ...] ...]
+# "probe": [<number> <number>], "electrodes": [[<number> <number>] ...] }
+
+def subplot_timeseries(fig, index, jdict, key):
+ data = jdict[key]
+ sub = fig.add_subplot(index, ylabel=data['unit'], title=key, xlabel='time (ms)')
+ ts = data['time']
+ vss = data['values'] if 'values' in data else [data['value']]
+
+ for vs in vss: sub.plot(ts, vs)
+
+def subplot_morphology(fig, index, jdict, key, xlim, ylim):
+ data = jdict[key]
+ unit = data['unit']
+ sub = fig.add_subplot(index, xlabel='x ('+unit+')', ylabel='y ('+unit+')', title=key, xlim=xlim, ylim=ylim)
+
+ for samples in data['samples']:
+ polys = [([x0-s0*dy, x0+s0*dy, x1+s1*dy, x1-s1*dy], [y0+s0*dx, y0-s0*dx, y1-s1*dx, y1+s1*dx])
+ for ((x0, y0, r0), (x1, y1, r1)) in zip(samples, samples[1:])
+ for dx, dy in [(x1-x0, y1-y0)]
+ for d in [math.sqrt(dx*dx+dy*dy)]
+ if d>0
+ for s0, s1 in [(r0/d, r1/d)]]
+
+ for xs, ys in polys: sub.fill(xs, ys, 'k')
+ sub.plot(*[u for x, y in data['electrodes'] for u in [[x], [y], 'o']])
+ sub.plot(*[u for x, y in [data['probe']] for u in [[x], [y], 'r*']])
+
+P = argparse.ArgumentParser(description='Plot results of LFP demo.')
+P.add_argument(
+ 'input', metavar='FILE', nargs='?', type=argparse.FileType('r'), default=sys.stdin,
+ help='LFP example output in JSON')
+P.add_argument(
+ '-o', '--output', metavar='FILE', dest='outfile',
+ help='save plot to file FILE')
+
+args = P.parse_args()
+j = json.load(args.input)
+
+fig = plt.figure(figsize=(9, 5))
+fig.subplots_adjust(wspace=0.6, hspace=0.9)
+
+subplot_morphology(fig, 131, j, 'morphology', xlim=[-100, 100], ylim=[-100, 600])
+subplot_timeseries(fig, 332, j, 'synaptic current')
+subplot_timeseries(fig, 335, j, 'membrane potential')
+subplot_timeseries(fig, 338, j, 'ionic current density')
+subplot_timeseries(fig, 133, j, 'extracellular potential')
+
+if args.outfile:
+ fig.savefig(args.outfile)
+else:
+ plt.show()
+
+plt.close(fig)