Merge pull request #4 from subhacom/master

Added test for KC model and updated README for NEURON

NEURON/README.md

@@ -1,17 +1,50 @@
 This directory contains scripts for simulating the GGN model and the
-mushroom body olfactory circuit around it in NEURON 7.4 with Python
-2.7 related to the article:
+mushroom body olfactory circuit around it (originally in NEURON 7.4 with Python
+2.7) related to the article:
 
 "Feedback inhibition and its control in an insect olfactory circuit".
 Subhasis Ray, Zane Aldworth, and Mark Stopfer; 2020. eLife.
 
+- April 2024: Updated to run with NEURON 8.2.4 on Python 3.12 during COMBINE/Harmony-2024 Hackathon.
 
 # Running the simulations
-  For testing, run simulation scripts from the current directory. Edit the nrn/nrninit.sh (or nrninit.bat on Windows) with the path for this directory. 
+  Required packages: 
+
+  * networkx (graph data structures)
+  * neuron (simulation)
+  * sklearn (spatial clustering based connections)
+  * pint (units)
+  * pyyaml (configuration)
+  * nsdf (saving data): `pip install git+https://github.com/nsdf/nsdf.git`
+  * vtk (for 3D visualization)
+  * numpy (numerics)
+  * matplotlib (visualization)
+
+  For testing, run simulation scripts from the current directory
+  (where this `README.md` file resides). Edit the `nrn/nrninit.sh`
+  (or `nrninit.bat` on Windows) with the path for this directory.
+
+  If you are running from another directory, edit `nrn/nrninit.sh` and
+  set `MODEL_DIR` to the path of that directory.
+  * To initialize the environment variables run the `nrn/nrninit.sh`
+    file. For Linux/OSX with bash shell: `source nrn/nrninit.sh`
   * You also have to build the NEURON mechanisms (.mod) files:
-  `nrnivmodl mb/mod`
-  * To simulate GGN in isolation to study signal attenuation with distance, run:
-    `mb/test_cell/ggn_voltage_attenuation_vclamp.py`
+
+    `nrnivmodl mb/mod`
+
+  * To simulate GGN in isolation to study signal attenuation with
+    distance, run:
+
+      `python mb/test_cell/ggn_voltage_attenuation_vclamp.py`
+
+  * To simulate the PN->KC<->GGN network model, create a directory
+    called `data` for dumping simulated data and run:
+
+    `python mb/network/pn_kc_ggn_network.py`
+
+    You may want to edit the configuration file
+    `mb/network/config.yaml` to reduce the number of cells in the
+    network and the duration of the simulation.
 
 
 # analysis
@@ -62,7 +95,9 @@ Subhasis Ray, Zane Aldworth, and Mark Stopfer; 2020. eLife.
 
 
 ## slurm : utility scripts for running simulations in batch mode under slurm (on NIH biowulf).
-
+    *NOTE: Many of these scripts have hard coded absolute paths which
+    need to be updated for running in a new environment.*
+
 -   `batch_run_remove_kcs_run.sh` : example script for running successive
     simulations after removing high spiking KCs.
 -   `circular_run` : scripts for running running successive simulations
@@ -116,11 +151,13 @@ Subhasis Ray, Zane Aldworth, and Mark Stopfer; 2020. eLife.
 
 # nrn
 
--   `nrnutils.py` : Utilities for handling NEURON model
+-   `nrnutils.py` : utilities for handling NEURON model
     -   convert a NEURON cell model into a networkx graph
     -   insert alpha synapses
     -   insert ion channel mechanisms
     -   set up recording of Vm
+- `display_celltemplate.py`: 3D visualization of the morphology of a
+  neuron specified in NEURON cell template
 -   `localized_input_output.py` : apply synaptic inputs at specified
     branches.  This scripts runs simulation with synchronous synaptic
     inputs at multiple compartments on specific branches.

NEURON/mb/test_cell/test_kc.py

@@ -0,0 +1,69 @@
+# test_kc.py --- 
+# 
+# Filename: test_kc.py
+# Description: 
+# Author: Subhasis Ray
+# Created: Wed Apr 10 16:56:23 2024 (+0530)
+# Last-Updated: Wed Apr 10 17:54:28 2024 (+0530)
+#           By: Subhasis Ray
+# 
+
+# Code:
+
+"""Create model from cell template file `filename`. The cell name in
+    the template should be specified in `cellname`.
+
+"""
+from config import Q_, ur, h
+import numpy as np
+from matplotlib import pyplot as plt
+import ephys
+
+# filename and cellname specify cell template file and the name of the
+# cell in the template respectively
+filename='mb/cell_templates/kc_1_comp.hoc'
+cellname='KC'
+
+
+Em = Q_('-70mV')        # Wustenberg, et al. 2004, TABLE 2 legend
+
+h.xopen(filename)
+kc = eval(f'h.{cellname}()')
+delay = Q_(100, 'ms')
+duration=Q_(500.0, 'ms')
+inj_current = Q_(16, 'pA')
+clamp = ephys.setup_current_clamp(kc.soma, delay=delay, duration=duration, amplitude=inj_current)
+
+# Recording data
+vvec = ephys.setup_sec_rec(kc.soma, 'v')[0]
+ik_vec = ephys.setup_sec_rec(kc.soma, 'ik')[0]
+ina_vec = ephys.setup_sec_rec(kc.soma, 'ina')[0]
+im_vec = h.Vector()
+im_vec.record(clamp._ref_i)
+tvec = h.Vector()
+tvec.record(h._ref_t)
+
+h.tstop = delay.to('ms').m + duration.to('ms').m
+h.v_init = Em.to('mV').m
+h.init()
+h.run()
+# Now plot the data
+t = Q_(np.asarray(tvec.x), 'ms')
+fig, axes = plt.subplots(nrows=2, sharex='all')
+axes[0].plot(t, np.array(vvec.x))
+axes[1].plot(t, np.array(im_vec.x))
+axes[0].set_ylabel('Vm (mV)')
+axes[1].set_ylabel('Im (nA)')
+
+data = np.array([t.to('s').m, np.array(vvec.x)], dtype=[('t', float), ('v', float)])
+data_file = 'kc_vm.npy'
+np.save(data_file, data)
+
+print(f'Simulated {cellname} from template file {filename}.')
+print(f'Settling time {delay} followed by {inj_current} current injection for {duration}')
+print(f'Saved Vm in {data_file}')
+plt.show()
+
+
+# 
+# test_kc.py ends here

NEURON/morphutils/neurograph.py

@@ -453,6 +453,26 @@ def remove_shorter_edges(G, n0=1, lim=0.1, verbose=False):
                     todo[n] = None
 
 
+def to_undirected_nrn(G):
+    """Workaround to convert graph with nodes whose attribute 'orig'
+    is neuron section"""
+    G_ = G.__class__()
+    G_.add_nodes_from(G)
+    for edge in G.edges():
+        print('#', edge)
+        G_.add_edge(*edge)
+    for node in G.nodes:
+        G_.nodes[node].update(G.nodes[node])
+        if G.nodes[node]['orig'] is not None:
+            G_.nodes[node]['orig'] = G.nodes[node]['orig'].name()
+        else:
+            G_.nodes[node]['orig'] = None
+    G2 = G_.to_undirected()
+    for node in G.nodes:
+        G2.nodes[node]['orig'] = G.nodes[node]['orig']
+    return G2
+
+
 def remove_longer_edges(G, lim=100.0, verbose=False):
     """Delete the edges which are longer than lim. Returns
 
@@ -476,10 +496,18 @@ def remove_longer_edges(G, lim=100.0, verbose=False):
             # edges are already sorted in descending order, skip shorter edges
             break
     components = []
+    # NEURON7.7/NetworkX 3.x sim incompatible - sections cannot be deepcopied as they cannot be pickled
+    G_ = G.__class__()
+    G_.add_nodes_from(G)
+    G_.add_edges_from(G)
+    for node in G.nodes:
+        G_.nodes[node]['orig'] = G.nodes[node]['orig'].name()
     if len(long_edges) > 0:
-        G2 = G.to_undirected()
+        G2 = G_.to_undirected()
         for (n0, n1) in long_edges:
             G2.remove_edge(n0, n1)
+        for node in G2:
+            G2.nodes[node]['orig'] = G.nodes[node]['orig']
         for sub in nx.connected_component_subgraphs(G2):
             components.append(nx.DiGraph(G.subgraph(sub.nodes())))
     else:
@@ -520,7 +548,7 @@ def renumber_nodes(G, start=1, undirected=True):
     ret = nx.DiGraph()
     node_map = {}
     if undirected:
-        G = G.to_undirected()
+        G = to_undirected_nrn(G)
     # The nodes are connected as child->parent, we are starting from root,
     # hence reverse    
     for ii, (n1, n2) in enumerate(nx.dfs_edges(G, start)):