This repository is associated with the following publication:
A unified model library maps how neuromodulation reshapes the excitability landscape of neurons across the brain
Domenico Guarino*, Ilaria Carannante*, and Alain Destexhe.
PLOS Computational Biology. 2025 Dec 1;21(12):e1013765.
https://doi.org/10.1371/journal.pcbi.1013765

The activity of neurons in the central nervous system changes drastically across different states like asleep or awake, attentive or drowsy, executing motor commands or resting. Neuromodulation supports these changes by affecting neuronal excitability.

Much work has been devoted to understanding the effects of neuromodulators on various neurons. However, we still lack a cohesive picture of how neuromodulation affects neuronal dynamics.
Here, we analyzed electrophysiological data from published papers and extracted features for seven types of neurons from five areas of the human and rodent central nervous system, under the effect of five neuromodulators.

We studied neuromodulation dynamics using the widespread, simple, yet biologically accurate, Adaptive Exponential Integrate-and-Fire (AdEx) model. We have found that different neuromodulators remap the feature space of neurons into non-overlapping clusters. This work organizes our knowledge about neuromodulation in a compact format, useful for modelers and theoreticians alike.

We systematically screened the available literaure using targeted PubMed searches. An overview of the recordings included in our dataset is here.

Key electrophysiological features are extracted from each trace.
For example by running the Jupyter Notebook Features_extraction.ipynb in the folder extracting_features (here), the user can extract key features from the voltage recording of a layer III cortical neuron.
The data was scraped using WebPlotDigitizer from Fig. 1B in McCormick and Williamson., 1989.
Running the notebook results in the generation of a .json file containing the extracted features.

A core mechanism of the algorithm is the simulation of the AdEx models, which consists in numerically solving the two coupled differential equations that describe the membrane potential $V(t)$ and adaptation variable $w(t)$ at each time step:

To do so, we implemented in Python the forward Euler numerical method Euler_current_protocol. To handle the spike discontinuity, if the membrane potential $V$ reaches or exceeds the peak voltage $V_\textit{peak}$, $V$ is reset to $V_\textit{reset}$ and $w$ is increased by $b$; then, for a refractory period (usually within $5$ ms), both quantities are kept constant.

The algorithm we developed takes as input the extracted features (reading them from the corresponding $\texttt{.json}$ files) and the user-defined ranges for the AdEx model parameters. Additionally, the user can specify the number of models to generate. We explore the parameter space by sampling each parameter within its specified bounds. To efficiently sample the parameter space we use Sobol's sequences via the python library $\texttt{scipy.stats.qmc.Sobol}$.
And example of such procedure can be fund here.
The user can run in parallel the AdEx_params_space.py file (the num_threads is set to 4 but can be easily changed), which creates as output a folder called test_cortex_control_v0 containing the best models and a pdf with a summary figure. To evaluate how well each AdEx model reproduces the extracted features and hence to rank the best models, an error function was implemented.

To visualize the output of the simulation, the user can run the notebook Cortical_human_pyramidal_neuron_control.ipynb.

We performed three main post-fitting analyses for each neuron type:

We compared the AdEx parameters distribution in control stage and under the influence of the different neuromodulators for each neuron type. Moreover, a jupyter notebook which collects all these plots can be found here. The output of this notebook constitutes Supplementary Material 2.

To visualize the structure of the parameter space and assess whether control and modulated models formed distinct groups, we performed a Principal Component Analysis.

For every parameter set, we extracted:

Following the analytical classification of Toubul and Brette (2008), the relative magnitudes of $R_a$ and $R_\tau$ uniquely determine the local bifurcation that terminates the resting state. In brief, a strong-and-fast adaptation brake ($R_a>R_\tau$) produces a sub-critical Hopf bifurcation (type-II onset), whereas a weak-or-slow brake ($R_a<R_\tau$) yields a saddle–node on invariant circle (type-I onset). Exact equality marks a codimension-2 Bogdanov–Takens (BT) point.

Since the AdEx model is described by two coupled differential equations, its dynamics can be easily visualized in a phase plane plot. In particular, we examine the nullclines:

namely the curves where the derivatives of $V$ and $w$ are zero.

Each "summary notebook" (for example this) contains all of these analyses.

The AdEx model allows us to predict the effect of neuromodulators across species. For example, to predict the effect of Dopamine on Human Cortincal pyramidal neurons, we used rodent data.
A notebook exemplifying this process is here.

To test the pipeline presented in this repository, the user can follow these steps:
0) Create a virtualenv (the files have been created using a virtualenv named 'nm') and use the requirements.txt file to install all you need to run the Jupyter notebooks;
1) Run this features extraction notebook;
2) Run this parameters exploration file. With the default setting the code runs in serial, but by setting parallel_running = True and defining the num_threads, the code runs in parallel;
3) Plot the results using this.