Simulating Biological Neural Networks with NetPyNE
A hands-on tutorial: Exploring the effect of neuromodulation on oscillatory patterns
This tutorial provides a glimpse into simulating and analyzing biological neural networks. Following the two previous posts (prepare fMRI data for analysis, and using LSTMs to analyze rs-fMRI) — this post takes a different approach. Where the two previous posts focused on investigating the bigger picture, analyzing imaging data, this post introduces the benefits of bottom-up inspections of neural dynamics. Neural simulations enable altering parameters at different levels of analysis, inspecting their effects on the dynamics of the network as a whole.
In this post, I present a brief overview of NEURON & NetPyNE and provide a docker image including the full installation of both software. Docker is a customized virtual environment that delivers an isolated, ready-to-run deployment model (read more about containers here). That means, that there is nothing else needed to be installed. Simple.
Second, I provide a well-explained example demonstrating one type of analysis NetPyNE offers. This example shows how modifying different neuromodulators impact the network’s oscillatory pattern. You can find the full code here.
Modeling biological neural networks
Biological neural networks commonly discuss neural dynamics in three; the macroscopic level (whole brain and regional dynamics), the mesoscopic level (population dynamics), and microscopic (cellular and biophysical dynamics). Elements of each layer mutually impact the dynamics at different scales. In the example below, neuromodulators at the biophysical level effect firing dynamics at the population level. In NEURON, one can configure the model using parameters at these levels. The example below demonstrates this.
NEURON and NetPyNE
The NEURON simulator provides a robust environment for neural models at different levels; it supports the biological properties of neurons and synaptic mechanisms, as well as neural architecture. Nicely, it abstracts most of the mathematical representations and offers the user to think in familiar neuroscience concepts.
However, NEURON is developed in hoc, which is (lightly said) a less familiar development environment. Thankfully, NetPyNE came to save the day. NetPyNE is a python package that facilitates developing simulation and analysis of biological neuronal networks using the NEURON simulator.
Installation
NEURON (and NetPyNE) requires a non-trivial installation process that may be a bit challenging. To help you skip this pain and delve right into modeling and analyzing, I have created a couple of ready-to-go Docker images. The first offers NetPyNE over jupyter-notebook and the second NetPyNE with python3. The first is suitable for simple models, exploring objects native to NetPyNE using plots, whereas the latter could support more extensive simulations.
Both Dockerfiles support multiprocessing via mpich. However, I would be careful running a big or long simulation locally (as the complexity, thus the running time and the memory required, rises super-linearly).
Exploring neuromodulation and oscillatory patterns
Theory (why is it even an exciting example?)
Briefly, neural oscillations (at a macroscopic level) reflect the synchronization of a group of neurons. That means that oscillatory activity commonly results from harmonized firing patterns. An abundance of evidence has linked GABA receptor regulation to oscillatory patterns (Kujala et al., 2015; Gonzalez-Burgos et al., 2008; Hentschke et al., 2009; Rovó et al., 2014; Lowet et al., 2016). Additionally, a few recent papers recognize the importance NMDA receptors have on synchronized activity (Molina et al., 2014; Placantonakis et al., 2001; Hunt et al., 2013) and some have associated such relation to the NMDA hypofunction hypothesis of schizophrenia (Jadi et al., 2016). Impaired oscillatory patterns are associated with various memory and perceptual deficits (Vishnoi et al. ,2015).
Simulation (So how to do it?)
I chose to demonstrate using the notebook environment so I could easily show the figures. To run the docker’s environment, follow the instructions here.
In this section, we firstly will build a model (configuring the network level, the cellular level, and the biophysical level). Second, we run it comparing three different neuromodulator configurations and inspecting the changes shown at the network firing rate.
Once you have the notebook running in the Docker environment, let’s import NetPyNE.
NetPyNE is organized via a complex set of dictionary-like configuration objects. The main two are the network parameters and the simulation parameters.
The network below describes 200 pyramidal cells of the Hodgkin–Huxley model. The Hodgkin–Huxley (HH) is a common, well-established model, which is native to NetPyNE (i.e., does not require any additional files).
Next, we delve into the cellular level and modify (homogenously) the cells themselves. Creating a CellRule dictionary, we focus on the cell’s soma. Notice that each cell model will have a slightly different set of parameters to set.
- The soma’s geometry is defined by its diameter (
diam), its length (L), and its axial resistivity (Ra). See more here. - The soma’s mechanisms (unique to the HH) describe the cell’s conductance. Its maximum sodium channel conductance is marked by
gnbar(g for conductance, Na denotes sodium, and bar reflects the max value), similarly, the cell’s maximum potassium conductance is defined by thegkbarparameter. Additionally, the cell’s leakage conductance is determined bygl, and its leakage reversal potential is.stated byel. Lastly, we set the cell’s initial potential to -71 mv.
After defining the population’s features on the cellular and the morphological levels, we move to the biophysical level and set the synaptic mechanisms. The synaptic mechanism dictates the electric relay between neurons. A simple model (which is, as well, native to NetPyNE) is Exp2Syn standing for a double exponential synaptic mechanism. This model has two kinetic states; the synapse exponential rise time is given by tau1, and its exponential decay time constant by tau2 . The parameter e denotates the synapse reversal potential. See more here.
The model comprises three types of synaptic mechanisms; GABA receptors, AMPA receptors, and NMDA receptors. NMDA receptors and AMPA receptors are both excitatory glutamatergic channels whereas GABA receptors are inhibitory.
**Pay attention that using the Exp2Syn model for all three reflects an aggressive simplification (find an NMDA specific model here). However, for this tutorial- such abstraction is tolerated.
The last step in setting up the network is defining its connectivity. Using the synaptic mechanisms described above, this section connects neurons within the population. Simply, preConds refers to the presynaptic neuron and postConds to the postsynaptic one. The weight of the connection indicates its strength. The delay, here (illustrated by a normal distribution), dictates the time delta between the excitation of the presynaptic and postsynaptic neurons. The synMech expresses the synaptic mechanism used.
This network, as defined above, is illustrated in the following figure, where the blue nodes reflect neurons, and the red edges illustrate their synaptic connections (of any kind). This plot, offered by NetPyNE, is a bit over-simplistic and messy but it gives a good intuition to the spatial organization of the network. Remember that not all synaptic connections are of the same kind or of the same strength (despite them all looking exactly the same in the plot).
An external stimulus is introduced to excite the network. The stimulus towards the population demonstrates an electric current of rate 30 Hz and starts after 1 ms. Its weight and delay regulate its magnitude.
To ‘bring the network to life,’ the configuration parameters of the simulations are set. Here you can find the duration of the simulation (duration ), the timestep length(dt), the randomization seeds(seeds), reporting verbose (verbose), and the initial potential (hparam[‘v_init’]). Additionally, the recording configuration is set. Here, all cells are recorded (by leaving recordCells empty), the voltage is being captured from the middle of the soma (set in recordTraces[‘Vsoma’]), and the recording timestep is 0.1 ms (recordStep). Setting the parameter recordStim to True provides an explicit confirmation to record the simulation based on the set configuration.
Finally, let’s run the simulation! This step creates the network using the network parameters netParams, runs the simulation and records it according to the simConfig parameters.
The console will report as it builds the network and runs the simulation. Pay attention that the details match the configuration above. When NetPyNE is fully done, it reports the total running time.
Analysis (So what is it good for?)
Let’s see how neuromodulators alter the firing rate synchrony. First, inspecting the activity with the parameters seen above (i.e., modeling ‘typical’ activity).
The raster plot illustrates the spiking activity of the whole population for each neuron (y-axis) over time (x-axis). Thus, each blue dot reflects a spike of a cell in time. The red lines indicate the spiking synchrony (which is quantified to a percentage at the top of the figure under ‘sync’). The other quantities reported are the number of cells, the synapses per cell (NOT the number of connections per cell as not every synapse is active), and the network rate in Hz.
Here, the network seems to spike relatively synchronously. If to be precise — 89% synchronous.
The spike histogram sheds more light on each harmonized spike. By default, it counts the spikes for every 5 ms bin and reports the rate of spiking for each bin. This histogram plots the mean firing rate of the population (y-axis) over time (x-axis) and demonstrates the magnitude of the population activation. Here, we observe that at each 5 ms, the rate of spiking is approximately between 20 to 35 Hz.
- Pay attention, both plots use rate (Hz) in their report. The rate given by the raster plot indicates the spiking rate of the whole network whereas the rate reported by the spike histogram reflects the count of spikes at each bin.
OK, but — will the model capture the link between modifying neuromodulators and the network spiking behavior ??
First, let us ‘introduce’ an NMDA antagonist to the network by silencing the NMDA connections (i.e., reducing the connectivity weight of the NMDA to 0).
After rerunning the simulation (using the same code as seem in runsim.py)- we see the following raster and spike histogram plots. The raster plot below indicates less synchronization (0.77 compared to 0.89 before). The synchronization shrinkage comes with an increased activity rate (7.5 Hz firing rate); the population fires more often. The spike histogram reveals a slight reduction in the firing rate. Considering the increase in the spikes’ frequency — such a reduction makes sense (since the refractoriness of each neuron remains the same, an increased frequency of the network’s spiking could result in less neurons spiking at each time). This result corresponds to the literature (Molina et al., 2014).
Lastly, let’s reset the NMDA to its original value and boost the GABA synapses. Based on literature — we would expect to see a much larger synchronization with a much smaller total activation (due to the increased inhibition in the network).
The plots below follow these expectations. We observe that the synchronization went up to 92% (!!), and the spikes lessened both in frequency (2.7Hz) and rate (15–30 Hz).
Summary
This tutorial achieves two main goals- the first is providing one with ready-to-go tools to develop and analyze using NEURON and NetPyNE. Secondly, it demonstrates a simple example of linking multiple layers of analysis via modeling of biological neural networks. Lastly — it offers a well-documents glimpse into NetPyNE, which hopefully helps to deal with mode complicated models.
Additional resources
The full documentation for using NEURON with python can be found here. A basic NetPyNE tutorial is offered on their website here.
An extensive database of models (in NEURON but also other languages) is freely available here. Similarly, a database of (ion) channels can be found here.
Refrence
Gonzalez-Burgos, G., & Lewis, D. A. (2008). GABA neurons and the mechanisms of network oscillations: implications for understanding cortical dysfunction in schizophrenia. Schizophrenia bulletin, 34(5), 944–961.
Hentschke, H., Benkwitz, C., Banks, M. I., Perkins, M. G., Homanics, G. E., & Pearce, R. A. (2009). Altered GABAA, slow inhibition and network oscillations in mice lacking the GABAA receptor β3 subunit. Journal of neurophysiology, 102(6), 3643–3655.
Hunt, M. J., & Kasicki, S. (2013). A systematic review of the effects of NMDA receptor antagonists on oscillatory activity recorded in vivo. Journal of Psychopharmacology, 27(11), 972–986.
Jadi, M. P., Behrens, M. M., & Sejnowski, T. J. (2016). Abnormal gamma oscillations in N-methyl-D-aspartate receptor hypofunction models of schizophrenia. Biological psychiatry, 79(9), 716–726.
Kujala, J., Jung, J., Bouvard, S., Lecaignard, F., Lothe, A., Bouet, R., … & Jerbi, K. (2015). Gamma oscillations in V1 are correlated with GABA A receptor density: A multi-modal MEG and Flumazenil-PET study. Scientific reports, 5, 16347.Chicago
Lowet, E., Roberts, M. J., Bonizzi, P., Karel, J., & De Weerd, P. (2016). Quantifying neural oscillatory synchronization: a comparison between spectral coherence and phase-locking value approaches. PloS one, 11(1), e0146443.
Molina, L. A., Skelin, I., & Gruber, A. J. (2014). Acute NMDA receptor antagonism disrupts synchronization of action potential firing in rat prefrontal cortex. PloS one, 9(1), e85842.
Placantonakis, D. G., & Welsh, J. P. (2001). Two distinct oscillatory states determined by the NMDA receptor in rat inferior olive. The Journal of physiology, 534(1), 123–140.
Rovó, Z., Mátyás, F., Barthó, P., Slézia, A., Lecci, S., Pellegrini, C. & Acsády, L. (2014). Phasic, nonsynaptic GABA-A receptor-mediated inhibition entrains thalamocortical oscillations. Journal of Neuroscience, 34(21), 7137–7147.
Vishnoi, S., Raisuddin, S., & Parvez, S. (2015). Modulatory effects of an NMDAR partial agonist in MK-801-induced memory impairment. Neuroscience, 311, 22–33.
