Code to run neuronal network simulation in Alfonsa et al (2020)

alfonsa_2020_neuronal_network_simulation.py

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

"""
Determine the effect of changes in the GABA reversal potential in
pyramidal neurons on their ability to synchronize.

Copyright (C) 2020 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 pathlib
import numpy as np
import matplotlib.pyplot as plt
import brian2 as b2
import pandas as pd

from functools import partial
from analytics import ProgressBar

import matplotlib as mpl
mpl.rcParams['axes.spines.top'] = False
mpl.rcParams['axes.spines.right'] = False


###########################################################
# PARSE COMMAND LINE ARGUMENTS

# We would like to run this script multiple times with different parameterization.
# However, the brian magic makes functionalisation and encpasulation of the simulation near impossible.
# Normally, we could still run multiple simulations within the same script using the 'store' and 'restore' brian functions.
# However, those functions do not reset the simulation clock.
# Hence the TimedArray class used below for input does not properly support running a simulation multiple times in a row.
# Therefor we use an external script to call this script multiple times, and pass in cmd arguments.

from argparse import ArgumentParser

parser = ArgumentParser()
parser.add_argument(
    '-o', '--output',
    help="/path/to/output.csv",
    default='output.csv')
parser.add_argument(
    '--hz',
    help="Dominant frequency of rhythmic input",
    default=1,
    type=int,
)
parser.add_argument(
    '-c', '--color',
    help="Plot color",
    default='#1f77b4',
    type=str,
)
args = parser.parse_args()

output_filepath = args.output
hz = args.hz
color = args.color

###########################################################
# DEFINE NETWORK PARAMETERS

b2.start_scope()

total_neurons = 500

# Define network architecture
params = dict()
params['exc'] = {'name' :'excitatory', 'count':int(4* total_neurons / 5), 'color':'crimson'}
params['inh'] = {'name' :'inhibitory', 'count':int(   total_neurons / 5), 'color':'cornflowerblue'}

# Define neuron's characteristics
E_L       = -60.   *b2.mV       # Resting potential
# er        = -80.   *b2.mV       # Inhibitory reversal potential
er_inh    = -80.   *b2.mV       # Inhibitory reversal potential for inhibitory neurons
V_th      = -50.   *b2.mV       # Spiking threshold
t_ref     =   5.   *b2.ms       # Refractory period
gl        =  10.   *b2.nsiemens # Leak conductance
C_m       = 200.   *b2.pfarad   # Membrane capacitance
tau_m     = C_m/gl              # Membrane time constant
gmax      = 100.   *b2.nsiemens # Maximum inhibitory weight
weight_ex =   0.3  *b2.nsiemens # Excitatory weight
weight_in =   3.   *b2.nsiemens # Inhibitory weight

# --------------------------------------------------------------------------------
# Main experimentally varied parameter:
# Inhibitory reversal potential for excitatory neurons

epoch_duration = 1000
er_exc_train = np.array([-80])
er_exc_test = np.array([-80, -77.5, -75, -72.5, -70, -67.5, -65, -62.5, -60])
er_exc    = b2.TimedArray(np.r_[er_exc_train, er_exc_test]*b2.mV, dt=epoch_duration*b2.ms)

# --------------------------------------------------------------------------------

eqs_excitatory_neurons='''
dv/dt    = -(v - E_L)/tau_m - I/C_m + I_ext/C_m : volt (unless refractory)
dg_ex/dt = -g_ex/tau_ex_decay                   : siemens
dg_in/dt = -g_in/tau_in_decay                   : siemens
I_ext    =  I_rythmic(t) + I_random(t, i)       : amp
I_ex     =  g_ex*v                              : amp
I_in     =  g_in*(v - er_exc(t))                : amp
I        =  I_ex + I_in                         : amp
'''

eqs_inhibitory_neurons='''
dv/dt    = -(v - E_L)/tau_m - I/C_m + I_ext/C_m : volt (unless refractory)
dg_ex/dt = -g_ex/tau_ex_decay                   : siemens
dg_in/dt = -g_in/tau_in_decay                   : siemens
I_ext    =  I_rythmic(t) + I_random(t, i)       : amp
I_ex     =  g_ex*v                              : amp
I_in     =  g_in*(v - er_inh)                   : amp
I        =  I_ex + I_in                         : amp
'''

net = b2.Network(b2.collect())

population = dict()
population['exc'] = b2.NeuronGroup(params['exc']['count'],
                                   model      = eqs_excitatory_neurons,
                                   threshold  = 'v > V_th',
                                   reset      = 'v=E_L',
                                   refractory = t_ref,
                                   method     = 'euler'
)
population['inh'] = b2.NeuronGroup(params['inh']['count'],
                                   model      = eqs_inhibitory_neurons,
                                   threshold  = 'v > V_th',
                                   reset      = 'v=E_L',
                                   refractory = t_ref,
                                   method     = 'euler'
)
net.add(population)

###########################################################
# DEFINE CONNECTIVITY

# Synapse parameters
tau_ex_decay =  5.*b2.ms      # Glutamatergic synaptic time constant
tau_in_decay = 10.*b2.ms      # GABAergic synaptic time constant
tau_stdp     = 20.*b2.ms      # STDP time constant
rho          =  3.*b2.Hz      # Target excitatory population rate
beta         = rho*tau_stdp*2 # Target rate parameter

# Connection probabilities
p_e = 0.05 # excitatory connections
p_i = 0.05 # inhibitory connections

# Define synapses
static_ex_model = dict(model='w : siemens', on_pre='g_ex_post += w')
static_in_model = dict(model='w : siemens', on_pre='g_in_post += w')
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*nS)*eta, 0, gmax)
                g_in_post += w''',
           on_post='''
                Apost += 1.*nsiemens
                w = clip(w+Apre*eta, 0, gmax)
                ''')

static_ex_synapse = partial(b2.Synapses, **static_ex_model)
static_in_synapse = partial(b2.Synapses, **static_in_model)
vogels_synapse =    partial(b2.Synapses, **vogels_model)

# Create connections
conn_params = dict()
conn_params[('exc','exc')] = dict(model=static_ex_synapse, p=p_e, w=weight_ex)
conn_params[('exc','inh')] = dict(model=static_ex_synapse, p=p_e, w=weight_ex)
conn_params[('inh','exc')] = dict(model=vogels_synapse,    p=p_e, w=1e-10*b2.nsiemens)
# conn_params[('inh','inh')] = dict(model=static_in_synapse, p=p_e, w=weight_in)

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

# Create spike, state and rate monitors
spike_monitor = dict()
state_monitor = dict()
rate_monitor  = dict()
for label in params.keys():
    spike_monitor[label] = b2.SpikeMonitor(population[label])
    state_monitor[label] = b2.StateMonitor(population[label], variables=True, record=True)
    rate_monitor[label]  = b2.PopulationRateMonitor(population[label])

net.add(spike_monitor)
net.add(state_monitor)
net.add(rate_monitor)

# ###########################################
# DEFINE EXTERNAL INPUT

total_training_epochs = len(er_exc_train)
total_testing_epochs = len(er_exc_test)
training_time = total_training_epochs * epoch_duration
testing_time  = total_testing_epochs * epoch_duration
total_time    = training_time + testing_time

# rythmic external input
up_down_cycle_duration = 1000. / hz
up_down_cycle_amplitudes = [50., 0.]
total_up_down_cycles  = int((total_time + up_down_cycle_duration)/up_down_cycle_duration / len(up_down_cycle_amplitudes))
arr_rythmic = np.tile(up_down_cycle_amplitudes, 2*total_up_down_cycles)
I_rythmic = b2.TimedArray(arr_rythmic*b2.pA, dt=up_down_cycle_duration/2*b2.ms)

# random background noise (time-varying)
noise_time_scale = 10.
arr_noise = 50 + 50 * np.random.randn(int(total_time / noise_time_scale), total_neurons)
I_random  = b2.TimedArray(arr_noise*b2.pA, dt=noise_time_scale*b2.ms)

# ###########################################
# RUN SIMULATION

for simulation_time, eta in zip((training_time, testing_time), [0.1, 0.]):
    net.run(simulation_time * b2.ms, report=ProgressBar(), report_period=10*b2.ms)

# ###########################################
# PLOT

fig, (ax1, ax2, ax3, ax4) = plt.subplots(4, 1, sharex=True, figsize=(6.85, 10))

# plot inhibitory reversal potential
x = np.arange(training_time, total_time)
y = er_exc(x * b2.ms)
ax1.plot(x, y / b2.mV, color='black', alpha=0.9)
ax1.set_ylabel('Inhibitory reversal potential [mV]')
ax1.set_xlim(training_time, total_time)

# plot neuronal spike times
ax2.plot(spike_monitor['exc'].t/b2.ms, spike_monitor['exc'].i, '.',
         markersize = 0.5,
         # color      = params['exc']['color'],
         color      = color,
         alpha      = 0.5,
         rasterized = True)
# ax2.set_ylabel(params['exc']['name'])
ax2.set_ylabel('Excitatory neurons')

# plot population firing rates
rate_time_points = rate_monitor['exc'].t/b2.ms
rate_values = rate_monitor['exc'].smooth_rate(window = 'gaussian', width = 10*b2.ms)/b2.Hz
ax3.plot(rate_time_points,
         rate_values,
         # label = params['exc']['name'],
         # color = params['exc']['color'],
         color = color,
         alpha = 0.9,
)
ax3.set_ylabel('Firing rate [Hz]')

# determine ON/DOWN ratios
up_down_cycle_boundaries = np.arange(training_time, total_time+1, up_down_cycle_duration)
up_down_cycle_midpoints = up_down_cycle_boundaries[:-1] + up_down_cycle_duration / 2
up_intervals = np.c_[up_down_cycle_boundaries[:-1], up_down_cycle_midpoints]
down_intervals = np.c_[up_down_cycle_midpoints, up_down_cycle_boundaries[1:]]

is_up = np.zeros_like(rate_time_points, dtype=np.bool)
for ii, (start, stop) in enumerate(up_intervals):
    mask = np.logical_and(rate_time_points >= start, rate_time_points < stop)
    is_up += mask

is_down = np.zeros_like(rate_time_points, dtype=np.bool)
for ii, (start, stop) in enumerate(down_intervals):
    mask = np.logical_and(rate_time_points >= start, rate_time_points < stop)
    is_down += mask

epoch_boundaries = np.arange(training_time, total_time+1, epoch_duration)
epoch_intervals = np.c_[epoch_boundaries[:-1], epoch_boundaries[1:]]

up_rates = np.zeros((total_testing_epochs))
down_rates = np.zeros((total_testing_epochs))
for ii, (start, stop) in enumerate(epoch_intervals):
    mask = np.logical_and(rate_time_points >= start, rate_time_points < stop)
    up_rates[ii]   = np.mean(rate_values[np.logical_and(is_up,   mask)])
    down_rates[ii] = np.mean(rate_values[np.logical_and(is_down, mask)])

up_down_ratio = up_rates / down_rates

ax4.bar(np.mean(epoch_intervals, axis=1), up_down_ratio, width=epoch_duration, facecolor='None', edgecolor=color)
ax4.set_ylabel(r'Firing rate $\frac{UP}{DOWN}$')
ax4.set_xlabel('Time [ms]')

fig.tight_layout()
fig.savefig(f'../figures/proof_of_principle--{hz}_hz.pdf')
fig.savefig(f'../figures/proof_of_principle--{hz}_hz.png')
fig.savefig(f'../figures/proof_of_principle--{hz}_hz.svg')

data = dict(
    up_rates = up_rates,
    down_rates = down_rates,
    up_down_ratio = up_down_ratio,
    gaba_reversal = er_exc_test,
    hz = hz * np.ones_like(er_exc_test),
)

output_path = pathlib.Path(output_filepath)
if output_path.exists():
    df1 = pd.read_csv(output_path)
    df2 = pd.DataFrame(data)
    df = pd.concat([df1, df2], sort=True)
else:
    df = pd.DataFrame(data)

df.to_csv(output_path, index=False)

plt.show()

Hi, I am trying to reproduce the code, but I get an error: "brian2.units.fundamentalunits.DimensionMismatchError: Argument number 2 for function clip was supposed to have the same units as argument number 1, but '0.0' has unit 1, while 'w + ((Apost - (beta * nS)) * eta)' has unit S
". Might you help me?

Hi, thanks for raising the issue. Apparently, some versions of brian complain about dimensionless zeros, others don't. Try substituting l. 173 with

w = clip(w+(Apost-beta*nS)*eta, 0*nS, gmax)

and l. 177 with

w = clip(w+Apre*eta, 0*nS, gmax)