import numpy as np from brian2.units import * from pypet import Environment, cartesian_product from pypet.brian2.parameter import Brian2MonitorResult from scripts.spatial_maps.spatial_network_layout import create_grid_of_excitatory_neurons, \ create_interneuron_sheet_entropy_max_orientation, get_excitatory_neurons_in_inhibitory_axonal_clouds from scripts.spatial_maps.perlin_map.uniform_perlin_map import UniformPerlinMap from scripts.spatial_network.spatial_network_setup import get_synaptic_weights, \ create_head_direction_input, ex_in_network from scripts.spatial_network.models import * POLARIZED = 'ellipsoid' CIRCULAR = 'circular' NO_SYNAPSES = 'no conn' DATA_FOLDER = "../../../data/" LOG_FOLDER = "../../../logs/" TRAJ_NAME = "spatial_network_perlin_map" def get_perlin_map(scale, seed, sheet_size, N_E): number_of_excitatory_neurons_per_row = int(np.sqrt(N_E)) # TODO: The +1 should not be there input_map = UniformPerlinMap(number_of_excitatory_neurons_per_row + 1, number_of_excitatory_neurons_per_row + 1, scale, sheet_size, sheet_size, seed) return input_map.get_tuning def get_input_head_directions(traj): directions = np.linspace(-np.pi, np.pi, traj.input.number_of_directions, endpoint=False) return directions def spatial_network_with_entropy_maximisation(traj): sheet_size = traj.input_map.sheet_size N_E = traj.network.N_E N_I = traj.network.N_I perlin_map = get_perlin_map(traj.input_map.scale, traj.input_map.seed, sheet_size, N_E) ex_positions, ex_tunings = create_grid_of_excitatory_neurons(sheet_size, int(np.sqrt(N_E)), perlin_map) inhibitory_axon_long_axis = traj.morphology.long_axis inhibitory_axon_short_axis = traj.morphology.short_axis entropy_maximisation_trial_orientations = traj.interneuron.entropy_maximisation.trial_orientations if \ inhibitory_axon_long_axis != inhibitory_axon_short_axis else 0 inhibitory_axonal_clouds, ellipse_single_trial_entropy = create_interneuron_sheet_entropy_max_orientation( ex_positions, ex_tunings, N_I, inhibitory_axon_long_axis, inhibitory_axon_short_axis, sheet_size, sheet_size, trial_orientations=entropy_maximisation_trial_orientations) ie_connections = get_excitatory_neurons_in_inhibitory_axonal_clouds(ex_positions, inhibitory_axonal_clouds) inhibitory_synapse_strength = traj.synapse.inhibitory * nS excitatory_synapse_strength = traj.synapse.excitatory * mV if inhibitory_synapse_strength != 0.0 * nS and excitatory_synapse_strength != 0.0 * mV \ and inhibitory_axon_long_axis == inhibitory_axon_short_axis: traj.f_add_derived_parameter("morphology.morph_label", CIRCULAR, comment="Interneuron morphology of this run is circular") elif inhibitory_synapse_strength != 0.0 * nS and excitatory_synapse_strength != 0.0 * mV: traj.f_add_derived_parameter("morphology.morph_label", POLARIZED, comment="Interneuron morphology of this run is ellipsoid") else: traj.f_add_derived_parameter("morphology.morph_label", NO_SYNAPSES, comment="There are no interneurons") ex_in_weights, in_ex_weights = get_synaptic_weights(N_E, N_I, ie_connections, excitatory_synapse_strength, inhibitory_synapse_strength) input_directions = get_input_head_directions(traj) for idx, direction in enumerate(input_directions): # We recreate the network here for every dir, which slows down the simulation quite considerably. Otherwise, # we get a problem with saving and restoring the spike times (0s spike for neuron 0) net = ex_in_network(N_E, N_I, hodgkin_huxley_eqs_with_synaptic_conductance, hodgkin_huxley_params, lif_interneuron_eqs, lif_interneuron_params, lif_interneuron_options, delta_synapse_model, delta_synapse_on_pre, delta_synapse_param, ex_in_weights, exponential_synapse, exponential_synapse_on_pre, exponential_synapse_params, in_ex_weights, random_seed=traj.input_map.seed) input_to_excitatory_population = create_head_direction_input(traj.input.baseline * nA, ex_tunings, traj.input.phase_dispersion, traj.input.amplitude * nA, direction) excitatory_neurons = net["excitatory_neurons"] excitatory_neurons.I = input_to_excitatory_population inhibitory_neurons = net["interneurons"] inhibitory_neurons.u_ext = traj.inh_input.baseline * mV inhibitory_neurons.tau = traj.interneuron.tau * ms net.run(traj.simulation.duration * ms) direction_id = 'dir{:d}'.format(idx) traj.f_add_result(Brian2MonitorResult, '{:s}.spikes.e'.format(direction_id), net["excitatory_spike_monitor"], comment='The spiketimes of the excitatory population') traj.f_add_result(Brian2MonitorResult, '{:s}.spikes.i'.format(direction_id), net["inhibitory_spike_monitor"], comment='The spiketimes of the inhibitory population') traj.f_add_result('ex_positions', np.array(ex_positions), comment='The positions of the excitatory neurons on the sheet') traj.f_add_result('ex_tunings', np.array(ex_tunings), comment='The input tunings of the excitatory neurons') ie_connections_save_array = np.zeros((N_I, N_E)) for i_idx, ie_conn in enumerate(ie_connections): for e_idx in ie_conn: ie_connections_save_array[i_idx, e_idx] = 1 traj.f_add_result('ie_adjacency', ie_connections_save_array, comment='Recurrent connection adjacency matrix') axon_cloud_save_list = [[p.x, p.y, p.phi] for p in inhibitory_axonal_clouds] axon_cloud_save_array = np.array(axon_cloud_save_list) traj.f_add_result('inhibitory_axonal_cloud_array', axon_cloud_save_array, comment='The inhibitory axonal clouds') return 1 def main(): # TODO: Set ncores to the desired number of processes to use by pypet env = Environment(trajectory=TRAJ_NAME, comment="Compare the head direction tuning for circular and polarized interneuron morphology, " "when tuning orientations to maximise entropy of connected excitatory tunings.", multiproc=True, filename=DATA_FOLDER, ncores=3, overwrite_file=True, log_folder=LOG_FOLDER) traj = env.trajectory traj.f_add_parameter_group("input_map") traj.f_add_parameter("input_map.scale", 200.0, comment="Scaling factor of the input map") traj.f_add_parameter("input_map.seed", 1, comment="Random seed for input map generation.") traj.f_add_parameter("input_map.sheet_size", 900, comment="Sheet size in um") traj.f_add_parameter_group("network") traj.f_add_parameter("network.N_E", 3600, comment="Number of excitatory neurons") traj.f_add_parameter("network.N_I", 400, comment="Number of inhibitory neurons") traj.f_add_parameter_group("interneuron") traj.f_add_parameter("interneuron.tau", 7., comment="Interneuron timescale in ms") traj.f_add_parameter("interneuron.entropy_maximisation.trial_orientations", 30, comment="Steps for entropy maximisation") traj.f_add_parameter_group("synapse") traj.f_add_parameter("synapse.inhibitory", 30.0, "Strength of conductance-based inhibitory synapse in nS.") traj.f_add_parameter("synapse.excitatory", 2.5, "Strength of conductance-based inhibitory synapse in mV.") traj.f_add_parameter_group("input") traj.f_add_parameter("input.phase_dispersion", 2.5, comment="Standard deviation of incoming head direction input.") traj.f_add_parameter("input.baseline", 0.05, comment="Head direction input baseline") traj.f_add_parameter("input.amplitude", 0.6, comment="Head direction input amplitude") traj.f_add_parameter("input.number_of_directions", 12, comment="Number of probed directions") traj.f_add_parameter_group("inh_input") traj.f_add_parameter("inh_input.baseline", -50., comment="Head direction input baseline") traj.f_add_parameter("inh_input.amplitude", 0., comment="Head direction input amplitude") traj.f_add_parameter_group("morphology") traj.f_add_parameter("morphology.long_axis", 100.0, comment="Long axis of axon ellipsoid") traj.f_add_parameter("morphology.short_axis", 25.0, comment="Short axis of axon ellipsoid") traj.f_add_parameter_group("simulation") traj.f_add_parameter("simulation.dt", 0.1, comment="Network simulation time step in ms") traj.f_add_parameter("simulation.duration", 1000, comment="Network simulation duration in ms") # To test if the simulation, analysis and plotting work as intended, use this setup to have a lower overall runtime # scale_range = [200.0] # seed_range = [1] # TODO: This is the parameter range used for the results shown in Fig 4. Uncomment for the full run. scale_range = np.linspace(1.0, 800.0, 24, endpoint=True).tolist() seed_range = range(10) ellipsoid_parameter_exploration = { "morphology.long_axis": [100.0], "morphology.short_axis": [25.0], "input_map.scale": scale_range, "input_map.seed": seed_range, "synapse.inhibitory": [30.], "synapse.excitatory": [2.5] } corresponding_circular_radius = float(np.sqrt(ellipsoid_parameter_exploration[ "morphology.long_axis"][0] * ellipsoid_parameter_exploration[ "morphology.short_axis"][0])) circle_parameter_exploration = { "morphology.long_axis": [corresponding_circular_radius], "morphology.short_axis": [corresponding_circular_radius], "input_map.scale": ellipsoid_parameter_exploration["input_map.scale"], "input_map.seed": ellipsoid_parameter_exploration["input_map.seed"], "synapse.inhibitory": ellipsoid_parameter_exploration["synapse.inhibitory"], "synapse.excitatory": ellipsoid_parameter_exploration["synapse.excitatory"] } no_conn_parameter_exploration = { "morphology.long_axis": [corresponding_circular_radius], "morphology.short_axis": [corresponding_circular_radius], "input_map.scale": ellipsoid_parameter_exploration["input_map.scale"], "input_map.seed": ellipsoid_parameter_exploration["input_map.seed"], "synapse.inhibitory": [0.], "synapse.excitatory": [0.] } expanded_dicts = [cartesian_product(dict) for dict in [ellipsoid_parameter_exploration, circle_parameter_exploration, no_conn_parameter_exploration]] final_dict = {} for key in expanded_dicts[0].keys(): list_of_parameter_lists = [dict[key] for dict in expanded_dicts] final_dict[key] = sum(list_of_parameter_lists, []) traj.f_explore(final_dict) env.run(spatial_network_with_entropy_maximisation) env.disable_logging() if __name__ == "__main__": main()