Python and NEURON scripts for running network simulations of reduced-morphology layer V pyramidal cells.
Tuomo Maki-Marttunen, 2015-2018
CC BY 4.0
HOC-commands for simulations including in vivo-like synaptic firing based on (Hay & Segev 2015, "Dendritic excitability and gain control
in recurrent cortical microcircuits", Cerebral Cortex
25(10): 3561-3571)
Library of variants based on (Maki-Marttunen et al. 2016, "Functional Effects of Schizophrenia-Linked Genetic Variants on Intrinsic
Single-Neuron Excitability: A Modeling Study", Biol Psychiatry Cogn Neurosci Neuroimaging. 2016 Jan
1;1(1):49-59.
Files included:
CaDynamics_E2.mod #mod file for Intracellular Ca dynamics
Ca_HVA.mod #mod file for HVA Ca currents
Ca_LVAst.mod #mod file for LVA Ca currents
Ih.mod #mod file for HCN current
Im.mod #mod file for M-current
K_Pst.mod #mod file for persistent K current
K_Tst.mod #mod file for transient K current
NaTa_t.mod #mod file for transient Na current
Nap_Et2.mod #mod file for persistent Na current
ProbAMPANMDA2.mod #mod file for probabilistic glutamatergic synapse
ProbAMPANMDA2group.mod #mod file for probabilistic glutamatergic synapse group
ProbAMPANMDA2groupdet.mod #mod file for probabilistic glutamatergic synapse group with pre-determined activation times
ProbUDFsyn2.mod #mod file for probabilistic GABAergic synapse
ProbUDFsyn2group.mod #mod file for probabilistic GABAergic synapse group
ProbUDFsyn2groupdet.mod #mod file for probabilistic GABAergic synapse group with pre-determined activation times
README.html #This file
SK_E2.mod #mod file for SK current
SKv3_1.mod #mod file for Kv3.1 potassium current
approxhaynetstuff.py #Python library for running simulations with reduced-morphology L5PCs
calcEEG.py #Python file for predicting the EEG signal. For spikes*.sav that have already been combined
calcEEG_uncombined.py #Python file for the EEG signal. For spikes*.sav that have not been combined (assumes they were run single-thread)
calcmutgains.py #Python file for running non-oscillatory simulations of single variants
calcmutgains_comb.py #Python file for running non-oscillatory simulations of variant combinations
calcmutoscs.py #Python file for running oscillatory simulations of single variants
calcmutoscs_comb.py #Python file for running oscillatory simulations of variant combinations
calcspectra.py #Python file for calculating the power spectrum of the population spike trains (single variants)
calcspectra_comb.py #Python file for calculating the power spectrum of the population spike trains (variant combinations)
combinemutgains_parallel.py #Python file for combining simulation results (single-variant, non-oscillatory)
combinemutgains_parallel_comb.py #Python file for combining simulation results (variant combination, non-oscillatory)
combinemutgains_parallel_withLFP.py #Python file for combining simulation results (single-variant, non-oscillatory, with recording of dipole moments)
combinemutoscs_parallel.py #Python file for combining simulation results (single-variant, oscillatory)
combinemutoscs_parallel_comb.py #Python file for combining simulation results (variant combination, oscillatory)
drawfig1ab.py #Python file for drawing simulated results (Fig. 1A--B)
drawfig1c.py #Python file for drawing simulated results (Fig. 1C)
drawfig2ab.py #Python file for drawing simulated results (Fig. 2A--B)
drawfig2c.py #Python file for drawing simulated results (Fig. 2C)
drawstationary_EEG.py #Python file for plotting single-cell simulation results with dipole moments and EEG signal (Fig. 4A--C)
drawstationary_EEG_pop.py #Python file for plotting network simulation results with dipole moments and EEG signal (Fig. 4D--F)
mutation_stuff.py #Python library for the effects of SNP-like SCZ-associated variants
mytools.py #Python library, general tools
runmanycellsLFP.py #Python file for running a single cell simulation with recording of dipole moments
runsinglecellLFP.py #Python file for running a network simulation with recording of dipole moments
saveisidists.py #Python file for saving the ISI distributions (needed for making isidist*.sav)
savespikesshufflephases.py #Python file for saving oscillatory Poisson ISIs (needed for making isidist*.sav)
simosc_parallel.py #Python library, main simulation code for oscillatory simulations (single variants)
simosc_parallel_comb_varconn.py #Python library, main simulation code for oscillatory simulations (variant combinations)
simseedburst_func.py #Python library, main simulation code for nonoscillatory simulations (single variants)
simseedburst_func_comb_varconn.py #Python library, main simulation code for nonoscillatory simulations (variant combinations)
simseedburst_func_withLFP.py #Python library, main simulation code for nonoscillatory simulations, including the recording of dipole moments
models #Directory with HOC files
pars_withmids_combfs_final.sav #Data file with parameters for the reduced-morphology L5PCs
scalings_cs.sav #Data file with scaling coefficients for each variant of mutation_stuff.py
To run NEURON simulation of the network model and plot the results (Figure 1), run the following commands:
################ Initialization ################
NUMP=1 #How many CPUs do you want to use in parallel? Change this to speed up simulations (in single thread, one
#network simulation took me approximately 2 hours).
LASTSEED=0 #In this example, only one repetition used in most simulations. Note that in the article normally 5-15
#repetitions used, and at times much larger. Change this to get more accurate predictions.
nrnivmodl #Compile the mechanisms
########## Simulations for Fig 1 A--B ##########
# If your NEURON installation does not support running scripts in the "nrniv -python" manner, try either of the following:
# 1) Replace "nrniv" by the full path to the binary (typically NEURON_INSTALLATION_DIRECTORY/nrn/x86_64/bin/nrniv
# 2) Replace "mpirun -np $NUMP nrniv -python -mpi PYTHON_SCRIPT.py" by a single-CPU run "python PYTHON_SCRIPT.py NULL NULL NULL".
# The three gibberish arguments are necessary as the scripts (calcmutgains*.py and calmutoscs*.py) consider the arguments from
# the 4th argument on. Thus, the first simulation below would be "python calcmutgains.py NULL NULL NULL 0 $irate 0"
rates=( 0.4 0.5 0.6 0.7 0.8 0.9 1.0 1.1 1.2 1.3 1.4 1.5 1.6 ) #Rate coefficient array
seeds=( `seq 1 200` ) #Seed array
for (( irate=0; irate<${#rates[@]}; irate++ )); #Go through the rate array
do
mpirun -np $NUMP nrniv -python -mpi calcmutgains.py 0 $irate 0 #Main simulation (control)
python combinemutgains_parallel.py 0 $irate 0 $NUMP #Combine results (needed if
# many CPUs were used)
echo "rm spikes_parallel_5_11000_0_0.00039_0.0006_${rates[irate]}_1.07_1.07_${seeds[iseed]}_*_of_150.sav" #Remove the uncombined
rm spikes_parallel_5_11000_0_0.00039_0.0006_${rates[irate]}_1.07_1.07_${seeds[iseed]}_*_of_150.sav # data files.
done
mpirun -np $NUMP nrniv -python -mpi calcmutgains.py 70 6 0 #Main simulation (variant)
python combinemutgains_parallel.py 70 6 0 $NUMP #Combine results
echo "rm spikes_parallel_5_11000_70_0.00039_0.0006_1.0_1.07_1.07_${seeds[iseed]}_*_of_150.sav" #Remove the uncombined
rm spikes_parallel_5_11000_70_0.00039_0.0006_1.0_1.07_1.07_${seeds[iseed]}_*_of_150.sav # data files.
python drawfig1ab.py #Plot the results (Fig1A,B)
############ Simulations for Fig 1 C ###########
# Note that simulations in ../haymod have to #
# be performed in advance #
################################################
MS=( 0 125 126 127 128 189 190 191 192 293 294 295 296 313 314 315 316 333 334 335 336 369 370 371 372 421 422 423 424 441 442 443 444 )
MCS=( 0 1 2 3 )
for iseed in `seq 0 $LASTSEED`;
do
#Single-variant simulations:
for (( iimut=0; iimut<${#MS[@]}; iimut++ ));
do
imut=${MS[iimut]}
for (( irate=0; irate<${#rates[@]}; irate++ ));
do
mpirun -np $NUMP nrniv -python -mpi calcmutgains.py $imut $irate $iseed
python combinemutgains_parallel.py $imut $irate $iseed $NUMP
echo "rm spikes_parallel_5_11000_${imut}_0.00039_0.0006_${rates[irate]}_1.07_1.07_${seeds[iseed]}_*_of_150.sav"
rm spikes_parallel_5_11000_${imut}_0.00039_0.0006_${rates[irate]}_1.07_1.07_${seeds[iseed]}_*_of_150.sav
done
done
#Combination simulations:
for (( iimut=0; iimut<${#MCS[@]}; iimut++ ));
do
imut=${MCS[iimut]}
for (( irate=0; irate<${#rates[@]}; irate++ ));
do
mpirun -np $NUMP nrniv -python -mpi calcmutgains_comb.py $imut $irate $iseed
python combinemutgains_parallel_comb.py $imut $irate $iseed $NUMP
echo "rm spikes_parallel_5_11000_comb${imut}_5.0_NsynE10000_NsynI2500_0.00039_0.0006_${rates[irate]}_1.07_1.07_${seeds[iseed]}_*_of_150.sav"
rm spikes_parallel_5_11000_comb${imut}_5.0_NsynE10000_NsynI2500_0.00039_0.0006_${rates[irate]}_1.07_1.07_${seeds[iseed]}_*_of_150.sav
done
done
done
python drawfig1c.py
########## Simulations for Fig 2 A--B ##########
# In order to model the L5PC network response to oscillations, we need the cumulative probability functions
# of intervals between synaptic activation times (these do not anymore follow exponential distribution as
# in the case of stationary Poisson inputs). The values of the cumulative probability functions are pre-saved
# in files isidists_XXX_0.25.sav (0.25 referring to the fact that the lambda term oscillated +-25%) for each
# phase of the oscillation; these files can be downloaded in Open Science Framework (https://osf.io/2k5ub/).
# Alternatively, the event times of oscillatory Poisson processes could be generated on the fly in the main
# simulation script using mytools.oscillatorypoissontimeseries(), but if the same oscillation frequency is
# used many times, it turns out to be more efficient to save the cumulative probability functions of the
# inter-event intervals for each phase of the oscillation and sample from that distribution.
wget https://osf.io/zv29k/download # zipped file size 200 MB
mv download isidists.tar.gz # This is the original name of the file
tar xvfz isidists.tar.gz # NB: This will need an extra 4 GB of space! If you only want to simulate
# panel B, only extract isidists_1.5_0.25.sav, and if you
# want to simulate panel A but not C, choose the isidist files of the 17
# frequencies needed and remove the rest
# Alternatively, these files can be generated from scratch as follows (this is a heavy operation):
FREQARRAY=( 0.5 0.625 0.75 0.875 1.0 1.25 1.5 1.75 2.0 2.5 3.0 3.5 4.0 5.0 7.5 10.0 15.0 )
for (( iosc=0; iosc<${#FREQARRAY[@]}; iosc++ )); #Go through the frequency array
do
for (( iphase=0; iphase<50; iphase++ )); #Go through 50 phases
do #For each of the 50 phases, simulate
python savespikesshufflephases.py $iphase 12000 ${FREQARRAY[iosc]} 0.25 #oscillatory Poisson trains of events.
done #
python saveisidists.py ${FREQARRAY[iosc]} 0.25 #Save the inter-event distributions
done #to files isidists_$freq_0.25.sav
# The simulation of L5PC network activity in the presence of oscillations in the backround synaptic firing
# starts here.
FREQARRAY=( 0.5 0.625 0.75 0.875 1.0 1.25 1.5 1.75 2.0 2.5 3.0 3.5 4.0 5.0 7.5 10.0 15.0 )
seeds=( `seq 1 200` ) #Seed array
for (( iosc=0; iosc<${#FREQARRAY[@]}; iosc++ )); #Go through the frequency array
do
mpirun -np $NUMP nrniv -python -mpi calcmutoscs.py 0 $iosc 0 #Main simulation (control)
python combinemutoscs_parallel.py 0 $iosc 0 $NUMP #Combine results (needed if
# many CPUs were used)
echo "rm spikes_parallel_osc${FREQARRAY[iosc]}_5_11000_0_0.00039_0.0006_1.0_1.07_1.07_${seeds[iseed]}_*_of_150.sav" #Remove the uncombined
rm spikes_parallel_osc${FREQARRAY[iosc]}_5_11000_0_0.00039_0.0006_1.0_1.07_1.07_${seeds[iseed]}_*_of_150.sav # data files.
done
mpirun -np $NUMP nrniv -python -mpi calcmutoscs.py 70 6 0 #Main simulation (variant)
python combinemutoscs_parallel.py 70 6 0 $NUMP #Combine results
echo "rm spikes_parallel_osc1.5_5_11000_70_0.00039_0.0006_1.0_1.07_1.07_${seeds[iseed]}_*_of_150.sav" #Remove the uncombined
rm spikes_parallel_osc1.5_5_11000_70_0.00039_0.0006_1.0_1.07_1.07_${seeds[iseed]}_*_of_150.sav # data files.
for (( iosc=0; iosc<${#FREQARRAY[@]}; iosc++ )); #Go through the frequency array
do
python calcspectra.py 0 ${FREQARRAY[iosc]} #Calculate the power spectrum of the control
#network's population spike train response
#Saves spectrum_freq*_0.sav
done
python calcspectra.py 70 1.5 #Calculate the spectra for variant network
#Saves spectrum_freq1.5_70.sav
python drawfig2ab.py #Plot the results (Fig1A,B)
############ Simulations for Fig 2 C ###########
MS=( 0 125 126 127 128 189 190 191 192 293 294 295 296 313 314 315 316 333 334 335 336 369 370 371 372 421 422 423 424 441 442 443 444 )
MCS=( 0 1 2 3 )
for iseed in `seq 0 $LASTSEED`;
do
#Single-variant simulations:
for (( iimut=0; iimut<${#MS[@]}; iimut++ ));
do
imut=${MS[iimut]}
for (( iosc=0; iosc<${#FREQARRAY[@]}; iosc++ ));
do
mpirun -np $NUMP nrniv -python -mpi calcmutoscs.py $imut $iosc $iseed
python combinemutoscs_parallel.py $imut $iosc $iseed $NUMP
echo "rm spikes_parallel_osc${FREQARRAY[iosc]}_5_11000_${imut}_0.00039_0.0006_1.0_1.07_1.07_${seeds[iseed]}_*_of_150.sav"
rm spikes_parallel_osc${FREQARRAY[iosc]}_5_11000_${imut}_0.00039_0.0006_1.0_1.07_1.07_${seeds[iseed]}_*_of_150.sav
done
done
#Combination simulations:
for (( iimut=0; iimut<${#MCS[@]}; iimut++ ));
do
imut=${MCS[iimut]}
for (( iosc=0; iosc<${#FREQARRAY[@]}; iosc++ ));
do
mpirun -np $NUMP nrniv -python -mpi calcmutoscs_comb.py $imut $iosc $iseed
python combinemutoscs_parallel_comb.py $imut $iosc $iseed $NUMP
echo "rm spikes_parallel_osc${FREQARRAY[iosc]}_5_11000_comb${imut}_5.0_NsynE10000_NsynI2500_0.00039_0.0006_1.0_1.07_1.07_${seeds[iseed]}_*_of_150.sav"
rm spikes_parallel_osc${FREQARRAY[iosc]}_5_11000_comb${imut}_5.0_NsynE10000_NsynI2500_0.00039_0.0006_1.0_1.07_1.07_${seeds[iseed]}_*_of_150.sav
done
done
done
#Calculate the spectra for each variant and each frequency:
for (( iosc=0; iosc<${#FREQARRAY[@]}; iosc++ ));
do
for (( iimut=0; iimut<${#MS[@]}; iimut++ ));
do
python calcspectra.py ${MS[iimut]} ${FREQARRAY[iosc]} #Saves spectrum_freq*_${MS[iimut]}.sav
done
python calcspectra_comb.py ${FREQARRAY[iosc]} #Saves spectrum_freq*_comb0.sav
done
python drawfig2c.py
########## Simulations for Fig 4 A--F ##########
# In order to model the EEG signature of the neuron population, we use a slightly modified simulation script
# (simseeburst_func_withLFP.py). This script uses the functionality of the LFPy toolbox, where the main
# NEURON simulation loop is replaced by our own step-by-step simulation loop that calculates and stores all
# dipole moments. This script does not need LFPy to be installed, but the parts of LFPy that were needed are
# hard-coded in the script. From the output of this script, we calculate the EEG signature of the population.
# This, by contrast, requires that LFPy be installed (see https://lfpy.github.io/)
#
# Note that while the above parts work fine with NEURON 7.4, the part below requires NEURON 7.5 due to the
# use_fast_imem() functionality of the CVode.
python runsinglecellLFP.py # This is a rather light operation as it is only one cell simulated for 11000 ms.
python calcEEG_uncombined.py spikes_parallel_5_11000_0_0.00039_0.0006_1.0_1.07_1.07_1_0_of_1_withLFP.sav
# Calculates the EEG signature using the head model of four concentric spheres.
# Saves to spikes_parallel_5_11000_0_0.00039_0.0006_1.0_1.07_1.07_1_0_of_1_withLFP_EEG.sav
python drawstationary_EEG.py 1 # Draws the results (Fig 4 A--C)
python runmanycellsLFP.py # This takes several hours on an average PC. To speed up use an alternative
# mpirun -np $NUMP nrniv -python -mpi runmanycellsLFP.py (but then the
# NUMP argument in the following line should also be set to $NUMP, i.e.
#python combinemutgains_parallel_withLFP.py 0 6 0 $NUMP
python combinemutgains_parallel_withLFP.py 0 6 0 1
# Saves spikes_parallel150_mutID0_1.0_gNoise1.07_gsyn1.07_seed1_withLFP.sav
python calcEEG.py spikes_parallel150_mutID0_1.0_gNoise1.07_gsyn1.07_seed1_withLFP.sav
# Calculates the EEG signature using the head model of four concentric spheres.
# Saves to spikes_parallel150_mutID0_1.0_gNoise1.07_gsyn1.07_seed1_withLFP_EEG.sav
python drawstationary_EEG_pop.py 1 # Draws the results (Fig 4 D--E)