Add documentation through project

Author: Daniel Incicau

README.md

@@ -1 +1,20 @@
-# Spiking Neural Network for the Visual Cortex (V1) Layer 5
+Spiking Neural Network for the Visual Cortex (V1) Layer 5
+============================
+
+### Project content
+
+    .
+    ├── data                                        # Contains hypertuning results from amplitude simulation run
+    ├── output                                      # Contains output results & plots from complete simulation run (.gitignored)
+    ├── notebooks                                   # Jupyter notebooks used for testing different scenarios
+    │   ├── layer5_CC_CS_connection.ipynb           # Scalar plot simulation for weighted CC->CS connection probability
+    │   ├── layer5_SST_Soma_selectivity.ipynb       # Scalar plot simulation for weighted SST->CC/CS Soma connection probability
+    │   ├── layer5_sandbox.ipynb                    # Sandbox notebook used when creating initial network topology
+    │   └── ...               
+    ├── layer5_CC_CS.py                             # Main network file that defines the topology of the simulation                   
+    ├── equations.py                                # Equations for different neuron types
+    ├── parameters.py                               # Default parameters for simulation
+    ├── helpers.py                                  # Helper methods for analysis
+    ├── plotting.py                                 # Helper methods for plotting
+    ├── run_amplitude_hypertuing.py                 # Script for running multiple simulations for amplitude hypertuning for different neurons
+    └── run_complete_simulation.py                  # Complete simulation run entrypoint

helpers.py

@@ -30,6 +30,10 @@ def bin(spiketime, dt):
 
 
 def count_spikes_for_neuron_type(spike_mon, from_t=None, to_t=None):
+    """
+    Computes spike count by neuron for all neurons monitored by a spike monitor between given interval of time
+    """
+
     spike_counts = {}
 
     for index in spike_mon.spike_trains():
@@ -45,13 +49,21 @@ def count_spikes_for_neuron_type(spike_mon, from_t=None, to_t=None):
 
 @check_units(sim_duration=ms)
 def compute_firing_rate_for_neuron_type(spike_mon, from_t, to_t):
+    """
+    Computes firing rate by neuron for all neurons monitored by a spike monitor between given interval of time
+    """
+
     spikes_for_i = count_spikes_for_neuron_type(spike_mon, from_t, to_t)
     duration = to_t - from_t
 
     return spikes_for_i / duration
 
 
 def compute_interspike_intervals(spike_mon, from_t, to_t):
+    """
+    Computes inter-spike intervals by neuron for all neurons monitored by a spike monitor between given interval of time
+    """
+
     by_neuron = []
 
     for neuron_index in spike_mon.spike_trains():
@@ -65,6 +77,10 @@ def compute_interspike_intervals(spike_mon, from_t, to_t):
 
 
 def compute_autocorr(spike_intervals):
+    """
+    Computes auto-correlation function on the inter-spike intervals
+    """
+
     if len(spike_intervals) == 0:
         return None, None
 
@@ -75,6 +91,10 @@ def compute_autocorr(spike_intervals):
 
 
 def find_minimum_autocorr(acorr):
+    """
+    Finds the minimum of the auto-correlation function given by acorr
+    """
+
     if acorr is None:
         return None
 
@@ -93,6 +113,10 @@ def find_minimum_autocorr(acorr):
 
 
 def compute_autocorr_struct(interspike_intervals, no_bins):
+    """
+    Helper method for computing all relevant auto-correlation metrics for inter-spike intervals given
+    """
+
     autocorr_sst = None
 
     sorted_isi = np.sort(interspike_intervals)
@@ -110,6 +134,12 @@ def compute_autocorr_struct(interspike_intervals, no_bins):
 
 
 def compute_equilibrium_for_neuron_type(spike_mon, time_frame=0.05, tol=0.001):
+    """
+    Computes equilibrium time for a group of neurons monitored by a spike monitor.
+    It does that by using a sliding window algorithm of `time_frame` and
+    checking if firing rate of neurons has changed between 2 consecutive time windows
+    """
+
     t = 0
     last_spike_t = spike_mon.t[-1] / second if len(spike_mon.t) > 0 else t
 
@@ -128,6 +158,12 @@ def compute_equilibrium_for_neuron_type(spike_mon, time_frame=0.05, tol=0.001):
 
 
 def compute_burst_mask(spikes, maxISI):
+    """
+    Computes a burst mask for the spikes of a given neuron.
+    A spike is said to be in a burst if the inter-spike interval time is less than maxISI
+    Returns a binary vector. 1 means spike is in a burst. 0 means is not
+    """
+
     isi = np.diff(spikes)
     mask = np.zeros(len(spikes))
 
@@ -140,6 +176,10 @@ def compute_burst_mask(spikes, maxISI):
 
 
 def compute_burst_trains(spike_mon, maxISI):
+    """
+    Computes a burst mask for each neuron monitored by the spike_mon
+    """
+
     burst_trains = {}
     for neuron_index in spike_mon.spike_trains():
         burst_mask = compute_burst_mask(spike_mon.spike_trains()[neuron_index], maxISI)
@@ -149,7 +189,11 @@ def compute_burst_trains(spike_mon, maxISI):
 
 
 def compute_burst_lengths(burst_mask):
-    "count how many series of consecutive occurances of 1 appear"
+    """
+    Given a burst mask which is iterated, a burst is defined as series of consecutive occurrences of 1.
+    Each of these series is considered a burst and has an associated burst length.
+    Returns a vector of bursts lengths.
+    """
 
     burst_lengths = []
     idx = 0
@@ -180,6 +224,10 @@ def compute_burst_lengths_by_neuron_group(burst_trains):
 
 
 def save_agg_results_to_folder(results, output_folder=None, file_name='results.json'):
+    """
+    Helper method for saving aggregated results for multiple simulations with different degree inputs to folder
+    """
+
     if output_folder is not None:
         if not os.path.exists(output_folder):
             os.makedirs(output_folder)
@@ -190,6 +238,10 @@ def save_agg_results_to_folder(results, output_folder=None, file_name='results.j
 
 
 def save_results_to_folder(results, output_folder=None, file_name='results.json'):
+    """
+    Helper method for saving individual simulation results to folder
+    """
+
     if output_folder is not None:
         if not os.path.exists(output_folder):
             os.makedirs(output_folder)

layer5_CC_CS.py

@@ -12,11 +12,16 @@ def __init__(self, **entries):
 
 
 def analyse_network_simulation(spike_monitors, state_monitors, synapses, p, output_folder=None):
+    """
+    Does an analysis of a simulation post-run. Saves plots to `output_folder`
+    Returns the results of simulation analysis
+    """
+
     spike_mon_sst, spike_mon_pv, spike_mon_cs, spike_mon_cc = spike_monitors
     state_mon_sst, state_mon_pv, state_mon_cs, state_mon_cc = state_monitors
 
     ################################################################################
-    # Analysis and plotting
+    # Compute equilibrium time of simulation
     ################################################################################
 
     if p.recompute_equilibrium:
@@ -41,6 +46,10 @@ def analyse_network_simulation(spike_monitors, state_monitors, synapses, p, outp
     from_t = equilibrium_t
     to_t = p.duration
 
+    ################################################################################
+    # Analysis and plotting
+    ################################################################################
+
     raster_from_t = from_t if p.plot_only_from_equilibrium else 0
     raster_to_t = min(raster_from_t + 3 * second, p.duration)
     plot_raster(spike_mon_cs, spike_mon_cc, spike_mon_sst, spike_mon_pv, raster_from_t, raster_to_t,
@@ -140,6 +149,11 @@ def analyse_network_simulation(spike_monitors, state_monitors, synapses, p, outp
 
 
 def run_simulation_without_exh_dendrite(network, neurons, synapses, p, base_output_folder, use_synaptic_probabilities):
+    """
+    Runs simulation for network topology with PYR cells having ONLY somas and NO dendrites.
+    Also analyses simulation run and returns results.
+    """
+
     sst_neurons, pv_neurons, cs_neurons, cc_neurons = neurons
 
     E_l = p.E_l  # leak reversal potential
@@ -218,6 +232,13 @@ def run_simulation_without_exh_dendrite(network, neurons, synapses, p, base_outp
 
 
 def run_simulation_for_weighted_sst_soma_dendrite(network, neurons, synapses, p, p_SST_CS_soma, p_SST_CC_soma, base_output_folder, use_synaptic_probabilities):
+    """
+    Runs simulation for network topology with PYR cells having BOTH somas and dendrites.
+    Connection probability of SST->CC/CS Soma is weighted through parameters `p_SST_CS_soma` and `p_SST_CC_soma`
+    Connection probability of SST->CC/CS Dendrite is given by `1 - p_SST_CS_soma` and `1 - p_SST_CC_soma`
+    Also analyses simulation run and returns results.
+    """
+
     sst_neurons, pv_neurons, cs_neurons, cc_neurons = neurons
 
     E_l = p.E_l  # leak reversal potential
@@ -315,13 +336,20 @@ def run_simulation_for_weighted_sst_soma_dendrite(network, neurons, synapses, p,
 
 
 def run_simulation_for_input(params, use_synaptic_probabilities, use_dendrite_model, seed_val, base_output_folder=None):
+    """
+    Given external input and network parameters (through `params`) simulations are run accordingly.
+    If `use_dendrite_model` multiple simulations are run for each (pSST_CS_soma, pSST_CC_soma) pairs
+    If NOT `use_dendrite_model` a single simulation is run.
+    Returns vector of results of simulations
+    """
+
     p = Struct(**params)
 
     start_scope()
     seed(seed_val)
 
     E_l = p.E_l  # leak reversal potential
-    V_t = p.V_t  # spiking threashold
+    V_t = p.V_t  # spiking threshold
 
     assert len(p.pSST_CS_soma) == len(p.pSST_CC_soma)  # since they are taken in pairs
 
@@ -410,7 +438,7 @@ def run_simulation_for_input(params, use_synaptic_probabilities, use_dendrite_mo
     neurons = [sst_neurons, pv_neurons, cs_neurons, cc_neurons]
 
     # ##############################################################################
-    # Define Synapses
+    # Define Synapses (common synapses for simulation)
     # ##############################################################################
 
     synapses = {}
@@ -510,6 +538,10 @@ def run_simulation_for_input(params, use_synaptic_probabilities, use_dendrite_mo
 
     defaultclock.dt = p.sim_dt
 
+    # ##############################################################################
+    # Continue defining specific simulation run
+    # ##############################################################################
+
     results = []
     if use_dendrite_model:
         for p_SST_CS_soma, p_SST_CC_soma in zip(p.pSST_CS_soma, p.pSST_CC_soma):
@@ -532,6 +564,12 @@ def run_simulation_for_input(params, use_synaptic_probabilities, use_dendrite_mo
 
 
 def run_complete_simulation(params, use_dendrite_model=True, use_synaptic_probabilities=True, seed_val=12345):
+    """
+    Runs a complete simulation for given parameters.
+    Computes both individual simulation results and aggregated simulation results.
+    Selectivity of network is calculated through varying external input to neurons.
+    """
+
     p = Struct(**params)
 
     N = [p.N_cs, p.N_cc, p.N_sst, p.N_pv]

parameters.py

@@ -32,7 +32,7 @@
     "E_e":    0 * mV,  # Excitatory synaptic reversal potential
     "E_i":  -80 * mV,  # Inhibitory synaptic reversal potential
 
-    "V_t":  -50 * mV,  # spiking threashold
+    "V_t":  -50 * mV,  # spiking threshold
     "V_r":  -70 * mV,  # reset potential ~ same as E_l ~
 
     "c_d": 2600 * pA,  # back-propagates somatic spikes to the dendrites

input_amplitude_test.py run_amplitude_hypertuning.pyinput_amplitude_test.py renamed to run_amplitude_hypertuning.py

@@ -64,6 +64,8 @@ def run_input_amplitude_simulation(params, input_amplitudes, seed_val=12345):
 
 
 def test_input_amplitudes(params, csv_writer=None, seed_val=12345):
+
+    # combinatorically run simulations for the following neuron type amplitudes
     input_cs_amplitudes = [100]
     input_cc_amplitudes = [100, 200, 300]
     input_sst_amplitudes = [10, 25, 50, 75, 100]