diff --git a/K_stable.npy b/K_stable.npy deleted file mode 100644 index 8edd2bf..0000000 Binary files a/K_stable.npy and /dev/null differ diff --git a/README.md b/README.md index 9528759..b715d6e 100644 --- a/README.md +++ b/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 diff --git a/multiarea_model/default_params.py b/multiarea_model/default_params.py index a09b449..2a42af8 100644 --- a/multiarea_model/default_params.py +++ b/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} diff --git a/multiarea_model/multiarea_model.py b/multiarea_model/multiarea_model.py index 235246e..4260044 100644 --- a/multiarea_model/multiarea_model.py +++ b/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'] diff --git a/multiarea_model/simulation.py b/multiarea_model/simulation.py index 4d474ba..9c0f5cf 100644 --- a/multiarea_model/simulation.py +++ b/multiarea_model/simulation.py @@ -331,6 +331,7 @@ 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, @@ -338,7 +339,8 @@ def logging(self): '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, diff --git a/run_example_downscaled.py b/run_example_downscaled.py index d5f2781..5f984dc 100644 --- a/run_example_downscaled.py +++ b/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]))) diff --git a/run_example_fullscale.py b/run_example_fullscale.py index 30e0c8b..c8c0036 100644 --- a/run_example_fullscale.py +++ b/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])))