Neuromodcell is freely available to be downloaded from its Github site¹, written in Python with requirements specified in requirements.txt including NEURON² and Jupyter Notebook. To install Neuromodcell use pip install neuromodcell. The software is compatible with Linux operating system and super-computer clusters (Cray XC40 system). The simulations of large networks use Snudda (Hjorth et al., 2020, 2021), a Python package for creating data-driven networks of neurons available from its Github site³. To simulate neuromodulation, within the Snudda framework, additional simulation classes and associated neuromodulation subpackages were created, snudda.neuromodulation. To use Snudda, install via pip install snudda and follow the instructions on its wiki page (see text footnote 3).

Neuromodcell and Snudda require multicompartmental models specified for NEURON simulator described by three files: a morphology file (SWC), a parameter file (JSON) and a mechanism file (JSON) including a directory with the ion channel model files (.mod files, NEURON model description language). The models used in this publication were optimized using BluePyOpt (Van Geit et al., 2016) in previous publications (Hjorth et al., 2020; Lindroos and Hellgren Kotaleski, 2021) with modification for muscarinic modulation; and include direct and indirect striatal projection neuron (dSPN and iSPN) and fast-spiking interneuron (FS) models from the striatum.

To modulate the models during the simulation, we introduce additional parameters for the specific neuromodulator within the .mod files. The modulation is implemented in a phenomenological manner, where a modulation parameter is multiplied with the target (for example the conductance of an ion channel). Hence, the conductance (in this case) can be regulated through the simulation. For example, for dopamine, the three parameters are maxModDA, modDA and levelDA (using the key, DA to indicate dopamine), based on previous formalism developed in Hjorth et al. (2020). In general, the convention is ‘maxMod^∗’, ‘level^∗’ and ‘mod^∗’ and the ^∗ should be replaced with a specific name for the neuromodulator. The maxMod^∗ defines the modulation degree (which can vary, with 0.6 being 40% reduction). The level^∗ defines the transient level of modulation throughout the simulation. The mod^∗ defines the activation of modulation (0 inactive, 1 fully active). Within the examples, dopaminergic and muscarinic modulation was included (with muscarinic modulation using key ACh). The list of modulated ion channels for each cell type is presented in Supplementary Table 1. For examples of .mod-files with these modulation parameters (see text footnote 1).

Neuromodcell is a software for optimizing and simulating neuromodulation in detailed multicompartmental neuron models. The software loads and simulates a population of models with user-defined parameter variations of maxMod^∗, which is how the program induces the changes associated with neuromodulation. The software also incorporates current clamp and voltage clamp parameters to simulate specific biological experiments. Hence, the user can define the simulated experiment in terms of both neuromodulator and the specific protocol to be used. Following the simulation, the parameter sets that reproduce the experimental data are selected based on a user-specified selection criteria. We provide several electrophysiological features by default, but custom features can be included.

Here, Neuromodcell was used to optimize dopaminergic modulation for dSPN, iSPN and FS, as well as muscarinic modulation for dSPN and iSPN. The optimizations for dSPN and iSPN were validated against data from Planert et al. (2013) and Lahiri and Bevan (2020) (see Supplementary Figure 1) for dopaminergic modulation and Lv et al. (2017) for muscarinic modulation. The dopaminergic modulation for FS was validated against data from Bracci et al. (2002).

Neuromodulation in Snudda is implemented as a separate simulation module called snudda.neuromodulation. Within the module, there are two simulation classes for neuromodulation, SnuddaSimulateNeuromodulation and SnuddaSimulateNeuromodulationSynapse. These classes represent the two strategies used to simulate neuromodulation, replay and adaptive. The replay mode plays a predetermined transient through the level^∗-parameters and hence modulates each component of the circuit in a predefined time dependent manner. On the other hand, the adaptive mode originates from the fact that a neuromodulatory neuron integrated within the microcircuit could not modulate via a predetermined transient. Instead, the spiking activity of such a neuron would have to be continuously translated into an instantaneous level of neuromodulation. Therefore, an alternative approach was developed where an intermediate mechanism would integrate the spiking activity. Adaptive simulations couple one or several presynaptic neurons to postsynaptic neurons via an integrating mechanism called conc^∗, which modifies the neuromodulatory parameters (level^∗) in the circuit (for example concDA to levelDA for dopamine). Hence, neuromodulation within the simulation can either follow a fixed recipe (‘replay’ mode), or dynamically change based on the activity in the network (‘adaptive’ mode). The mode used within a particular simulation will therefore depend on the network structure. Using either approach, the level of neuromodulation would be updated every time step via the level^∗-parameter. The replay mode loads an array which has set values for the level^∗ for each time step of the simulation. This would suffice if the neuromodulation is synchronous and the neuromodulatory neurons are not affected by activity of the circuit. On the other hand, if the activity of the microcircuit could for example inhibit the neuromodulatory neurons, like striosomal dSPNs in the striatum, replay mode does not suffice and adaptive mode would be used instead, as the presynaptic neuromodulatory neurons can interact with the microcircuit throughout the simulation.

For replay, the maxMod^∗, mod^∗ and level^∗ are introduced as RANGE variables. In the adaptive mode, which is based on pointers, the parameters are introduced with the keyword POINTERS⁴ and ‘ptr’ is added to the filename (i.e., na.mod to na_ptr.mod).

Lastly, both types of simulations require a configuration file (JSON) which is created using either Neuromodulation class (snudda.neuromodulation.modulation_network, for replay simulations) or NeuromodulationSynapse (snudda.neuromodulation.modulation_synapse, for adaptive simulations). For examples of how to use these classes (see text footnote 4).

To demonstrate the features of Snudda.neuromodulation, we simulated several networks using both replay and adaptive mode. Firstly, using the Snudda framework, a network of 10 000 neurons was created and simulated with dopamine and acetylcholine transients; and with a cortical activation at 1 s, with cortical and thalamic background activity. The network is based on previous publications (Hjorth et al., 2020) and the connection probabilities follows the diagram in Supplementary Figure 2. The transients consisted of a tonic level of ACh accompanied with a burst or pause of ACh and a DA burst. The simulations were performed on a super-computer cluster (Cray XC40 system at PDC Center for High Performance Computing, KTH). A smaller network with dSPNs was simulated to show the effect of DA and ACh without network interactions from iSPNs and FS. Furthermore, several dopamine transients with different start times were created using the Neuromodulation class and applied to a network of 20 dSPNs with a cortical activation lasting for 500 ms. Lastly, a network of dSPNs was simulated with presynaptic spiking neurons using the adaptive mode. The presynaptic neurons represented dopaminergic neurons and were simulated with and without activation (i.e., bursting). The activity of the neurons within the simulations were measured by the percentage of spiking neurons and the mean firing frequency during the cortical activations using custom Python code and Electrophysiology Analysis Toolkit⁵.

A tutorial on Neuromodcell and an example of dopaminergic modulation of dSPN are available (see text footnote 1).

Examples of simulations using Snudda.neuromodulation module are available (see text footnote 4).

Neuromodcell performs parameter variation of ion channels and receptor models within multicompartmental models to simulate neuromodulation. The components of the software are presented in Figure 1 with the main class, DefineModulation (Figure 2) and its associated class methods, which are used to define the parameters, protocols, selection criteria and modulation for the specific cell and neuromodulator. The user specifies the range of parameter variation for each ion channel and/or receptor as well as the experimental protocol to be replicated in the simulation (Figures 3, 4). The models with parameter variations are simulated using the optimise module, which uses mpi4py for parallelization.

Following the simulation, a selection process is performed, where the models are compared to the experimental data provided in the setup (using DefineModulation). The parameter variations of the models which pass validation are saved in a separate file. The code for defining the selection criteria is freely available (module selection_criteria and selection_function) and can be modified to measure a specific electrophysiological feature of the simulation, for example firing frequency or the number of action potentials.

Below, we describe the process of setting up, optimizing and selecting the models that pass validation. Following the completion of the optimization, the results are transferred to Snudda and simulated using snudda.neuromodulation module. A step-by-step tutorial is available (see text footnote 4).

We performed optimizations of dopaminergic and muscarinic modulation for dSPN, iSPN, and FS. In the following examples, the dopaminergic optimization of dSPNs is used to demonstrate the features of Neuromodcell. We instantiated the DefineModulation class as shown in Figure 2. The class requires several arguments to define the optimization, including the model directory and the number of model variations to simulate (population). The next step is parameter variation. In Figure 3, the ion channels which are modulated by DA within dSPN (based on previous publications Hjorth et al., 2020 and Lindroos et al., 2018) were set using the define_modulation_parameter method within the class. The modulation can also be defined to be somatic, dendritic (basal or apical, due to the SWC naming convention) or axonal.

The experimental data used in the example optimization was taken from Figure 6C in Planert et al. (2013) for dSPNs. The experiment used patch clamp recordings and bath application of dopamine to measure the effect of DA on dSPNs. Following the application of DA, there was an increase in the intrinsic excitability of dSPNs measured by the number of action potentials. Hence, an experimental setup was introduced into modset, with a current clamp protocol and bath application of dopamine (Figure 4). Furthermore, the method define_selection_criteria was used to introduce the measurement used in the validation step (i.e., the change in the number of action potentials). Both the define_modulation_function and define_selection_criteria can select several pre-defined functions but they can also be customized to accommodate specific transients or criteria, respectively.

The optimization uses the OptimiseModulation, from optimise module in Neuromodcell. The optimization for dSPN used a custom Python script called optimise_dspn.py, which loaded the path_to_model and the specific seed for the optimization (Figure 5). The seed is used to randomize the modulation parameters defined using define_modulation_parameter. The modulation parameters sets (i.e., population) are simulated in parallel with the specific model (defined in DefineModulation), as illustrated in Figure 1.

The result of the optimization was visualized using dSPNanalysis (inheritance from OptimiseAnalysis class in Neuromodcell) (Figure 6). The class loads the output files following the optimization and enables the user to plot and/or analyze the results. In Figure 7, the result of dopaminergic modulation of dSPN is shown with the model variations (modulated) and control simulation (in black). The dopaminergic modulation increases the intrinsic excitability of dSPNs via D1 receptors. The result of the Neuromodcell optimization produced model variations which reproduce the experimental data from Planert et al. (2013). The control simulation (in black, Figure 7A) has four spikes following current injection. The simulated dopamine modulation increases the number of spikes and the change in the number of action potentials are within the range extracted from the experimental data (Figure 7B).

Following the optimization, the parameter sets which reproduce the neuromodulation are saved as modulation.json and moved to the model directory (i.e., the model is now defined by parameters.json, morphology.swc, mechanisms.json and modulation.json). The modified .mod files should also be transferred to a common Snudda mechanism directory.

Using the replay mode in Snudda.neuromodulation, we simulated a network of 10 000 (Figure 8A) multicompartmental neuron models of dSPNs, iSPNs, and FS to exemplify the effect of the dopamine modulation. The simulation included a cortical stimulation of dSPNs, iSPNs and FS at 1 s, lasting for 0.5 s. In Figure 8, the effect of dopamine on network activity is shown. Following cortical activation, 6% of dSPNs within the network were spiking. By modulating the ion channels defined in the optimization result of Neuromodcell, dopamine modulation increases the percentage of spiking dSPNs to the double (Figures 8B,C). iSPNs are modulated by D2 receptors, which reduces the intrinsic excitability. Hence, in Figure 8B, the percentage of spiking iSPNs within the network decreases by approximately 3%. The middle panel in Figure 8C shows an example trace of an iSPN, where the dopamine modulation reduces the number of spikes. FS are modulated by D1-like receptors, which depolarizes the membrane potential, as seen in Figure 8C. Following, an optimization using Neuromodcell, the parameter sets for FS dopamine modulation increase the number of spiking neurons within the network (Figure 8B). The cortical activation occurs via activating glutamate ionotropic receptor models. These receptors can also be modulated by the same formalism as previously described which leads to changes in the amplitude and release probability (Supplementary Figure 3).

The neuromodulation extension to Snudda can also simulate multiple neuromodulators. To demonstrate this, we simulated five scenarios of dopamine and acetylcholine modulation with the same network presented in Figure 8A. As previously described, dSPNs and FS are modulated by D1 receptors while iSPNs are modulated by D2 receptors. In addition, we used Neuromodcell to optimize for muscarinic modulation for dSPNs and iSPNs.

In Figure 9A, we show the transients used in the simulations, which consisted of an acetylcholine burst or pause and a dopamine burst. In the striatum, acetylcholine is released continuously by cholinergic interneurons due to their spontaneous activity. Hence, the acetylcholine transients had a tonic modulation level of 0.5 with burst and pause response causing an increase or decrease (by 0.5), respectively. These transients were simulated individually but also simultaneously (DA and ACh burst; DA and ACh pause). As shown in Figure 8, dopamine increased the percentage of spiking dSPNs within the network. The acetylcholine burst also slightly increases the percentage of spiking dSPNs, although in the simultaneous simulation (Figure 9B, light purple), the effect is not additive. Instead, the dopamine burst and acetylcholine pause produces a larger response. This is contrary to previous results in Lindroos and Hellgren Kotaleski (2021), where DA and ACh had additive effects when single neurons were simulated. Hence, we performed a simulation with only dSPNs to investigate the effect of DA and ACh. Without network interactions, the effect of DA and ACh was indeed additive (Supplementary Figures 4A,C). Hence, the decrease in the percentage of spiking dSPNs could be attributed to the network interactions and the changes seen in iSPNs and FS. As described in Figure 8, the intrinsic excitability of iSPNs is decreased following dopamine modulation, while ACh modulates iSPNs via M1 receptor, which lead to the opposite effect. Compared to dSPNs, iSPNs have a larger increase in the number of spiking neurons following an ACh burst (Figure 9C). During the ACh pause, the percentage of spiking iSPNs is expected to decreased, due to the reduction in ACh levels. In contrast, the percentage of spiking iSPNs was larger than in the control. In Zucca et al. (2018), SPNs were recorded in vivo and cholinergic interneurons were inhibited optogenetically. They observed a reduction in the spontaneous activity of SPNs. The control simulation in Figure 9 was, however, run without any neuromodulation, which could explain the inconsistency. Hence, we simulated a small network of SPNs with a cholinergic pause, but we modified the control simulation to contain a continuous tonic level of ACh. By comparing the mean firing frequency of SPNs with and without the pause, we observed the expected reduction in SPN activity (Supplementary Figure 4). Therefore, the tonic level of ACh has an influence on the activity of SPNs. For FS, the main effect is seen within the dopamine burst modulation (Figure 9D), while for ACh pause and burst there is no visible effect (due to the lack of direct muscarinic modulation of FS).

The user can construct different transients to investigate the effect of timing on network activity. We created a smaller network of 30 dSPNs to illustrate the effect such an investigation could have on network activity. We defined six transients which were identical, in terms of time constants and modulation level (Figure 10A). We then adjusted the start time in relation to cortical stimulation at 1s (of a 3 s simulation). The different delays caused the peak modulation to occur at different times during the cortical stimulation. Figure 10B shows that the transient −100 ms results in the largest increase in mean firing frequency. While −/+500 ms transients are not different from control. This is due to the dopamine level being too low and the peak dopamine response occurring much earlier and later compared to the cortical stimulation. The +100 ms delay produced a lower increase in the mean firing frequency but the percentage of spiking neurons were similar to the −100 ms case (Figure 10C). The user could conduct similar investigations into other features of the transients such as the time constants, tonic levels or more complex transients (Lahiri and Bevan, 2020).

The previous simulations have used a transient(s) to modulate the network. In Figure 11, we simulate neuromodulation using the adaptive mode in Snudda. In adaptive mode, the spiking activity of the presynaptic neurons are translated to modulation levels as described in Figure 11A. The top panel in Figure 11A shows the spike times of two presynaptic neurons (green and red). These neurons are connected to an integrating mechanisms called concDA (for dopamine, in orange in Figure 11A). The concDA integrates the spiking activity of the two cells into modulation levels (as seen in Figure 11A in left most box). The concDA sends the modulation level from both neurons into the ion channel model (right most box in Figure 11A). The level of current follows the modulation as the conductance of the ion channel is modulated and hence affecting the neuronal activity.

We created a small network of dSPNs to illustrate the effect of changing the firing frequency of the presynaptic neurons. We placed spontaneously active neurons to modulate the dSPNs via D1 receptors. In control, 50% of the dSPNs within the network were spiking and the spontaneous activity of the presynaptic neurons did not change the percentage (Figure 11D). Although in Figure 11B, the spontaneous activity has an effect on the membrane potential of the dSPNs. By stimulation the presynaptic spontaneous neurons, we increased the firing frequency (Figure 11C). As expected the level of neuromodulation increased as the concDA mechanism integrates the number of spikes into neuromodulation. This leads to a depolarization of the dSPN membrane potential. On a network scale, the increased firing frequency leads to simulated D1 receptor modulation and an increase in the percentage of spiking dSPNs.