Python & NEURON
Interfacing NEURON with the Python API
Table of Contents
1. Getting Started:
1.1. Installing NEURON
To get started first one must install NEURON. To do so please follow the instructions based on your operating system 1.1.1, 1.1.2, 1.1.3.
1.1.1. Windows
First download and install NEURON from the MSWin Setup.exe:
Setup.exe, or can be downloaded from Yale's website at: https://www.neuron.yale.edu/neuron/download#mswin
After having done so move onto the Python Section
1.1.2. MacOS
For Intel Macs use the following link:
For M1/ARM Macs use this:
Both of the above links may be downloaded directly from Yale's page at https://www.neuron.yale.edu/neuron/download#macos
Now carry on to the Python Section
1.1.3. Linux
On Linux in your python environment of choice simply run
pip install neuron
Then you may easily carry on to the Python Section.
1.2. Setting up Python
Once you have NEURON installed, using it is as simple as:
from neuron import h,gui
Where h is the interpreter for hoc files, and gui ensures that
NEURONs graphical interface is started, allowing for quick model
building, viewing, parameter changes etc. On MacOS that interface
looks like the following:
When working with models like those from ModelDB one will often have to perform some extra work prior to importing the downloaded model.
Depending on your platform in the commandline you will have to type one of the following:
Windows:
nrnivmodl
MacOS:
mknrndll || nrnivmodl
Linux:
nrnivmodl
After which a file will be generated in a folder labelled either
x86_64 or arm64, entitled nrnmech.dll on windows,
libnrnmech.dylib on MacOS, or libnrnmech.so. After this file is
generated from python you must load the file prior to the model.
2. A Working Example
As an example we will load and use the model for a STN neuron, from
ModelDB access 151460. After having downloaded and extracted the zip file, and having
generated the nrnmech file as in the previous paragraph you can load
and view the STN neuron as follows (as done on MacOS, replace file
paths as necessary):
from neuron import h,gui h.nrn_load_dll("./Miocinovic/MiocinovicEtAl2006/arm64/libnrnmech.dylib") h.load_file("./Miocinovic/MiocinovicEtAl2006/STN.hoc") h.Shape()
We shall work from a modified version of this model Download from: HERE
Having run the above code you should now have a 3D view of the modelled STN neuron.
To access the properties of the loaded neuron, we can access them from
the hoc interpreter h, simply by using the dot operator. For
example:
h.dend[0] # Grabs the first dendrite compartment for s in h.soma: # Iterates through all soma print(s) print(h.somaelements) # Prints the number of soma elements
A problem with loading models this way (with h.load_file(...)) is
that unless the model was constructed with a template all the values
are loaded in as global variables, meaning if we loaded a different
neuron, or even tried to make another of the same neuron, they would
overwrite one another.
To fix this is somewhat easy, though tedious. We must encapsulate the
.hoc model of choice with begintemplate NAME and endtemplate
NAME, then we must declare all of the variables used in the model as
public. Then we must add a function/procedure named init that
executes all the previous code in order; doing so may require
encapsulating regions of the previous code in other function/procedure
blocks for the STN neuron model that looks like the following:
begintemplate STN
public total,dendelements,somaelements,iselements,\
axonnodes,paranodes1,paranodes2,axoninter,\
celsius,v_init,fiberD,paralength1,paralength2,\
nodelength,space_p1,space_p2,space_i,rhoa,mycm,\
mygm,axonD,nodeD,paraD1,paraD2,deltax,nl,Rpn0,Rpn1,\
Rpn2,Rpx,interlength,s
.
.
.
proc composeCell() {
connect soma[1](0), soma[0](1)
connect soma[2](0), soma[1](1)
.
.
.
node[27]{
pt3dadd(-6745.186944,2709.982462,-2317.434852,1.400000)
pt3dadd(-6746.142167,2709.686826,-2317.447023,1.400000)
}
}
proc init() {
model_globals()
dependent_var()
initcell()
composeCell()
}
endtemplate STN
After encapsulating the model in a template, once loaded into Python
with the load_file mechanism we may create as many instances of the
neuron type as we would like by calling neuron = h.NAME(), the
assignment is necessary as otherwise it is deleted by NEURON. In the
case of the STN neuron example an STN neuron may be created by calling
stn_neuron = h.STN().
The variable stn_neuron would then have all the properties defined
in the public statement in our template. As such we can acess the
soma stn_neuron.soma[?] where there are stn_neuron.somaelements
total elements, and similar for dendrites (.dend[?]), axon nodes
(.node[?]) etc.
2.1. Stimulation
An important feature of simulations is the ability to apply some kind of stimulation, be it either a constant voltage, current or otherwise. Thankfully with our newly created templates for neurons we have easy access to all the components necessary to start stimulation.
The most basic form of interaction would be the current & voltage clamps, according to the documentation they are of the form:
from neuron import h, gui h.nrn_load_dll("./Miocinovic/MiocinovicEtAl2006/arm64/libnrnmech.dylib") h.load_file("./Miocinovic/MiocinovicEtAl2006/STN.hoc") # Modified to have # begin and end # templates etc. neuron = h.STN() stimulator = h.IClamp(neuron.soma[0](0.5)) stimulator.delay = 10 # ms stimulator.dur = 0.1 # ms stimulator.amp = 15 # nA
Teh above code sample creates an STN neuron and a current clamp that
is attached to the midpoint (0.5) of the first [0] somatic
segment. Where the stimulator will activate with a 1ms delay for
0.1ms with a current of 2nA.
2.2. Recording
Now that we have access to current stimulation you may be wanting to recording the membrane potential to watch the action potential develop. Thankfully this is rather easy:
... node_v = h.Vector().record(stn.node[0](0.5)._ref_v) from neuron.units import ms,mV h.finitialize(-77 * mV) h.continuerun(500 * ms) import matplotlib.pyplot as plt plt.plot(node_v) plt.show()
Where h.Vector() creates a vector object of arbitrary length
.record records a point voltage of a given source throughout the
simulation.
Combining the IClamp with the recording code above should result in the following graph:
2.3. Connecting Neurons
Of course to construct a network between neurons we need to be able to communicate between them. In NEURON this requires two things:
- A synapse constructed on a given segment
- A NetCon from a Segment to a Synapse
2.3.1. Creating a Synapse
Creating a synapse is thankfully rather easy, once you've selected your target section, for example:
target = stn_neuron.dend[0]
You can create a synapse with one of the following functions:
- ExpSyn
Creates a simple synapse where the current is maximal upon the delay and decays with an exponent of a given \(\tau\) - Exp2Syn
A synapse whose current is constructed between two exponentials, one for the rise in current, another for the fall - AlphaSynapse
Creates a current defined by the alpha function: \(i = g*(V-E)\) where g is 0 prior to onset and \(g=g_{\text{max}}\frac{t-\text{onset}}{\tau} e^{-\frac{t-\text{onset}-\tau}{\tau}}\) after onset.
As an example we can create a ExpSyn synapse on our target with
the following code:
synapse = h.ExpSyn(target(0.5))
2.3.2. Connecting to a Synapse
Once you have a synapse to connect to, denoting what the presynaptic
element is is rather easy.
For example let us consider an axon node that we want to use as the
presynaptic element to out previous target. We may construt a
connection with the following:
node = stn_neuron.node[0] connection = h.NetCon(node(0.5)._ref_v,synapse,sec=node) # When connecting segments, as compared to point # processes, you're required to specify "sec" connection.weight[0] = 1
When creating a connection to a synapse you are required to specify the weight of said synapse, as otherwise it defaults to 0 and you will not be able to observe the connection in action.
2.4. Putting it all Together
Now that we have all the necessary components for constructing a synapse, injecting current, and connecting neurons we can put together a simple model wherein we create a neuron that synapses upon itself.
from neuron import h, gui from neuron.units import ms,mV h.nrn_load_dll("./Miocinovic/MiocinovicEtAl2006/arm64/libnrnmech.dylib") h.load_file("./STN_Template.hoc") # Modified to have # begin and end # templates etc. stn_neuron = h.STN() stimulator = h.IClamp(stn_neuron.soma[0](0.5)) stimulator.delay = 10 # ms stimulator.dur = 0.1 # ms stimulator.amp = 15 # nA target = stn_neuron.dend[0] synapse = h.ExpSyn(target(0.5)) synapse.tau = 1.5 node = stn_neuron.node[0] connection = h.NetCon(node(0.5)._ref_v,synapse,sec=node) # When connecting segments, as compared to point # processes, you're required to specify "sec" connection.weight[0] = 50 connection.delay = 100 * ms node_v = h.Vector().record(stn_neuron.node[0](0.5)._ref_v) h.finitialize(-77 * mV) h.continuerun(500 * ms) import matplotlib.pyplot as plt plt.plot(node_v) plt.show()
Doing so should produce the following plot:
3. Running HOC Code directly
Once you have import h from neuron you can run any arbitrary hoc
code through the TopLevelInterpreter h like the following:
from neuron import h h(""" NUM_THING = 5 neuron = STN() syn = ExpSyn(neuron.dend[0](0.5)) """)
This allows you to interact with scripts from python rather easily,
allowing you to declare and modify variables for use in the
imported/loaded .hoc files.
4. Simulating DBS
With regards to how one might simulate DBS, thankfully NEURON
incporates a prcess called extracellular that many models use,
including our STN neuron model. This process incorporates the dynamics
surrounding the extracellular space of a neuron, allowing us to change
its field potential and have that quickly and easily incoporated into
the dynamics of the membrane.
To simulate a DBS pulse using this extracellular process we can do
the following:
DBS_Current = 3 # mA sigma = 0.3 # S/m neuron = STN() for s in neuron.s: s.sec.e_extracellular = (DBS_Current/(4*Math.pi*sigma/1e6))/DISTANCE # DISTANCE is calculated distance from electrode
5. Approximating LFP
To approximate a LFP recording is rather simple, albeit tedious for to implement when using many neuron models. All that is required is to calculate the distance of the neuronal segment to the electrode weighted by distance and summate that for every segment in every neuron.
An example with the STN Neuron may be done as follows:
neuron = STN() field_rec = 0 for s in neuron.s: field_rec = field_rec + s.sec.i_membrane/DISTANCE # DISTANCE to Electrode
Doing so for every element allows you to construct a basic LFP.
Further processing is necessary to make the signal match more closely to physiological recordings.
One of the easiest things to do is to excise recordings from neurons that are within 500 \(\mu m\) of the electrode, such that if any segment of the model is within that range, do not account for it when summing the LFP.