123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193 |
- import numpy as np
- from brian2.units import *
- from pypet import Environment
- import matplotlib.pyplot as plt
- from scripts.interneuron_placement import create_grid_of_excitatory_neurons, \
- create_interneuron_sheet_entropy_max_orientation, get_excitatory_neurons_in_inhibitory_axonal_clouds, Pickle, \
- get_position_mesh
- from scripts.ring_network.head_direction import ex_in_network
- from scripts.spatial_network.head_direction_index_over_noise_scale import excitatory_eqs, excitatory_params, \
- lif_interneuron_eqs, lif_interneuron_params, lif_interneuron_options, ei_synapse_model, ei_synapse_on_pre, \
- ei_synapse_param, ie_synapse_model, ie_synapse_on_pre, ie_synapse_param, get_synaptic_weights, \
- create_head_direction_input, calculate_rates
- from scripts.spatial_network.run_entropy_maximisation_orientation_map import get_orientation_map
- DATA_FOLDER = "../../data/"
- LOG_FOLDER = "../../logs/"
- TRAJ_NAME = "hdi_maximisation_by_differential_evolution"
- def plot_and_save_sheet(ex_positions, ex_tunings, interneuron_positions, interneuron_orientations, mean_hdi, a, b, label):
- axonal_clouds = [Pickle(x, y, a, b, phi)
- for x, y, phi in zip(*zip(*interneuron_positions), interneuron_orientations)]
- X, Y = get_position_mesh(ex_positions)
- n_ex = int(np.sqrt(len(ex_positions)))
- head_dir_preference = np.array(ex_tunings).reshape((n_ex, n_ex))
- fig = plt.figure()
- ax = fig.add_subplot(111)
- plt.set_cmap('twilight')
- c = ax.pcolor(X, Y, head_dir_preference.T, vmin=-np.pi, vmax=np.pi)
- fig.colorbar(c, ax=ax, label="Tuning")
- if axonal_clouds is not None:
- for i, p in enumerate(axonal_clouds):
- ell = p.get_ellipse()
- ax.add_artist(ell)
- plt.title('Mean HDI: {}'.format(mean_hdi))
- plt.savefig(DATA_FOLDER + 'optimization_figure_' + label + '.png')
- def optimize_interneuron_orientation_by_hdi(traj):
- sheet_size = traj.orientation_map.sheet_size
- N_E = traj.network.N_E
- N_I = traj.network.N_I
- orientation_map = get_orientation_map(traj.orientation_map.correlation_length, traj.orientation_map.seed,
- sheet_size, N_E)
- ex_positions, ex_tunings = create_grid_of_excitatory_neurons(sheet_size,
- sheet_size,
- int(np.sqrt(N_E)), orientation_map)
- inhibitory_axon_long_axis = traj.morphology.long_axis
- inhibitory_axon_short_axis = traj.morphology.short_axis
- # Einmal die Entropy Maximierung
- inhibitory_axonal_clouds, _ = 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=30)
- initial_interneuron_orientations = [p.phi for p in inhibitory_axonal_clouds]
- interneuron_positions = [(p.x, p.y) for p in inhibitory_axonal_clouds]
- interneuron_orientation_bounds = [(-np.pi,np.pi) for x in interneuron_positions]
- mean_hdi = 1. - spatial_network_hdi_by_orientations(initial_interneuron_orientations, interneuron_positions, ex_positions, ex_tunings, traj)
- plot_and_save_sheet(ex_positions, ex_tunings, interneuron_positions, initial_interneuron_orientations, mean_hdi, inhibitory_axon_long_axis, inhibitory_axon_short_axis, 'before')
- # Und dann die optimierte Variante
- axon_cloud_save_array = np.load(DATA_FOLDER + 'current_cloud_save_array.npy')
- interneuron_positions = [(p[0], p[1]) for p in axon_cloud_save_array]
- initial_interneuron_orientations = [p[2] for p in axon_cloud_save_array]
- mean_hdi = 1. - spatial_network_hdi_by_orientations(initial_interneuron_orientations, interneuron_positions, ex_positions, ex_tunings, traj)
- plot_and_save_sheet(ex_positions, ex_tunings, interneuron_positions, initial_interneuron_orientations, mean_hdi, inhibitory_axon_long_axis, inhibitory_axon_short_axis, 'before')
- plt.show()
- def spatial_network_hdi_by_orientations(interneuron_orientations, interneuron_positions, ex_positions, ex_tunings, traj):
- N_E = traj.network.N_E
- N_I = traj.network.N_I
- inhibitory_axon_long_axis = traj.morphology.long_axis
- inhibitory_axon_short_axis = traj.morphology.short_axis
- inhibitory_axonal_clouds = [Pickle(x, y, inhibitory_axon_long_axis, inhibitory_axon_short_axis, phi)
- for x, y, phi in zip(*zip(*interneuron_positions), interneuron_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
- ex_in_weights, in_ex_weights = get_synaptic_weights(N_E, N_I, ie_connections, excitatory_synapse_strength,
- inhibitory_synapse_strength)
- sharpness = 1.0 / (traj.input.width) ** 2
- directions = np.linspace(-np.pi, np.pi, traj.input.number_of_directions, endpoint=False)
- firing_rate_array = np.ndarray((traj.N_E, traj.input.number_of_directions))
- net = ex_in_network(N_E, N_I, excitatory_eqs, excitatory_params, lif_interneuron_eqs,
- lif_interneuron_params,
- lif_interneuron_options, ei_synapse_model, ei_synapse_on_pre,
- ei_synapse_param,
- ex_in_weights, ie_synapse_model, ie_synapse_on_pre,
- ie_synapse_param, in_ex_weights, random_seed=2)
- for dir_idx, dir in enumerate(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)
- input_to_excitatory_population = create_head_direction_input(traj.input.baseline * nA, ex_tunings,
- sharpness,
- traj.input.amplitude * nA, dir)
- excitatory_neurons = net["excitatory_neurons"]
- excitatory_neurons.I = input_to_excitatory_population
- net.run(traj.simulation.duration * ms)
- exc_spike_mon = net["excitatory_spike_monitor"]
- ex_spike_trains = exc_spike_mon.spike_trains()
- ex_spike_rates = calculate_rates(ex_spike_trains.values())
- for n_idx, spike_rate in enumerate(ex_spike_rates):
- firing_rate_array[n_idx, dir_idx] = spike_rate
- #TODO: Get head direction indices
- tuning_vectors = np.zeros((N_E, traj.input.number_of_directions, 2))
- for ex_id, ex_rates in enumerate(firing_rate_array):
- rate_sum = 0.
- for dir_id, dir in enumerate(directions):
- tuning_vectors[ex_id, dir_id] = np.array([np.cos(dir), np.sin(dir)]) * ex_rates[dir_id]
- rate_sum += ex_rates[dir_id]
- if rate_sum != 0.:
- tuning_vectors[ex_id, :, :] /= rate_sum
- tuning_vectors_summed = np.sum(tuning_vectors, axis=1)
- head_direction_indices = np.array([np.linalg.norm(v) for v in tuning_vectors_summed])
- mean_hdi = np.mean(head_direction_indices)
- return 1. - mean_hdi
- if __name__ == "__main__":
- print(np.linspace(0.0, 400.0, 30).tolist()[14])
- env = Environment(trajectory=TRAJ_NAME,
- comment="Compare the head direction tuning for circular and ellipsoid 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("orientation_map")
- traj.f_add_parameter("orientation_map.correlation_length", np.linspace(0.0, 400.0, 30).tolist()[14],
- comment="Correlation length of orientations in um")
- traj.f_add_parameter("orientation_map.seed", 1, comment="Random seed for map generation.")
- traj.f_add_parameter("orientation_map.sheet_size", 450, comment="Sheet size in um")
- traj.f_add_parameter_group("network")
- traj.f_add_parameter("network.N_E", 900, comment="Number of excitatory neurons")
- traj.f_add_parameter("network.N_I", 100, comment="Number of inhibitory neurons")
- traj.f_add_parameter_group("synapse")
- traj.f_add_parameter("synapse.inhibitory", 30, "Strength of conductance-based inhibitory synapse in nS.")
- traj.f_add_parameter("synapse.excitatory", 1, "Strength of conductance-based inhibitory synapse in mV.")
- traj.f_add_parameter_group("input")
- traj.f_add_parameter("input.width", np.pi / 3.0, comment="Standard deviation of incoming head direction input.")
- traj.f_add_parameter("input.baseline", 0.2, comment="Head direction input baseline")
- traj.f_add_parameter("input.amplitude", 0.5, comment="Head direction input amplitude")
- traj.f_add_parameter("input.number_of_directions", 12, comment="Number of probed directions")
- 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.entropy_maximisation.steps", 30, comment="Steps for entropy maximisation")
- 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")
- optimize_interneuron_orientation_by_hdi(traj)
|