refactored simulation of delayed inhibition into a separate class

scripts/prc_hodkin_huxley_h_current.py

@@ -1,18 +1,9 @@
-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))
+from scripts.timed_inhibition import TimedInhibition
 
 
 def get_mean_period(spike_train):
@@ -212,58 +203,24 @@ 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', 'g_syn']
+record_variables = ['v', 'r', 'g_syn']
 integration_method = 'exponential_euler'
 
 initial_states = {
     "v": hodgkin_huxley_params["El"]
 }
+current_drive = 0.5 * nA
 
-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()
-
+timed_inhibition_experiment = TimedInhibition(neuron_eqs, synapse_model, synapse_on_pre, neuron_params,
+                                              synapse_params, initial_states, spike_threshold, integration_method,
+                                              current_drive, record_variables)
 
-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)
+spike_recorder, neuron_state_recorder = timed_inhibition_experiment.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
 g_state = neuron_state_recorder[1].g_syn
 spike_train = spike_recorder.spike_trains()[0]
 
@@ -284,9 +241,8 @@ inhibitory_delay = np.linspace(1 * ms, period_unperturbed - 1 * ms, n_simulation
 spike_trains_dict = {}
 periods = []
 
-
 for delay in inhibitory_delay:
-    run_sim(delay, inhibition_on, record_states=False, duration=100 * ms)
+    timed_inhibition_experiment.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
@@ -294,16 +250,14 @@ for delay in inhibitory_delay:
 
 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.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])
 
@@ -311,8 +265,6 @@ 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)
@@ -330,7 +282,6 @@ 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)

scripts/timed_inhibition.py

@@ -0,0 +1,58 @@
+import copy
+import warnings
+
+import brian2 as br
+from brian2 import ms
+
+
+class TimedInhibition:
+    def __init__(self, neuron_model, synapse_model, synapse_on_pre, neuron_parameters, synapse_parameters,
+                 initial_states, spike_threshold, integration_method, current_drive, record_variables):
+        threshold_eqs = """
+        v_threshold: volt
+        """
+
+        self.neuron = br.NeuronGroup(N=1, \
+                                     model=neuron_model + threshold_eqs, \
+                                     threshold='v > v_threshold', \
+                                     refractory='v > v_threshold', \
+                                     method=integration_method)
+
+        self.neuron.v_threshold = spike_threshold
+        set_parameters_from_dict(self.neuron, initial_states)
+        self.neuron.I = current_drive
+
+        self.spike_recorder = br.SpikeMonitor(source=self.neuron)
+        self.neuron_state_recorder = br.StateMonitor(self.neuron, record_variables, record=True)
+
+        self.auto_synapse = br.Synapses(source=self.neuron, target=self.neuron, model=synapse_model,
+                                        on_pre=synapse_on_pre,
+                                        delay=0.0 * ms)
+        self.auto_synapse.connect()
+        set_parameters_from_dict(self.auto_synapse, synapse_parameters)
+        self.net = br.Network(self.neuron)
+        self.net.add(self.auto_synapse)
+        self.net.add(self.spike_recorder)
+        self.net.add(self.neuron_state_recorder)
+        self.net.store()
+
+        self.network_params = copy.deepcopy(neuron_parameters)
+        self.network_params.update(synapse_parameters)
+
+    def run_sim(self, delay=0.0, pulse_strength=0.0, record_states=False, duration=50 * ms):
+        self.net.restore()
+        self.neuron_state_recorder.record = record_states
+        self.auto_synapse.delay = delay
+        self.auto_synapse.synaptic_strength = pulse_strength
+
+        self.net.run(duration=duration, namespace=self.network_params)
+
+        return self.spike_recorder, self.neuron_state_recorder
+
+
+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))