123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277 |
- import numpy as np
- from brian2 import BrianLogger
- from brian2.units import *
- import warnings
- warnings.simplefilter(action='ignore', category=FutureWarning) # Otherwise pypet's FutureWarning spams the console output
- from pypet import Environment, cartesian_product, Trajectory
- 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_network.spatial_network_setup import get_synaptic_weights, \
- create_head_direction_input, ex_in_network
- from scripts.spatial_network.models import *
- from scripts.spatial_maps.supplement_pinwheel_map.pinwheel_map import PinwheelMap
- POLARIZED = 'ellipsoid'
- CIRCULAR = 'circular'
- NO_SYNAPSES = 'no conn'
- DATA_FOLDER = "../../../data/"
- LOG_FOLDER = "../../../logs/"
- TRAJ_NAME = "spatial_network_pinwheel_map"
- TRAJ_NAME_PINWHEEL_MAPS = "precalculated_pinwheel_maps"
- # SCALE_RANGE = [200.0]
- # SEED_RANGE = [1]
- SCALE_RANGE = np.linspace(1.0, 650.0, 24, endpoint=True).tolist()
- SEED_RANGE = range(10)
- def get_uniform_pinwheel_map(scale, seed, sheet_size, N_E, data_folder=None):
- if data_folder is None:
- data_folder = DATA_FOLDER
- traj = Trajectory(filename=data_folder + TRAJ_NAME_PINWHEEL_MAPS + ".hdf5")
- traj.f_load(index=-1, load_parameters=2, load_results=2)
- available_scales = sorted(list(set(traj.f_get("scale").f_get_range())))
- closest_scale = available_scales[np.argmin(np.abs(np.array(available_scales)-scale))]
- if closest_scale != scale:
- print("Warning: desired correlation length {:.1f} not available. Taking {:.1f} instead".format(
- scale, closest_scale))
- map_by_params = lambda x, y: x == scale and y == seed
- idx_iterator = traj.f_find_idx(['scale', 'seed'], map_by_params)
- for idx in idx_iterator:
- traj.v_idx = idx
- map_angle_grid = traj.crun.pinwheel_map
- number_of_excitatory_neurons_per_row = int(np.sqrt(N_E))
- map = PinwheelMap(number_of_excitatory_neurons_per_row, number_of_excitatory_neurons_per_row,
- scale, sheet_size, sheet_size, seed)
- # Uniformize pinwheel map
- nrow = number_of_excitatory_neurons_per_row
- n = map_angle_grid / np.pi
- m = np.concatenate(n)
- sorted_idx = np.argsort(m)
- max_val = nrow * 2
- idx = len(m) // max_val
- for ii, val in enumerate(range(max_val)):
- m[sorted_idx[ii * idx:(ii + 1) * idx]] = val
- p_map = (m - nrow) / nrow
- map.angle_grid = p_map.reshape(nrow, -1) * np.pi
- return map.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
- pinwheel_map = get_uniform_pinwheel_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)), pinwheel_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():
- BrianLogger.suppress_name('method_choice')
- # 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=24, 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")
- scale_range = SCALE_RANGE
- seed_range = SEED_RANGE
- 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()
|