Add LFP example/demo. (#1073)

* Add C++ analogue to simple NEURON demo provided by LFPy project under example/lfp.
* Add plot script for rending example output to match that from NEURON demo.

example/CMakeLists.txt

@@ -10,3 +10,4 @@ add_subdirectory(ring)
 add_subdirectory(gap_junctions)
 add_subdirectory(single)
 add_subdirectory(probe-demo)
+add_subdirectory(lfp)

example/lfp/CMakeLists.txt

+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}")

example/lfp/README.md

+# 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
+```
+

example/lfp/lfp.cpp

+#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";
+}

example/lfp/neuron_lfp_example.py

+#!/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)
+

example/lfp/plot-lfp.py

+#!/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)