123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347 |
- import copy
- import warnings
- import brian2 as br
- import matplotlib.pyplot as plt
- import numpy as np
- from brian2.units import *
- 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))
- def get_mean_period(spike_train):
- return (np.max(spike_train) - np.min(spike_train)) / (spike_train.shape[0] - 1)
- def get_difference_in_periods(spike_trains_dict, t_isi):
- mean_response_dict = {}
- # t_isi = get_mean_period(spike_trains_dict[0.0 * ms])
- for delay, spike_train in spike_trains_dict.items():
- mean_response_dict[delay] = get_mean_period(spike_train) - t_isi
- return mean_response_dict
- def plot_spiking(spike_trains_dict):
- fig = plt.figure()
- ax = fig.add_subplot(111)
- for key, times in spike_trains_dict.items():
- ax.plot(times / ms, key * np.ones(times.shape), 'b.')
- # ax.plot(spike_train[0]/ms, np.ones(spike_train[0].shape), 'b|')
- ax.grid(axis='x')
- # ax.set_ylim(-0.1, 1.1)
- ax.set_xlabel("Time(ms)")
- def plot_mean_period(mean_period_dict):
- fig = plt.figure()
- ax = fig.add_subplot(111)
- ax.plot(list(mean_period_dict.keys()), list(mean_period_dict.values()), 'b')
- ax.grid(axis='x')
- ax.set_xlabel("Delay (ms)")
- def plot_mean_response(mean_response_dict):
- fig = plt.figure()
- ax = fig.add_subplot(111)
- ax.plot(np.array(list(mean_response_dict.keys())), 1e3 * np.array(list(mean_response_dict.values())), 'b')
- ax.plot(list(mean_response_dict.keys()), list(mean_response_dict.keys()), 'k--')
- ax.grid(axis='x')
- ax.set_xlabel("t_inh (ms)")
- ax.set_ylabel("exc_spike_delay (ms)")
- def plot_delay_minus_response(mean_response_dict):
- fig = plt.figure()
- ax = fig.add_subplot(111)
- print(np.array(list(mean_response_dict.keys())))
- print(np.array(list(mean_response_dict.values())))
- print(np.array(list(mean_response_dict.keys())) - np.array(list(mean_response_dict.values())))
- ax.plot(list(mean_response_dict.keys()),
- np.array(list(mean_response_dict.keys())) - 1e3 * np.array(list(mean_response_dict.values())), 'b')
- ax.grid(axis='x')
- ax.set_xlabel("Delay (ms)");
- def plot_phase_response_curve(mean_response_dict, t_isi):
- phi_inh = np.array(list(mean_response_dict.keys())) * ms / t_isi
- delta_phi = -np.array(list(mean_response_dict.values())) * second / t_isi
- fig = plt.figure()
- ax = fig.add_subplot(111)
- ax.plot(phi_inh, delta_phi, 'b')
- ax.plot(phi_inh, -phi_inh, 'k--')
- ax.grid(axis='x')
- ax.set_xlabel("phase of inhibition (ms)")
- ax.set_ylabel("phase shift (ms)");
- # Hodgkin Huxley model from Brian2 documentation
- area = 20000 * umetre ** 2
- hodgkin_huxley_params = {
- "Cm": 1 * ufarad * cm ** -2 * area,
- "gl": 5e-5 * siemens * cm ** -2 * area,
- "El": -65 * mV,
- "EK": -90 * mV,
- "ENa": 50 * mV,
- "g_na": 100 * msiemens * cm ** -2 * area,
- "g_kd": 30 * msiemens * cm ** -2 * area,
- "VT": -63 * mV
- }
- # The model
- hodgkin_huxley_eqs = br.Equations('''
- dv/dt = (gl*(El-v) - g_na*(m*m*m)*h*(v-ENa) - g_kd*(n*n*n*n)*(v-EK) + I + ih)/Cm : volt
- dm/dt = 0.32*(mV**-1)*(13.*mV-v+VT)/
- (exp((13.*mV-v+VT)/(4.*mV))-1.)/ms*(1-m)-0.28*(mV**-1)*(v-VT-40.*mV)/
- (exp((v-VT-40.*mV)/(5.*mV))-1.)/ms*m : 1
- dn/dt = 0.032*(mV**-1)*(15.*mV-v+VT)/
- (exp((15.*mV-v+VT)/(5.*mV))-1.)/ms*(1.-n)-.5*exp((10.*mV-v+VT)/(40.*mV))/ms*n : 1
- dh/dt = 0.128*exp((17.*mV-v+VT)/(18.*mV))/ms*(1.-h)-4./(1+exp((40.*mV-v+VT)/(5.*mV)))/ms*h : 1
- I : amp
- ''')
- hodgkin_huxley_eqs_with_synaptic_conductance = br.Equations('''
- dv/dt = (gl*(El-v) - g_na*(m*m*m)*h*(v-ENa) - g_kd*(n*n*n*n)*(v-EK) + I + ih - g_syn*(v-E_i))/Cm : volt
- dm/dt = 0.32*(mV**-1)*(13.*mV-v+VT)/
- (exp((13.*mV-v+VT)/(4.*mV))-1.)/ms*(1-m)-0.28*(mV**-1)*(v-VT-40.*mV)/
- (exp((v-VT-40.*mV)/(5.*mV))-1.)/ms*m : 1
- dn/dt = 0.032*(mV**-1)*(15.*mV-v+VT)/
- (exp((15.*mV-v+VT)/(5.*mV))-1.)/ms*(1.-n)-.5*exp((10.*mV-v+VT)/(40.*mV))/ms*n : 1
- dh/dt = 0.128*exp((17.*mV-v+VT)/(18.*mV))/ms*(1.-h)-4./(1+exp((40.*mV-v+VT)/(5.*mV)))/ms*h : 1
- I : amp
- g_syn: siemens
- ''')
- # # First H-current from Izhikevich p. 48
- # ghbar = 40. * nS #Other values would be 0.5, 2, 3.5, 20 depending on neuron type (Rothman, Manis)
- # Eh = -43*mV
- # V_half_h = -75.0
- # k_h = -5.5
- # V_max_h = -75.
- # sigma_h = 15.
- # C_amp_h = 1000.
- # C_base_h = 100.
- # eqs_ih = """
- # ih = ghbar*r*(Eh-v) : amp
- # dr/dt= (rinf-r)/rtau : 1
- # rinf = 1. / (1+exp((V_half_h - v/mV) / k_h)) : 1
- # rtau = (C_base_h + C_amp_h * exp(-(V_max_h-v/mV)**2./sigma_h**2)) * ms : second
- # """
- # Second H-current from Izhikevich p. 48
- ih_params = {
- "ghbar": 60. * nS, # Other values would be 0.5, 2, 3.5, 20 depending on neuron type (Rothman, Manis)
- "Eh": -1. * mV,
- "V_half_h": -98.0,
- "k_h": -5.5,
- "V_max_h": -75.,
- "sigma_h": 25.,
- "C_amp_h": 60.,
- "C_base_h": 40.,
- }
- eqs_ih = """
- ih = ghbar*r*(Eh-v) : amp
- dr/dt= (rinf-r)/rtau : 1
- rinf = 1. / (1+exp((V_half_h - v/mV) / k_h)) : 1
- rtau = 0.1*(C_base_h + C_amp_h * exp(-(V_max_h-v/mV)**2./sigma_h**2)) * ms : second
- """
- # eqs_ih = """
- # ih = ghbar*r*(Eh-v) : amp
- # dr/dt=(rinf-r)/rtau : 1
- # rinf = 1. / (1+exp((v/mV + 90.) / 7.)) : 1
- # # rtau = ((100000. / (237.*exp((v/mV+60.) / 12.) + 17.*exp(-(v/mV+60.) / 14.))) + 25.)*ms : second
- # rtau = 0.01*(100. + 1000. * exp(-(-76.-v/mV)**2./15.**2)) * ms : second #From Izhikevich
- # """
- delta_synapse_model = 'synaptic_strength: volt'
- delta_synapse_on_pre = 'v+=synaptic_strength'
- delta_synapse_param = {}
- exponential_synapse = """
- dg/dt = -g/tau_syn : siemens
- g_syn_post = g :siemens (summed)
- tau_syn : second
- synaptic_strength : siemens
- """
- exponential_synapse_on_pre = "g+=synaptic_strength"
- exponential_synapse_params = {
- "tau_syn": 5 * ms
- }
- spike_threshold = +40 * mV
- ## With voltage delta synapse
- neuron_eqs = hodgkin_huxley_eqs+eqs_ih
- synapse_model = delta_synapse_model
- synapse_on_pre = delta_synapse_on_pre
- synapse_params = delta_synapse_param
- neuron_params = hodgkin_huxley_params
- neuron_params.update(ih_params)
- inhibition_off = 0.0 * mV
- inhibition_on = -30*mV
- # ### 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(neuron_params)
- network_params.update(synapse_params)
- record_variables = ['v', 'ih']
- integration_method = 'exponential_euler'
- initial_states = {
- "v": hodgkin_huxley_params["El"]
- }
- threshold_eqs = """
- v_threshold: volt
- """
- neuron = br.NeuronGroup(N=1, \
- model=neuron_eqs + threshold_eqs, \
- threshold='v > v_threshold', \
- refractory='v > v_threshold', \
- method=integration_method)
- neuron.v_threshold = spike_threshold
- set_parameters_from_dict(neuron, initial_states)
- neuron.I = 0.5 * nA
- spike_recorder = br.SpikeMonitor(source=neuron)
- neuron_state_recorder = br.StateMonitor(neuron, record_variables, record=True)
- auto_synapse = br.Synapses(source=neuron, target=neuron, model=synapse_model, on_pre=synapse_on_pre,
- delay=0.0 * ms)
- auto_synapse.connect()
- set_parameters_from_dict(auto_synapse, synapse_params)
- net = br.Network(neuron)
- net.add(auto_synapse)
- net.add(spike_recorder)
- net.add(neuron_state_recorder)
- net.store()
- def run_sim(delay=0.0, pulse_strength=0.0, record_states=False, duration=50 * ms):
- net.restore()
- neuron_state_recorder.record = record_states
- auto_synapse.delay = delay
- auto_synapse.synaptic_strength = pulse_strength
- net.run(duration=duration, namespace=network_params)
- run_sim(delay=3 * ms, record_states=True, pulse_strength=inhibition_off, duration=200 * ms)
- t_state = neuron_state_recorder.t
- v_state = neuron_state_recorder[1].v
- ih_state = neuron_state_recorder[1].ih
- spike_train = spike_recorder.spike_trains()[0]
- fig, ax1 = plt.subplots()
- ax1.plot(t_state / ms, v_state / ms, 'C2')
- plt.show()
- period_unperturbed = get_mean_period(spike_train)
- print(period_unperturbed)
- n_simulations = 40
- inhibitory_delay = np.linspace(1 * ms, period_unperturbed - 1 * ms, n_simulations)
- spike_trains_dict = {}
- periods = []
- for delay in inhibitory_delay:
- run_sim(delay, inhibition_on, record_states=False, duration=100 * ms)
- spike_train = spike_recorder.spike_trains()[0]
- # print(spike_train)
- spike_trains_dict[delay / ms] = spike_train # Careful, delay is not a Quantity anymore because in must be hashable
- periods.append(get_mean_period(spike_train))
- periods = np.array(periods)
- plt.plot(inhibitory_delay/ms, periods/ms)
- plt.hlines(period_unperturbed/ms, inhibitory_delay[0]/ms, inhibitory_delay[-1]/ms, label="No inhibition")
- plt.xlabel("Delay spike to inhibition (ms)")
- plt.ylabel("Period (ms")
- plt.legend()
- t_diff_dict = get_difference_in_periods(spike_trains_dict, period_unperturbed)
- phases = np.array([delay / period_unperturbed for delay in inhibitory_delay])
- phase_diff = np.array([-t_diff_dict[delay / ms] / period_unperturbed for delay in inhibitory_delay])
- optimal_inhibitory_delay_idx = np.argmax(phase_diff)
- optimal_inhibitory_delay = list(t_diff_dict.keys())[optimal_inhibitory_delay_idx] # assumes a single maximum
- optimal_inhibitory_delay_phase = phases[optimal_inhibitory_delay_idx] # assumes a single maximum
- plt.figure()
- plt.plot(phases, phase_diff)
- plt.vlines(optimal_inhibitory_delay, -1, 1)
- plt.title("Phase response curve")
- plt.ylim(-1, 1)
- plt.figure()
- if optimal_inhibitory_delay_idx == 0:
- delta_i = phases
- else:
- delta_i = phases[:-optimal_inhibitory_delay_idx]
- print(delta_i.shape)
- delta_iplus1 = delta_i + 1 + phase_diff[optimal_inhibitory_delay_idx:] - phase_diff[optimal_inhibitory_delay_idx]
- plt.plot(delta_i, delta_iplus1)
- plt.plot(delta_i, delta_i, '--')
- plt.hlines(1, 0, 1)
- plt.xlabel("Delta i")
- plt.ylabel("Delta i+1")
- plt.title("Iterative map")
- plt.ylim(0, 1.2)
- plt.xlim(0, 1)
- # plot_mean_response(t_diff_dict)
- # plot_phase_response_curve(t_diff_dict, t_isi)
- # plot_spiking(spike_trains_dict)
- plt.show()
|