www.neuron.yale.edu

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()
*/

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()