Merge pull request #57 from mlober/reproducibility

Reproducibility

README.md

@@ -225,6 +225,12 @@ The multi-area model can be run in different modes.
    combining this mode with the previous mode 'Subset of the
    network').
 
+5. Dynamical regimes
+
+   As described in Schmidt et al. (2018) (https://doi.org/10.1371/journal.pcbi.1006359),
+   the model can be run in two dynamical regimes - the Ground state and the Metastable state.
+   The state is controlled by the value of the ```cc_weights_factor``` and ```cc_weights_I_factor``` parameters.
+
 ## Test suite
 
 The `tests/` folder holds a test suite that tests different aspects of

multiarea_model/default_params.py

@@ -75,16 +75,6 @@
 Single-neuron parameters
 """
 
-sim_params.update(
-    {
-        'initial_state': {
-            # mean of initial membrane potential (in mV)
-            'V_m_mean': -58.0,
-            # std of initial membrane potential (in mV)
-            'V_m_std': 10.0
-        }
-    })
-
 # dictionary defining single-cell parameters
 single_neuron_dict = {
     # Leak potential of the neurons (in mV).
@@ -109,9 +99,11 @@
     'neuron_model': 'iaf_psc_exp',
     # neuron parameters
     'single_neuron_dict': single_neuron_dict,
-    # Mean and standard deviation for the
-    # distribution of initial membrane potentials
-    'V0_mean': -100.,
+    # The initial membrane potential distribution is chosen 
+    # to be centered on a hyperpolarized state to limit synchronous 
+    # initial firing, which could bring the model into a high-activity 
+    # state when the parameters are set according to the metastable condition.
+    'V0_mean': -150.,
     'V0_sd': 50.}
 
 network_params.update({'neuron_params': neuron_params})
@@ -122,11 +114,12 @@
 """
 connection_params = {
     # Whether to apply the stabilization method of
-    # Schuecker, Schmidt et al. (2017). Default is None.
-    # Options are True to perform the stabilization or
+    # Schuecker, Schmidt et al. (2017).
+    # Options are True to perform the stabilization,
+    # None to not perform the stabilization or
     # a string that specifies the name of a binary
-    # numpy file containing the connectivity matrix
-    'K_stable': None,
+    # numpy file containing the connectivity matrix.
+    'K_stable': os.path.join(base_path+'/figures/SchueckerSchmidt2017', 'K_prime_original.npy'),
 
     # Whether to replace all cortico-cortical connections by stationary
     # Poisson input with population-specific rates (het_poisson_stat)
@@ -158,7 +151,7 @@
     'E_specificity': True,
 
     # Relative inhibitory synaptic strength (in relative units).
-    'g': -16.,
+    'g': -11.,
 
     # compute average indegree in V1 from data
     'av_indegree_V1': np.mean([av_indegree_Cragg, av_indegree_OKusky]),
@@ -169,10 +162,10 @@
     'rho_syn': 'constant',
 
     # Increase the external Poisson indegree onto 5E and 6E
-    'fac_nu_ext_5E': 1.,
-    'fac_nu_ext_6E': 1.,
+    'fac_nu_ext_5E': 1.125,
+    'fac_nu_ext_6E': 1.41666667,
     # to increase the ext. input to 23E and 5E in area TH
-    'fac_nu_ext_TH': 1.,
+    'fac_nu_ext_TH': 1.2,
 
     # synapse weight parameters for current-based neurons
     # excitatory intracortical synaptic weight (mV)
@@ -187,10 +180,14 @@
     'PSC_rel_sd_lognormal': 3.0,
 
     # scaling factor for cortico-cortical connections (chi)
-    'cc_weights_factor': 1.,
-    # factor to scale cortico-cortical inh. weights in relation
-    # to exc. weights (chi_I)
-    'cc_weights_I_factor': 1.,
+    # Default value is 1.9 which reproduces the "Metastable" state
+    # activity described in Schmidt et al. (2018). 
+    # A weight factor of 1.0 produces Ground state activity.
+    'cc_weights_factor': 1.9,
+    # Factor to scale cortico-cortical inh. weights in relation
+    # to exc. weights (chi_I). Default is 2.0 to reproduce metastable state
+    # activity. For ground state activity, set to 1.0.
+    'cc_weights_I_factor': 2.,
 
     # 'switch whether to distribute weights lognormally
     'lognormal_weights': False,
@@ -290,7 +287,7 @@
                  # The simulation time of the mean-field theory integration
                  'T': 50.,
                  # The time step of the mean-field theory integration
-                 'dt': 0.01,
+                 'dt': 0.1,
                  # Time interval for recording the trajectory of the mean-field calcuation
                  # If None, then the interval is set to dt
                  'rec_interval': None}

multiarea_model/multiarea_model.py

@@ -120,8 +120,10 @@ def __init__(self, network_spec, theory=False, simulation=False,
             dat = json.load(f)
 
         self.structure = OrderedDict()
+        self.structure_reversed = OrderedDict()
         for area in dat['area_list']:
             self.structure[area] = dat['structure'][area]
+            self.structure_reversed[area] = self.structure[area][::-1]
         self.N = dat['neuron_numbers']
         self.synapses = dat['synapses']
         self.W = dat['synapse_weights_mean']

multiarea_model/simulation.py

@@ -331,14 +331,16 @@ def logging(self):
         Write runtime and memory for the first 30 MPI processes
         to file.
         """
+        local_spike_counter = nest.GetKernelStatus('local_spike_counter')
         if nest.Rank() < 30:
             d = {'time_prepare': self.time_prepare,
                  'time_network_local': self.time_network_local,
                  'time_network_global': self.time_network_global,
                  'time_simulate': self.time_simulate,
                  'base_memory': self.base_memory,
                  'network_memory': self.network_memory,
-                 'total_memory': self.total_memory}
+                 'total_memory': self.total_memory, 
+                 'local_spike_counter': local_spike_counter}
             fn = os.path.join(self.data_dir,
                               'recordings',
                               '_'.join((self.label,

run_example_downscaled.py

@@ -14,33 +14,20 @@
 """
 d = {}
 conn_params = {'replace_non_simulated_areas': 'het_poisson_stat',
-               'g': -11.,
-               'K_stable': 'K_stable.npy',
-               'fac_nu_ext_TH': 1.2,
-               'fac_nu_ext_5E': 1.125,
-               'fac_nu_ext_6E': 1.41666667,
-               'av_indegree_V1': 3950.}
-input_params = {'rate_ext': 10.}
-neuron_params = {'V0_mean': -150.,
-                 'V0_sd': 50.}
+               'cc_weights_factor': 1.0, # run model in Ground State
+               'cc_weights_I_factor': 1.0}
 network_params = {'N_scaling': 0.01,
                   'K_scaling': 0.01,
-                  'fullscale_rates': os.path.join(base_path, 'tests/fullscale_rates.json'),
-                  'input_params': input_params,
-                  'connection_params': conn_params,
-                  'neuron_params': neuron_params}
+                  'fullscale_rates': os.path.join(base_path, 'tests/fullscale_rates.json')}
 
 sim_params = {'t_sim': 2000.,
               'num_processes': 1,
-              'local_num_threads': 1,
-              'recording_dict': {'record_vm': False}}
-
-theory_params = {'dt': 0.1}
+              'local_num_threads': 1}
 
 M = MultiAreaModel(network_params, simulation=True,
                    sim_spec=sim_params,
-                   theory=True,
-                   theory_spec=theory_params)
+                   theory=True)
+
 p, r = M.theory.integrate_siegert()
 print("Mean-field theory predicts an average "
       "rate of {0:.3f} spikes/s across all populations.".format(np.mean(r[:, -1])))

run_example_fullscale.py

@@ -20,32 +20,17 @@
 resources, for instance on a compute cluster.
 """
 d = {}
-conn_params = {'g': -11.,
-               'K_stable': os.path.join(base_path, 'K_stable.npy'),
-               'fac_nu_ext_TH': 1.2,
-               'fac_nu_ext_5E': 1.125,
-               'fac_nu_ext_6E': 1.41666667,
-               'av_indegree_V1': 3950.}
-input_params = {'rate_ext': 10.}
-neuron_params = {'V0_mean': -150.,
-                 'V0_sd': 50.}
+
 network_params = {'N_scaling': 1.,
-                  'K_scaling': 1.,
-                  'connection_params': conn_params,
-                  'input_params': input_params,
-                  'neuron_params': neuron_params}
+                  'K_scaling': 1.}
 
 sim_params = {'t_sim': 2000.,
               'num_processes': 720,
-              'local_num_threads': 1,
-              'recording_dict': {'record_vm': False}}
-
-theory_params = {'dt': 0.1}
+              'local_num_threads': 1}
 
 M = MultiAreaModel(network_params, simulation=True,
                    sim_spec=sim_params,
-                   theory=True,
-                   theory_spec=theory_params)
+                   theory=True)
 p, r = M.theory.integrate_siegert()
 print("Mean-field theory predicts an average "
       "rate of {0:.3f} spikes/s across all populations.".format(np.mean(r[:, -1])))