structured_inhibition_spatial_network_model/scripts/competition_gain/parameter_scan_merged_micro_circuit_empirical.py

import copy
import warnings

import brian2 as br
import matplotlib.pyplot as plt
import numpy as np
from brian2.units import *
from brian2 import network_operation

from scripts.models import hodgkin_huxley_eqs_with_synaptic_conductance, exponential_synapse, \
    exponential_synapse_on_pre, exponential_synapse_params, eqs_ih, hodgkin_huxley_params, \
    ih_params, delta_synapse_model, delta_synapse_on_pre


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" : 'exponential_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",
}



'''
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_eqs_with_synaptic_conductance + 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["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', '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=2, \
                                    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"])
# 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)
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])

'''
External Stimulus
'''

pulse_length = 50*ms
pulse_1_start = 200*ms
pulse_2_start = 400*ms
pulse_3_start = 600*ms
pulse_strength = 0.6 * nA
@network_operation(dt=1*ms)
def update_input(t):
    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
    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
    elif t > pulse_3_start and t < pulse_3_start + 2 * ms:
        excitatory_neurons[0].I = pulse_strength
    elif t > pulse_3_start + pulse_length and t < pulse_3_start + pulse_length + 2 * ms:
        excitatory_neurons[0].I = input_current




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

    # neuron_state_recorder.record = record_states

    net.run(duration=run_time, namespace=network_params)

run_time = 1000. * ms
exc_strength = lif_interneuron_params["v_threshold"]-lif_interneuron_params["v_reset"] + 5*mV


'''
Gain in merged microcircuit
'''

input_current = 0.15*nA

n_ghbar = 5
ghbar_range = linspace(0.0,50.0,n_ghbar)*nS
n_syn_strength = 5
synaptic_strength_range = linspace(0.0,160.0,n_syn_strength)*nS

bistable_array = ndarray((n_ghbar,n_syn_strength), dtype=float)

# Incredibly, enumerate doesn't work with Quantity objects
syn_str_id = 0
for synaptic_strength_val in synaptic_strength_range:
    ghbar_id = 0
    for ghbar_val in ghbar_range:
        net.restore()
        network_params['ghbar'] = ghbar_val
        net.store()

        run_sim(inh_strength=synaptic_strength_val, exc_strength=exc_strength, record_states=True, run_time=run_time)

        period_unperturbed, inhibitory_delay, periods = timed_inhibition_experiment.measure_response_curve(n_simulations,
                                                                                                           inhibition_on,
                                                                                                           offset_before_spike_peak=before_spike,
                                                                                                           offset_after_spike_peak=after_spike)

        isbistable = bistability_from_iterative_map_unidirectional(period_unperturbed, inhibitory_delay, periods)
        print(synaptic_strength_val,ghbar_val,isbistable)

        bistable_array[ghbar_id,syn_str_id] = isbistable
        ghbar_id = ghbar_id + 1
    syn_str_id = syn_str_id + 1



rates_disconnected = np.array(rates_disconnected)
rates_merged = np.array(rates_merged)
rates_unidirectional = np.array(rates_unidirectional)

fig, axs = plt.subplots(2, sharex=True, sharey=True)
axs[0].plot(ex_drive_range / nA, rate_diff_disconnected * ms, 'C1',label='disconnected')
axs[0].plot(ex_drive_range / nA, rate_diff_merged * ms, 'C2', label='merged')
axs[0].plot(ex_drive_range / nA, rate_diff_unidirectional * ms, 'C3', label='unidirectional')

axs[1].plot(ex_drive_range / nA,rates_disconnected[:,0] * ms,'C1')
axs[1].plot(ex_drive_range / nA,rates_disconnected[:,1] * ms,'C1',linestyle='--')

axs[1].plot(ex_drive_range / nA,rates_merged[:,0] * ms,'C2')
axs[1].plot(ex_drive_range / nA,rates_merged[:,1] * ms,'C2',linestyle='--')

axs[1].plot(ex_drive_range / nA,rates_unidirectional[:,0] * ms,'C3')
axs[1].plot(ex_drive_range / nA,rates_unidirectional[:,1] * ms,'C3',linestyle='--')

axs[1].set_xlabel('Input Current [nA]')
axs[0].set_ylabel('Firing Rate [Hz]')
axs[1].set_ylabel('Firing Rate [Hz]')
axs[0].legend()
# run_sim(inh_strength=0*nS, exc_strength=10.*mV, record_states=True, run_time=run_time)
# plot_spiking()
#
# plt.show()





plt.show()