Digital Drosophila Archive
374 lines
12 KiB
Python
374 lines
12 KiB
Python
import pandas as pd
|
|
from textwrap import dedent
|
|
|
|
# brian 2
|
|
from brian2 import NeuronGroup, Synapses, PoissonInput, SpikeMonitor, Network
|
|
from brian2 import mV, ms, Hz
|
|
|
|
# file handling
|
|
from pathlib import Path
|
|
|
|
# parallelization
|
|
from joblib import Parallel, delayed, parallel_backend
|
|
from time import time
|
|
|
|
default_params = {
|
|
# trials
|
|
't_run' : 1000 * ms, # duration of trial
|
|
'n_run' : 30, # number of runs
|
|
|
|
# network constants
|
|
# Kakaria and de Bivort 2017 https://doi.org/10.3389/fnbeh.2017.00008
|
|
'v_0' : -52 * mV, # resting potential
|
|
'v_rst' : -52 * mV, # reset potential after spike
|
|
'v_th' : -45 * mV, # threshold for spiking
|
|
't_mbr' : 20 * ms, # membrane time scale (capacitance * resistance = .002 * uF * 10. * Mohm)
|
|
|
|
# Jürgensen et al https://doi.org/10.1088/2634-4386/ac3ba6
|
|
'tau' : 5 * ms, # time constant
|
|
|
|
# Lazar et al https://doi.org/10.7554/eLife.62362
|
|
't_rfc' : 2.2 * ms, # refractory period
|
|
|
|
# Paul et al 2015 doi: 10.3389/fncel.2015.00029
|
|
't_dly' : 1.8*ms, # delay for changes in post-synaptic neuron
|
|
|
|
# Free parameter
|
|
'w_syn' : .275 * mV, # weight per synapse (note: modulated by exponential decay)
|
|
# Default activation rates
|
|
'r_poi' : 150*Hz, # default rate of the Poisson inputs
|
|
'r_poi2' : 0*Hz, # default rate of a 2nd class of Poisson inputs
|
|
'f_poi' : 250, # scaling factor for Poisson synapse; 250 is sufficient to cause spiking
|
|
|
|
# equations for neurons # alpha synapse https://doi.org/10.1017/CBO9780511815706; See https://brian2.readthedocs.io/en/stable/user/converting_from_integrated_form.html
|
|
'eqs' : dedent('''
|
|
dv/dt = (v_0 - v + g) / t_mbr : volt (unless refractory)
|
|
dg/dt = -g / tau : volt (unless refractory)
|
|
rfc : second
|
|
'''),
|
|
# condition for spike
|
|
'eq_th' : 'v > v_th',
|
|
# rules for spike
|
|
'eq_rst' : 'v = v_rst; w = 0; g = 0 * mV',
|
|
}
|
|
|
|
|
|
#######################
|
|
# brian2 model setup
|
|
def poi(neu, exc, exc2, params):
|
|
'''Create PoissonInput for neurons.
|
|
|
|
For each neuron in 'names' a PoissonInput is generated and
|
|
the refractory period of that neuron is set to 0 in NeuronGroup.
|
|
|
|
Parameters
|
|
----------
|
|
neu : NeuronGroup
|
|
Defined brian2.NeuronGroup object
|
|
exc : list
|
|
Indices of neurons for which to create Poisson input
|
|
exc2 : list
|
|
Indices of neurons for which to create Poisson input with `r_poi`
|
|
params : dict
|
|
Constants and equations that are used to construct the brian2 network model
|
|
|
|
Returns
|
|
-------
|
|
pois : list
|
|
PoissonInput objects for each neuron in 'exc'
|
|
neu : NeuronGroup
|
|
NeuronGroup with adjusted refractory periods
|
|
'''
|
|
|
|
pois = []
|
|
for i in exc:
|
|
p = PoissonInput(
|
|
target=neu[i],
|
|
target_var='v',
|
|
N=1,
|
|
rate=params['r_poi'],
|
|
weight=params['w_syn']*params['f_poi']
|
|
)
|
|
neu[i].rfc = 0 * ms # no refractory period for Poisson targets
|
|
pois.append(p)
|
|
|
|
for i in exc2:
|
|
p = PoissonInput(
|
|
target=neu[i],
|
|
target_var='v',
|
|
N=1,
|
|
rate=params['r_poi2'],
|
|
weight=params['w_syn']*params['f_poi']
|
|
)
|
|
neu[i].rfc = 0 * ms # no refractory period for Poisson targets
|
|
pois.append(p)
|
|
|
|
return pois, neu
|
|
|
|
def silence(slnc, syn):
|
|
'''Silence neuron by setting weights of all synapses from it to 0
|
|
|
|
Parameters
|
|
----------
|
|
slnc : list
|
|
List of neuron indices to silence
|
|
syn : brian2.Synapses
|
|
Defined synapses object
|
|
|
|
Returns
|
|
-------
|
|
syn : brian2.Synapses
|
|
Synapses with modified weights
|
|
'''
|
|
|
|
for i in slnc:
|
|
syn.w[' {} == i'.format(i)] = 0*mV
|
|
|
|
return syn
|
|
|
|
def create_model(path_comp, path_con, params):
|
|
'''Create default network model.
|
|
|
|
Convert the "completeness materialization" and "connectivity" dataframes
|
|
into a brian2 neural network model. Network constants and equations
|
|
are defined at the beginning of this file.
|
|
|
|
Parameters
|
|
----------
|
|
path_comp : str
|
|
path to "completeness materialization" dataframe
|
|
path_con : str
|
|
path to "connectivity" dataframe
|
|
params : dict
|
|
Constants and equations that are used to construct the brian2 network model
|
|
|
|
|
|
Returns
|
|
-------
|
|
neu : NeuronGroup
|
|
brian2.NeuronGroup object with neurons as in 'path_comp'
|
|
syn : Synapses
|
|
brian2.Synapses object with connections as in 'path_con'
|
|
spk_mon : SpikeMonitor
|
|
brian2.SpikeMonitor object, which records time of spike events
|
|
'''
|
|
|
|
# load neuron connectivity dataframes
|
|
df_comp = pd.read_csv(path_comp, index_col=0)
|
|
df_con = pd.read_parquet(path_con)
|
|
|
|
neu = NeuronGroup( # create neurons
|
|
N=len(df_comp),
|
|
model=params['eqs'],
|
|
method='linear',
|
|
threshold=params['eq_th'],
|
|
reset=params['eq_rst'],
|
|
refractory='rfc',
|
|
name='default_neurons',
|
|
namespace=params,
|
|
)
|
|
neu.v = params['v_0'] # initialize values
|
|
neu.g = 0
|
|
neu.rfc = params['t_rfc']
|
|
|
|
# create synapses
|
|
syn = Synapses(neu, neu, 'w : volt', on_pre='g += w', delay=params['t_dly'], name='default_synapses')
|
|
|
|
# connect synapses
|
|
i_pre = df_con.loc[:, 'Presynaptic_Index'].values
|
|
i_post = df_con.loc[:, 'Postsynaptic_Index'].values
|
|
syn.connect(i=i_pre, j=i_post)
|
|
|
|
# define connection weight
|
|
syn.w = df_con.loc[:,'Excitatory x Connectivity'].values * params['w_syn']
|
|
|
|
# object to record spikes
|
|
spk_mon = SpikeMonitor(neu)
|
|
|
|
return neu, syn, spk_mon
|
|
|
|
#####################
|
|
# running simulations
|
|
def get_spk_trn(spk_mon):
|
|
'''Extracts spike times from 'spk_mon'
|
|
|
|
The spike times recorded in the SpikeMonitor object during
|
|
simulation are converted to a list of times for each neurons.
|
|
Returns dict with "brian ID": "list of spike times".
|
|
|
|
Parameters
|
|
----------
|
|
spk_mon : SpikeMonitor
|
|
Contains recorded spike times
|
|
|
|
Returns
|
|
-------
|
|
spk_trn : dict
|
|
Mapping between brian neuron IDs and spike times
|
|
'''
|
|
|
|
spk_trn = {k: v for k, v in spk_mon.spike_trains().items() if len(v)}
|
|
|
|
return spk_trn
|
|
|
|
def construct_dataframe(res, exp_name, i2flyid):
|
|
'''Take spike time dict and collects spikes in pandas dataframe
|
|
|
|
Parameters
|
|
----------
|
|
res : list
|
|
List with spike time dicts for each trial
|
|
exp_name : str
|
|
Name of the experiment
|
|
i2flyid : dict
|
|
Mapping between Brian IDs and flywire IDs
|
|
|
|
Returns
|
|
-------
|
|
df : pandas.DataFrame
|
|
Dataframe where each row is one spike
|
|
'''
|
|
|
|
ids, ts, nrun = [], [], []
|
|
|
|
for n, i in enumerate(res):
|
|
for j, k in i.items():
|
|
ids.extend([j for _ in k])
|
|
nrun.extend([n for _ in k])
|
|
ts.extend([float(l) for l in k])
|
|
|
|
d = {
|
|
't': ts,
|
|
'trial': nrun,
|
|
'flywire_id': [i2flyid[i] for i in ids],
|
|
'exp_name': exp_name
|
|
}
|
|
df = pd.DataFrame(d)
|
|
|
|
return df
|
|
|
|
def run_trial(exc, exc2, slnc, path_comp, path_con, params):
|
|
'''Run single trial of coactivation/silencng experiment
|
|
|
|
During the coactivation experiment, the neurons in 'exc' are
|
|
Poisson inputs. The simulation runs for 't_run'.
|
|
|
|
|
|
Parameters
|
|
----------
|
|
exc: list
|
|
contains indices of neurons for PoissonInput
|
|
exc2: list
|
|
contains indices of neurons for PoissonInput of different frequency `r_poi2`
|
|
slnc: list
|
|
contains indices of neurons to silence
|
|
path_comp: Path
|
|
path to "completeness materialization" dataframe
|
|
path_con: Path
|
|
path to "connectivity" dataframe
|
|
params : dict
|
|
Constants and equations that are used to construct the brian2 network model
|
|
|
|
Returns
|
|
-------
|
|
spk_trn : dict
|
|
Mapping between brian neuron IDs and spike times
|
|
'''
|
|
|
|
# get default network
|
|
neu, syn, spk_mon = create_model(path_comp, path_con, params)
|
|
# define Poisson input for excitation
|
|
poi_inp, neu = poi(neu, exc, exc2, params)
|
|
# silence neurons
|
|
syn = silence(slnc, syn)
|
|
# collect in Network object
|
|
net = Network(neu, syn, spk_mon, *poi_inp)
|
|
|
|
# run simulation
|
|
net.run(duration=params['t_run'])
|
|
|
|
# spike times
|
|
spk_trn = get_spk_trn(spk_mon)
|
|
|
|
return spk_trn
|
|
|
|
|
|
def run_exp(exp_name, neu_exc, path_res, path_comp, path_con,
|
|
params=default_params, neu_slnc=[], neu_exc2=[],
|
|
n_proc=-1, force_overwrite=False):
|
|
'''
|
|
Run default network experiment
|
|
Neurons in `neu_exc` are Poisson external inputs
|
|
Neurons in `neu_slnc` are silenced
|
|
|
|
|
|
Parameters
|
|
----------
|
|
exp_name: str
|
|
name of the experiment
|
|
neu_exc: list
|
|
contains custom names or flywire IDs of neurons to be excited
|
|
path_res: str
|
|
path to the output folder where spike data is stored
|
|
path_comp: str
|
|
path to "completeness materialization" dataframe
|
|
path_con: str
|
|
path to "connectivity" dataframe
|
|
params : dict
|
|
Constants and equations that are used to construct the brian2 network model
|
|
neu_slnc : list (optional)
|
|
contains custom names or flywire IDs of neurons to be silenced
|
|
neu_exc2 : list (optional)
|
|
contains custom names or flywire IDs of neurons to be excited with different frequency
|
|
requires `r_poi2` to be defined in `params`
|
|
n_proc: int (optional)
|
|
number of cores to be used for parallel runs
|
|
default: -1 (use all available cores)
|
|
n_proc=1 is equivalent serial code
|
|
force_overwrite : bool (optional)
|
|
If True, overwrite output files, else skip simulation, by default False
|
|
'''
|
|
|
|
# convert to Path objects
|
|
path_res, path_comp, path_con = [ Path(i) for i in [path_res, path_comp, path_con] ]
|
|
|
|
# define output files
|
|
path_save = path_res / '{}.parquet'.format(exp_name)
|
|
if path_save.is_file() and not force_overwrite:
|
|
print('>>> Skipping experiment {} because {} exists and force_overwrite = {}'.format(exp_name, path_save, force_overwrite))
|
|
return
|
|
|
|
# load name/id mappings
|
|
df_comp = pd.read_csv(path_comp, index_col=0) # load completeness dataframe
|
|
flyid2i = {j: i for i, j in enumerate(df_comp.index)} # flywire id: biran ID
|
|
i2flyid = {j: i for i, j in flyid2i.items()} # brian ID: flywire ID
|
|
|
|
# print info
|
|
print('>>> Experiment: {}'.format(exp_name))
|
|
print(' Output file: {}'.format(path_save))
|
|
print(' Excited neurons: {}'.format(len(neu_exc + neu_exc2)))
|
|
if neu_slnc:
|
|
print(' Silenced neurons: {}'.format(len(neu_slnc)))
|
|
|
|
# start parallel calculation
|
|
n_run = params['n_run']
|
|
start = time()
|
|
with parallel_backend('loky', n_jobs=n_proc):
|
|
exc = [ flyid2i[n] for n in neu_exc ]
|
|
exc2 = [ flyid2i[n] for n in neu_exc2 ]
|
|
slnc = [ flyid2i[n] for n in neu_slnc ]
|
|
res = Parallel()(
|
|
delayed(
|
|
run_trial)(exc, exc2, slnc, path_comp, path_con, params) for _ in range(n_run))
|
|
|
|
|
|
# print simulation time
|
|
walltime = time() - start
|
|
print(' Elapsed time: {} s'.format(int(walltime)))
|
|
|
|
# dataframe with spike times
|
|
df = construct_dataframe(res, exp_name, i2flyid)
|
|
|
|
# store spike data
|
|
df.to_parquet(path_save, compression='brotli')
|