Help for adding multiple neurons and connecting them.

[Csmith3WPI]

Hello, I've recently been working on a project to simulate how signals pass from neuron to neuron and we'd like to use some pre-collected data. We were able to use a single SWC file but when we attempt to use more I'm unable to properly connect/add all of the neurons. Not sure if I've explained that correctly. Anyhow, ill share the code and see if yall have any ideas. Much appreciated in advance!

Code: Select all

// Modified by ~
// Load topology from six different .swc files and connect them

objref cell[6], swc[6], i3d[6], null, dend_indices, axon_indices, dend_list[6], axon_list[6]

// File names for SWC morphologies
objref filenames[6]
filenames[0] = new String("864691135356398031_0.swc")
filenames[1] = new String("864691135404857070_1.swc")
filenames[2] = new String("864691136389878647_2.swc")
filenames[3] = new String("864691135777479597_3.swc")
filenames[4] = new String("864691135119056861_4.swc")
filenames[5] = new String("864691135577012766_5.swc")

// Procedure to import SWC files into NEURON
proc import_all_swc() {
    for i = 0, 5 {
        swc[i] = new Import3d_SWC_read()
        swc[i].input(filenames[i].s)
        i3d[i] = new Import3d_GUI(swc[i], 0)
        cell[i] = new SectionList() // dummy init to avoid nil error
        i3d[i].instantiate(null)

// Initialize lists
        dend_list[i] = new List()
        axon_list[i] = new List()

        // Classify sections
        forsec cell[i] {
            if (strm(secname(), "1")) {
                dend_list[i].append(new SectionRef())
            } else if (strm(secname(), "2")) {
                axon_list[i].append(new SectionRef())
            }
        }
    }
}

import_all_swc()
/*
// Set up synapses on the postsynaptic cell (cell[0])
objref exc_syn[5], nc_exc[5]

// Define synapse parameters and connections
proc connect_synapses() {
    // dendrite indices and axon segment indices for each presynaptic neuron
    dend_indices = new Vector(5)
    axon_indices = new Vector(5)
    
    dend_indices.x[0] = 12
    dend_indices.x[1] = 27
    dend_indices.x[2] = 44
    dend_indices.x[3] = 49
    dend_indices.x[4] = 4

    axon_indices.x[0] = 43
    axon_indices.x[1] = 62
    axon_indices.x[2] = 39
    axon_indices.x[3] = 65
    axon_indices.x[4] = 55

    for i = 0, 4 {
        dend_list[0].o(dend_indices.x[i]).sec exc_syn[i] = new ExpSyn(0.5)
        exc_syn[i].tau = 2
        exc_syn[i].e = 0

        axon_list[i+1].o(axon_indices.x[i]).sec nc_exc[i] = new NetCon(&v(0.5), exc_syn[i])
        nc_exc[i].threshold = -20
        nc_exc[i].weight = 0.05
    }
}

connect_synapses()
*/

[ted]

Nice clean code. Looks like somebody has experience with C or Matlab.

General suggestions when developing network models: First define and test each cell class all by itself. Then try creating and connecting just two cells.

Cell classes are important for several reasons. First is avoiding collisions in name space, second is code readability. Every section created by top-level code (statements that are not inside a class definition) exists in the same name space. How many cell instances will your model net have? Will each one have a section called soma, dend, axon? How are you going to prevent each
create soma
statement from stepping on a pre-existing section called soma? I bet that's happening with your otherwise very nice example. I'm not sure if the hoc interpreter will generate error messages every time there's a name collision, so you might have to run a test with just two morphology files and then check the topology() printout or a ShapePlot.

"Oh, I'll prevent destructive interactions by using unique strings. Or maybe I'll use arrays of names so each soma will be called soma[j][k][n] where j is the index of the cell type, k is the index of the cell instance, and n is the index of however many pieces the soma is built from (depending on the whim or judgement of whoever did the morphological measurements)." Ow. That's going to be hard to write, and hard to read/debug/maintain.

Far better to define cell classes that have meaningful names e.g. L5pyr or Renshaw.

You can do this with hoc or Python. Python is best for new code development for reasons that include (but are not limited to): its expressive power, the availability of useful 3rd party libraries, and the fact that so many people (including undergrads and grad students) already know at least some Python.

I'll post a Python example from the NEURON course we ran at ASU in Tempe in March 2025.

[ted]

More code than I like to paste as a single lump in one of these threads.

Example of Python code that creates a cell class with morphology read from an SWC file. This was developed by Robert McDougal as one possible answer to one of the in-class hands-on exercises. It uses matplotlib to display simulation results. With a few changes, the cell class could be made suitable for use with NEURON's built-in Interviews-based GUI (which would be much more interactive than use of matplotlib).

Code: Select all

import requests
from neuron import h
from neuron.units import mV, ms, µm
import pathlib
import sys
import math
import matplotlib.pyplot as plt
h.load_file("import3d.hoc")
h.load_file("stdrun.hoc")

if not pathlib.Path("c91662.CNG.swc").exists():
    print("Morphology file c91662.CNG.swc not found.")
    print("Download morphology from:")
    print("    https://neuromorpho.org/dableFiles/amaral/CNG%20version/c91662.CNG.swc")
    print("and place it in the current directory.")
    sys.exit(-1)

class Pyramidal:
    def __init__(self):
        # morphology
        cell = h.Import3d_SWC_read()
        cell.input("c91662.CNG.swc")
        i3d = h.Import3d_GUI(cell, False)
        i3d.instantiate(self)

        # biophysics
        h.hh.insert(self.soma)
        h.hh.insert(self.axon)
        # NOTE: hh includes a leak mechanism, so adding pas to
        #       the soma and axon is redundant, but it was the
        #       assignment
        h.pas.insert(self.all)
        for sec in self.all:
            for seg in sec:
                seg.pas.g = 1e-6
        
        # discretize
        self._apply_d_lambda()
        
    def _apply_d_lambda(self, d_lambda=0.1, freq=100):
        # NOTE: if you're not using import3d, this requires stdlib.hoc
        for sec in self.all:
            sec.nseg = 1 + 2 * math.ceil(
                (sec.L / (d_lambda * h.lambda_f(freq, sec=sec))) / 2)

p = Pyramidal()

# drive the cell with a current clamp
stim = h.IClamp(p.soma[0](0.5))
stim.delay = 5 * ms
stim.dur = 0.2 * ms
stim.amp = 10  # nA

# setup recording (not strictly necessary for the assignment)
v = h.Vector().record(p.soma[0](0.5)._ref_v)
t = h.Vector().record(h._ref_t)

# initialize and run
h.finitialize(-65 * mV)
h.continuerun(6 * ms)

# plot the distribution of the membrane potential
fig = plt.figure()
ps = h.PlotShape(False)
ps.variable("v")
ps.plot(fig)
fig.suptitle(f"Membrane potential distribution at t = {h.t} ms")

# continue the simulation
h.continuerun(20 * ms)

# plot the membrane potential at the center of the soma over time
plt.figure()
plt.plot(t, v)
plt.title("Membrane potential at the center of the soma")
plt.xlabel("Time (ms)")
plt.ylabel("Membrane potential (mV)")

plt.show()

[ted]

Final comments: have you verified that your morphology files are suitable for computational modeling with NEURON? Are there orphaned branches, bottlenecks or other oddities affecting neurite diameters (e.g. minimum diameter larger than 1 um, "favorite numbers")? Is there evidence of z axis backlash (xz or yz view of cell shows abrupt jumps parallel to the x axis)?