123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313 |
- 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()
|