import copy import warnings import brian2 as br import matplotlib.pyplot as plt import numpy as np from brian2.units import * from brian2 import implementation from brian2 import network_operation from brian2 import defaultclock from scripts.models import hodgkin_huxley_cond_synapse_weak_noise_eqs, exponential_synapse, \ exponential_synapse_on_pre, exponential_synapse_params, eqs_ih, hodgkin_huxley_params, \ ih_params, delta_synapse_model, delta_synapse_on_pre, noise_params def set_parameters_from_dict(neurongroup, dictionary_of_parameters): for param_key, param_value in dictionary_of_parameters.items(): try: neurongroup.__setattr__(param_key, param_value) except AttributeError as err: warnings.warn("{:s} has no parameter {:s}".format(neurongroup.name, param_key)) ''' Model for the excitatory neurons ''' # Hodgkin Huxley model from Brian2 documentation spike_threshold = 40*mV excitatory_neuron_options = { "threshold" : "v > spike_threshold", "refractory" : "v > spike_threshold", "method" : 'euler' } ''' Model for the interneuron, currently only the inhibition timing is of interest thus the lif model ''' lif_interneuron_eqs = """ dv/dt =1.0/tau* (-v + u_ext) :volt (unless refractory) u_ext = u_ext_const : volt """ lif_interneuron_params = { "tau": 7*ms, "v_threshold": -40*mV, "v_reset": -50*mV, "tau_refractory": 0.0*ms, "u_ext_const" : - 41 * mV } lif_interneuron_options = { "threshold": "v>v_threshold", "reset": "v=v_reset", "refractory": "tau_refractory", 'method' : 'euler' } ''' Synapse Models ''' delta_synapse_delay = 2.0 * ms ''' Setup of the Model ''' ## With voltage delta synapse ei_synapse_options = { 'model' : delta_synapse_model, 'on_pre' : delta_synapse_on_pre, } # # ie_synapse_options = { # 'model' : delta_synapse_model, # 'on_pre' : delta_synapse_on_pre, # 'delay' : delta_synapse_delay # } excitatory_neuron_eqs = hodgkin_huxley_cond_synapse_weak_noise_eqs + eqs_ih excitatory_neuron_params = hodgkin_huxley_params excitatory_neuron_params.update(ih_params) excitatory_neuron_params.update(excitatory_neuron_params) excitatory_neuron_params.update({"E_i": -120 * mV}) excitatory_neuron_params.update(noise_params) # excitatory_neuron_params["stimulus"] = stimulus_array # excitatory_neuron_params["stimulus"] = stimulus_fun # ### With conductance based delta synapse # neuron_eqs = hodgkin_huxley_eqs_with_synaptic_conductance + eqs_ih # # synapse_model = exponential_synapse # synapse_on_pre = exponential_synapse_on_pre # synapse_params = exponential_synapse_params # # neuron_params = hodgkin_huxley_params # neuron_params.update(ih_params) # neuron_params.update({"E_i": -80 * mV}) # # inhibition_off = 0.0 * nS # inhibition_on = 100 * nS network_params = copy.deepcopy(excitatory_neuron_params) network_params.update(lif_interneuron_params) network_params["spike_threshold"] = spike_threshold record_variables = ['v', 'ih', 'I'] integration_method = 'exponential_euler' initial_states = { "v": hodgkin_huxley_params["El"] } threshold_eqs = """ spike_threshold: volt """ ''' Setup of the neuron groups, synapses and ultimately the network ''' excitatory_neurons = br.NeuronGroup(N=3, \ model=excitatory_neuron_eqs, \ threshold=excitatory_neuron_options["threshold"], \ refractory=excitatory_neuron_options["refractory"], \ method=excitatory_neuron_options["method"]) # set_parameters_from_dict(excitatory_neurons, excitatory_neuron_params) Why doesnt this work? set_parameters_from_dict(excitatory_neurons, initial_states) # excitatory_neurons.I_ext_const = 0.15*nA input_current = 0.15*nA excitatory_neurons.I = input_current inhibitory_neurons = br.NeuronGroup(N=1, \ model=lif_interneuron_eqs, \ threshold=lif_interneuron_options["threshold"], \ refractory=lif_interneuron_options["refractory"], \ reset=lif_interneuron_options["reset"], \ method = lif_interneuron_options['method']) # set_parameters_from_dict(inhibitory_neurons, lif_interneuron_params) Same here e_to_i_synapses = br.Synapses(source=excitatory_neurons, \ target=inhibitory_neurons, \ model=ei_synapse_options["model"], \ on_pre=ei_synapse_options["on_pre"]) e_to_i_synapses.connect() i_to_e_synapses = br.Synapses(source=inhibitory_neurons, \ target=excitatory_neurons, \ model=exponential_synapse, \ on_pre=exponential_synapse_on_pre, \ delay=2.0*ms, \ method='euler') i_to_e_synapses.connect() set_parameters_from_dict(i_to_e_synapses,exponential_synapse_params) exc_spike_recorder = br.SpikeMonitor(source=excitatory_neurons) inh_spike_recorder = br.SpikeMonitor(source=inhibitory_neurons) neuron_state_recorder = br.StateMonitor(excitatory_neurons, record_variables, record=[0,1,2], dt=0.1*ms) ''' External Stimulus ''' #First Try # @implementation('cython', ''' # cdef double stimulus_fun(double t): # return 6.25*1e-9 * sin(4. * 3.14 /(1000. 1e-3) * t) # ''') # @check_units(t=second, result=amp) # def stimulus_fun(t): # # return 4 * nA *t /ms # return 6.25* nA * sin(4. * 3.14 /(1000 * ms) * t) #Second Try # pulse = 0.25 # stimulus_array = br.TimedArray([0.,0.,0.,0.,0.,0.,pulse,0.,0.,0.,0.,0.,0.,-pulse,0.,0.,0.,0.,0.,0.] * nA,dt=100.*ms) #Third Try pulse_length = 40*ms pulse_1_start = 200*ms pulse_2_start = 400*ms pulse_3_start = 600*ms pulse_1_and_2_start = 800*ms pulse_strength = 0.6 * nA @network_operation(dt=1*ms) def update_input(t): # Push 1 if t > pulse_1_start and t < pulse_1_start + 2*ms: excitatory_neurons[0].I = pulse_strength elif t > pulse_1_start + pulse_length and t < pulse_1_start + pulse_length + 2*ms: excitatory_neurons[0].I = input_current # Push 2 elif t > pulse_2_start and t < pulse_2_start + 2*ms: excitatory_neurons[1].I = pulse_strength elif t > pulse_2_start + pulse_length and t < pulse_2_start + pulse_length + 2*ms: excitatory_neurons[1].I = input_current # Push 3 elif t > pulse_3_start and t < pulse_3_start + 2 * ms: excitatory_neurons[2].I = pulse_strength elif t > pulse_3_start + pulse_length and t < pulse_3_start + pulse_length + 2 * ms: excitatory_neurons[2].I = input_current # Push 1 and 2 elif t > pulse_1_and_2_start and t < pulse_1_and_2_start + 2 * ms: excitatory_neurons[0].I = pulse_strength excitatory_neurons[1].I = pulse_strength elif t > pulse_1_and_2_start + pulse_length and t < pulse_1_and_2_start + pulse_length + 2 * ms: excitatory_neurons[0].I = input_current excitatory_neurons[1].I = input_current defaultclock.dt = 0.005*ms net = br.Network(update_input) net.add(excitatory_neurons) net.add(inhibitory_neurons) net.add(e_to_i_synapses) net.add(i_to_e_synapses) net.add(exc_spike_recorder) net.add(inh_spike_recorder) net.add(neuron_state_recorder) net.store() ''' Run of the simulation ''' def plot_spiking(): ex_spike_trains = exc_spike_recorder.spike_trains() in_spike_trains = inh_spike_recorder.spike_trains() fig = plt.figure() ax = fig.add_subplot(111) for key, times in ex_spike_trains.items(): ax.plot(times / ms, key / 2.0 * np.ones(times.shape), 'b|') offset = 2 for key, times in in_spike_trains.items(): ax.plot(times / ms, (key + offset) / 2.0 * np.ones(times.shape), 'r|') ax.grid(axis='x') ax.set_ylim(-0.1, 1.1) ax.set_xlabel("Time(ms)"); def run_sim(inh_strength, exc_strength, record_states=True, run_time=50 * ms): net.restore() e_to_i_synapses.synaptic_strength = exc_strength i_to_e_synapses.synaptic_strength = inh_strength # excitatory_neurons.stimulus = get_sinoid_stimulus(0.25 * nA, 4 * np.pi / run_time) # network_params["omega"] = 4 * np.pi / run_time # excitatory_neurons.I_ext_diff = [0.25, -0.25] * nA # neuron_state_recorder.record = record_states net.run(duration=run_time, namespace=network_params) run_time = 1000. * ms inh_strength = 150 * nS exc_strength = lif_interneuron_params["v_threshold"]-lif_interneuron_params["v_reset"] + 5*mV run_sim(inh_strength=inh_strength, exc_strength=exc_strength, record_states=True, run_time=run_time) plot_spiking() print(neuron_state_recorder.v.shape) t_state = neuron_state_recorder.t v1_state = neuron_state_recorder[0].v v2_state = neuron_state_recorder[1].v v3_state = neuron_state_recorder[2].v I1_state = neuron_state_recorder[0].I I2_state = neuron_state_recorder[1].I I3_state = neuron_state_recorder[2].I fig, axs = plt.subplots(4, sharex=True) axs[0].plot(t_state / ms, v1_state / mV, 'C1') axs[1].plot(t_state / ms, v2_state / mV, 'C2') axs[2].plot(t_state / ms, v3_state / mV, 'C3') axs[3].plot(t_state / ms, I1_state / nA, 'C1') axs[3].plot(t_state / ms, I2_state / nA + 0.01, 'C2') axs[3].plot(t_state / ms, I3_state / nA + 0.02, 'C3') run_sim(inh_strength=0*nS, exc_strength=10.*mV, record_states=True, run_time=run_time) plot_spiking() plt.show() plt.show()