Script used to run neuronal network simulations in Burman et al. (2025) "Combined action of depolarising GABA and spike threshold adaptation underlies network dynamics in sleep-deprived cortex"

burman_et_al_2025_neuronal_network_simulation.py

#!/usr/bin/env python
# -*- coding: utf-8 -*

"""Simulate a LIF neuronal network under several conditions:

1) rested;
2) sleep deprived;
3) sleep deprived without spike threshold adaptation;
4) sleep deprived without E_GABA adaptation.

E_GABA, E_L, and V_th are varied between the conditions according to
experimentally observed variations in mouse pyramidal neurons in vivo:

|                      | Rested | Sleep deprived | SD without spike threshold adaptation | SD without E_GABA adaptation |
|----------------------+--------+----------------+---------------------------------------+------------------------------|
| Spike threshold (mV) |  -52.3 |          -45.5 |                                 -52.3 |                        -45.5 |
| E_GABA (mV)          |  -62.2 |          -52.1 |                                 -52.1 |                        -62.2 |
| E_leak (mV)          |  -63.2 |          -58.7 |                                 -58.7 |                        -58.7 |

Copyright (C) 2023 by Paul Brodersen.

This program is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or (at
your option) any later version. This program is distributed in the
hope that it will be useful, but WITHOUT ANY WARRANTY; without even
the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
PURPOSE.  See the GNU General Public License for more details. You
should have received a copy of the GNU General Public License along
with this program. If not, see <https://www.gnu.org/licenses/.

"""

import time
import pathlib
import numpy as np
import brian2 as b2
import pandas as pd

from argparse import ArgumentParser
from functools import partial
from uuid import uuid4
from scipy.signal import (
    welch,
    spectrogram as get_spectrogram,
    sosfilt,
    iirfilter,
)
from colorednoise import powerlaw_psd_gaussian # pip install colorednoise

parser = ArgumentParser()
parser.add_argument(
    '--balancing_time',
    help="Duration of balancing in seconds.",
    default=30.,
    type=float,
)
parser.add_argument(
    '--condition_time',
    help="Duration of each condition in seconds.",
    default=300.,
    type=float,
)
parser.add_argument(
    '--size',
    help="The total number of neurons.",
    default=1250,
    type=int,
)
parser.add_argument(
    '--stimulus',
    help="The mean magnitude of the stimulus current in pA.",
    default=50.,
    type=float,
)
parser.add_argument(
    '--noise',
    help="The mean magnitude of the noise current in pA.",
    default=50.,
    type=float,
)
parser.add_argument(
    '-o', '--output',
    help="/path/to/output.csv",
    default=None
)
parser.add_argument(
    '-d', '--data_directory',
    help="/path/to/data/directory",
    default=None
)

args = parser.parse_args()

################################################################################
# experimental conditions

balancing_time = int(args.balancing_time * 1000) # ms
condition_time = int(args.condition_time * 1000) # ms

# condition, spike threshold (mV), E_GABA (mV), E_L (mV), duration
parameters = [
    ("E/I balancing"                , -45.5,-52.1,-58.7, balancing_time, 0.1),
    ("Rested"                       , -52.3,-62.2,-63.2, condition_time, 0.),
    ("Sleep deprived"               , -45.5,-52.1,-58.7, condition_time, 0.),
    ("SD without spike adaptation"  , -52.3,-52.1,-58.7, condition_time, 0.),
    ("SD without E_GABA adaptation" , -45.5,-62.2,-58.7, condition_time, 0.),
]

conditions, vth, egaba , el, durations, plasticity = zip(*parameters)

################################################################################
# output summary data structure

uid = uuid4()

data = pd.DataFrame(
    {
        "uid" : uid,
        "Condition" : conditions,
        "Spike threshold [mV]" : vth,
        "E_GABA [mv]" : egaba,
        "Resting membrane potential [mV]" : el,
        "Duration [ms]" : durations,
        "Start [ms]" : np.r_[0, np.cumsum(durations)[:-1]],
        "Stop [ms]" : np.cumsum(durations),
        "Plasticity" : plasticity,
        "Network size" : args.size,
        "Stimulus [pA]" : args.stimulus,
        "Noise [pA]" : args.noise,
    }
)

################################################################################
# external input

b2.start_scope()

total_time = np.sum(durations)

# Construct a 1/f external input current.
# Use a strong stimulus during balancing to drive strong plasticity.
# Use a weaker stimulus during testing that does not cause population events.
stimulus_time_scale = 10 # ms
balancing_stimulus = 2 * args.stimulus * (1 + powerlaw_psd_gaussian(1, int(balancing_time                / stimulus_time_scale)))
condition_stimulus =     args.stimulus * (1 + powerlaw_psd_gaussian(1, int((total_time - balancing_time) / stimulus_time_scale)))

# Add random independent noise with the same mean magnitude as the stimulus.
noise_time_scale = 10 # ms
total_neurons = args.size
balancing_noise = 4 * args.stimulus * np.random.lognormal(size=(int(balancing_time                / noise_time_scale), total_neurons))
condition_noise = 2 * args.stimulus * np.random.lognormal(size=(int((total_time - balancing_time) / noise_time_scale), total_neurons))

I_stimulus = b2.TimedArray(-np.concatenate([balancing_stimulus, condition_stimulus], axis=0) * b2.pA, dt=stimulus_time_scale * b2.ms)
I_noise    = b2.TimedArray(-np.concatenate([balancing_noise,    condition_noise],    axis=0) * b2.pA, dt=noise_time_scale    * b2.ms)

################################################################################
# neuron model

params = dict()
params['Excitatory'] = {'name' :'Excitatory', 'count':int(0.8 * total_neurons), 'color':'crimson'}
params['Inhibitory'] = {'name' :'Inhibitory', 'count':int(0.2 * total_neurons), 'color':'cornflowerblue'}

dt = 1000 # ms
V_th_exc   = b2.TimedArray(np.repeat(vth,   durations) * b2.mV, dt=1 * b2.ms) # Spike threshold (excitatory)
E_GABA_exc = b2.TimedArray(np.repeat(egaba, durations) * b2.mV, dt=1 * b2.ms) # Inhibitory reversal potential (excitatory)
E_L_exc    = b2.TimedArray(np.repeat(el,    durations) * b2.mV, dt=1 * b2.ms) # Resting membrane potential (excitatory)

V_th_inh   = -50. * b2.mV        # Spiking threshold (inhibitory)
E_GABA_inh = -60. * b2.mV        # Inhibitory reversal potential (inhibitory)
E_L_inh    = -60. * b2.mV        # Resting membrane potential / reversal of the leak current (inhibitory)

E_GLUT     =   0. * b2.mV        # Excitatory reversal potential
tau_GLUT   =   5. * b2.ms        # Glutamatergic synaptic time constant
tau_GABA   =  10. * b2.ms        # GABAergic synaptic time constant
C_m        = 200. * b2.pfarad    # Membrane capacitance
g_L        =  10. * b2.nsiemens  # Leak conductance
t_ref_inh  =   1. * b2.ms        # Refractory period (inhibitory)
t_ref_exc  =   5. * b2.ms        # Refractory period (excitatory)

exc_neuron_eqs = '''
dv/dt      = -(I_L + I_GLUT + I_GABA + I_ext)/C_m : volt (unless refractory)
I_ext      =  I_stimulus(t) + I_noise(t, i)       : amp
I_L        =  g_L    * (v - E_L_exc(t))           : amp
I_GLUT     =  g_GLUT * (v - E_GLUT)               : amp
I_GABA     =  g_GABA * (v - E_GABA_exc(t))        : amp
dg_GLUT/dt = -g_GLUT / tau_GLUT                   : siemens
dg_GABA/dt = -g_GABA / tau_GABA                   : siemens
'''

pv_neuron_eqs = '''
dv/dt      = -(I_L + I_GLUT + I_GABA + I_ext)/C_m : volt (unless refractory)
I_ext      =  I_stimulus(t) + I_noise(t, i)       : amp
I_L        =  g_L    * (v - E_L_inh)              : amp
I_GLUT     =  g_GLUT * (v - E_GLUT)               : amp
I_GABA     =  g_GABA * (v - E_GABA_inh)           : amp
dg_GLUT/dt = -g_GLUT / tau_GLUT                   : siemens
dg_GABA/dt = -g_GABA / tau_GABA                   : siemens
'''

net = b2.Network(b2.collect())
population = dict()
population['Excitatory'] = b2.NeuronGroup(params['Excitatory']['count'],
                                          model      = exc_neuron_eqs,
                                          threshold  = 'v>V_th_exc(t)',
                                          reset      = 'v=E_L_exc(t)',
                                          refractory = t_ref_exc,
                                          method     = 'euler')
population['Inhibitory'] = b2.NeuronGroup(params['Inhibitory']['count'],
                                          model      = pv_neuron_eqs,
                                          threshold  = 'v>V_th_inh',
                                          reset      = 'v=E_L_inh',
                                          refractory = t_ref_inh,
                                          method     = 'euler')
population['Excitatory'].v = E_L_exc(0 * b2.ms)
population['Inhibitory'].v = E_L_inh
net.add(population)

################################################################################
# connectivity

static_exc_model = dict(model='w : siemens', on_pre='g_GLUT_post += w')
static_inh_model = dict(model='w : siemens', on_pre='g_GABA_post += w')

gmax      = 100.  * b2.nsiemens # Maximum inhibitory weight
tau_stdp  =  20.  * b2.ms       # STDP time constant
rho       =  20.  * b2.Hz       # Target excitatory population rate
beta      = rho * tau_stdp * 2  # Target rate parameter

vogels_model = dict(
    model='''
    w : siemens
    dApre/dt  = -Apre  / tau_stdp : siemens (event-driven)
    dApost/dt = -Apost / tau_stdp : siemens (event-driven)
    ''',
    on_pre='''
    Apre += 1. * nsiemens
    w = clip(w + (Apost - beta * nsiemens) * eta, 0 * nsiemens, gmax)
    g_GABA_post += w''',
    on_post='''
    Apost += 1. * nsiemens
    w = clip(w + Apre * eta, 0 * nsiemens, gmax)
    ''')

static_exc_synapse = partial(b2.Synapses, **static_exc_model)
static_inh_synapse = partial(b2.Synapses, **static_inh_model)
vogels_synapse     = partial(b2.Synapses, **vogels_model)

conn_params = dict()
conn_params[('Excitatory','Excitatory')] = dict(model=static_exc_synapse, p=0.1, w=1.*b2.nsiemens)
conn_params[('Excitatory','Inhibitory')] = dict(model=static_exc_synapse, p=0.6, w=1.*b2.nsiemens)
conn_params[('Inhibitory','Excitatory')] = dict(model=vogels_synapse,     p=0.6, w=1.*b2.nsiemens)
conn_params[('Inhibitory','Inhibitory')] = dict(model=static_inh_synapse, p=0.6, w=1.*b2.nsiemens)

connectivity = dict()
for connection in conn_params:
    connectivity[connection] = conn_params[connection]['model'](population[connection[0]], population[connection[1]])
    connectivity[connection].connect(p=conn_params[connection]['p'])
    connectivity[connection].w = conn_params[connection]['w']

net.add(connectivity)

################################################################################
# setup monitors

spike_monitor = dict()
for label in params.keys():
    spike_monitor[label] = b2.SpikeMonitor(population[label])
net.add(spike_monitor)

lfp_proxy = b2.NeuronGroup(1, 'synaptic_current : amp')
net.add(lfp_proxy)
lfp_aggregator = b2.Synapses(population["Excitatory"], lfp_proxy, 'synaptic_current_post = I_GLUT_pre / N_pre : amp (summed)')
lfp_aggregator.connect()
net.add(lfp_aggregator)
lfp_monitor = b2.StateMonitor(lfp_proxy, 'synaptic_current', record=True, dt=1*b2.ms)
net.add(lfp_monitor)

################################################################################
# run simulation

def get_human_readable_time(seconds):
    days    = seconds // 86400
    hours   = seconds // 3600 % 24
    minutes = seconds // 60 % 60
    seconds = seconds % 60

    msg = ""
    msg += f"{days} days, " if days else ""
    msg += f"{hours} hours, " if (days or hours) else ""
    msg += f"{minutes} minutes, " if (days or hours or minutes) else ""
    msg += f"{seconds:.1f} seconds"
    return msg

print("Simulation start:", time.strftime("%H:%M:%S", time.localtime()))
for condition, duration, eta in zip(conditions, durations, plasticity):
    print(f"{condition}...")
    tic = time.time()
    net.run(duration * b2.ms)
    toc = time.time()
    print(f"Time elapsed:", get_human_readable_time(toc-tic))

################################################################################
# save out results

if args.output: # save out summary data

    for population in params.keys():

        spike_times = spike_monitor[population].t / b2.ms
        spike_ids = spike_monitor[population].i
        total_neurons = params[population]["count"]

        def get_firing_rate(interval):
            start, stop = interval
            total_spikes = np.sum(np.logical_and(spike_times >= start, spike_times < stop))
            total_time_in_seconds = (stop - start) / 1000
            return total_spikes / total_time_in_seconds / total_neurons

        data[f"{population} neuron firing rate [Hz]"] = data[["Start [ms]", "Stop [ms]"]].apply(get_firing_rate, axis=1)

    print(data)

    output_path = pathlib.Path(args.output)
    if output_path.exists():
        old = pd.read_csv(output_path)
        df = pd.concat([old, data], axis=0)
    else:
        df = data
    df.to_csv(output_path, index=False)


if args.data_directory: # save out excitatory time series data

    lfp = np.c_[lfp_monitor.t / b2.ms, np.squeeze(lfp_monitor.synaptic_current / b2.pamp)]
    spikes = np.c_[spike_monitor["Excitatory"].t / b2.ms, spike_monitor["Excitatory"].i]
    np.savez(pathlib.Path(args.data_directory) / f"{uid}.npz", spikes=spikes, lfp=lfp)