from stimuli.AbstractStimulus import AbstractStimulus from models.AbstractModel import AbstractModel import numpy as np import functions as fu from numba import jit import time class LifacNoiseModel(AbstractModel): # all times in milliseconds # possible mem_res: 100 * 1000000 exact value unknown in p-units DEFAULT_VALUES = {"mem_tau": 20, "v_base": 0, "v_zero": 0, "threshold": 1, "v_offset": 50, "input_scaling": 1, "delta_a": 0.4, "tau_a": 40, "a_zero": 0, "noise_strength": 3, "step_size": 0.01} def __init__(self, params: dict = None): super().__init__(params) self.voltage_trace = [] self.adaption_trace = [] self.spiketimes = [] self.stimulus = None # self.frequency_trace = [] def simulate(self, stimulus: AbstractStimulus, total_time_s): self.stimulus = stimulus output_voltage = [] adaption = [] spiketimes = [] current_v = self.parameters["v_zero"] current_a = self.parameters["a_zero"] for time_point in np.arange(0, total_time_s*1000, self.parameters["step_size"]): # rectified input: stimulus_strength = fu.rectify(stimulus.value_at_time_in_ms(time_point)) v_next = self._calculate_voltage_step(current_v, stimulus_strength - current_a) a_next = self._calculate_adaption_step(current_a) if v_next > self.parameters["threshold"]: v_next = self.parameters["v_base"] spiketimes.append(time_point/1000) a_next += self.parameters["delta_a"] / (self.parameters["tau_a"] / 1000) output_voltage.append(v_next) adaption.append(a_next) current_v = v_next current_a = a_next self.voltage_trace = output_voltage self.adaption_trace = adaption self.spiketimes = spiketimes return output_voltage, spiketimes def simulate_fast(self, stimulus: AbstractStimulus, total_time_s): v_zero = self.parameters["v_zero"] a_zero = self.parameters["a_zero"] step_size = self.parameters["step_size"] threshold = self.parameters["threshold"] v_base = self.parameters["v_base"] delta_a = self.parameters["delta_a"] tau_a = self.parameters["tau_a"] v_offset = self.parameters["v_offset"] mem_tau = self.parameters["mem_tau"] noise_strength = self.parameters["noise_strength"] stimulus_array = stimulus.as_array(total_time_s, step_size) parameters = np.array([v_zero, a_zero, step_size, threshold, v_base, delta_a, tau_a, v_offset, mem_tau, noise_strength]) voltage_trace, adaption, spiketimes = simulate_fast(stimulus_array, total_time_s, parameters) self.stimulus = stimulus self.voltage_trace = voltage_trace self.adaption_trace = adaption self.spiketimes = spiketimes return voltage_trace, spiketimes def _calculate_voltage_step(self, current_v, input_v): v_base = self.parameters["v_base"] step_size = self.parameters["step_size"] v_offset = self.parameters["v_offset"] mem_tau = self.parameters["mem_tau"] noise_strength = self.parameters["noise_strength"] noise_value = np.random.normal() noise = noise_strength * noise_value / np.sqrt(step_size) return current_v + step_size * ((v_base - current_v + v_offset + input_v + noise) / mem_tau) def _calculate_adaption_step(self, current_a): step_size = self.parameters["step_size"] return current_a + (step_size * (-current_a)) / self.parameters["tau_a"] def min_stimulus_strength_to_spike(self): return self.parameters["threshold"] - self.parameters["v_base"] def get_sampling_interval(self): return self.parameters["step_size"] def get_frequency(self): # TODO also change simulates_frequency() if any calculation is added! raise NotImplementedError("No calculation implemented yet for the frequency.") def get_spiketimes(self): return self.spiketimes def get_voltage_trace(self): return self.voltage_trace def get_adaption_trace(self): return self.adaption_trace def simulates_frequency(self) -> bool: return False def simulates_spiketimes(self) -> bool: return True def simulates_voltage_trace(self) -> bool: return True def get_recording_times(self): # [delay, stimulus_start, stimulus_duration, time_to_end] self.stimulus = AbstractStimulus() delay = 0 start = self.stimulus.get_stimulus_start_s() duration = self.stimulus.get_stimulus_duration_s() total_time = len(self.voltage_trace) / self.parameters["step_size"] return [delay, start, duration, total_time] def get_model_copy(self): return LifacNoiseModel(self.parameters) def calculate_baseline_markers(self, stimulus_freq=750): """ calculates the baseline markers baseline frequency, vector strength and serial correlation based on simulated 30 seconds with a standard Sinusoidal stimulus with the given frequency :return: baseline_freq, vs, sc """ pass def calculate_fi_markers(self, contrasts, ): """ calculates the fi markers f_infinity, f_infinity_slope for given contrasts based on simulated 2 seconds for each contrast :return: """ def stimulus_to_numpy_array(stimulus: AbstractStimulus, total_time_s, step_size): total_time_points = int(total_time_s * 1000 / step_size) stimulus_values = np.zeros(total_time_points) for idx in range(len(stimulus_values)): # rectified input: stimulus_values[idx] = fu.rectify(stimulus.value_at_time_in_ms(step_size*idx)) return stimulus_values @jit(nopython=True) def simulate_fast(rectified_stimulus_array, total_time_s, parameters: np.ndarray): v_zero = parameters[0] a_zero = parameters[1] step_size = parameters[2] threshold = parameters[3] v_base = parameters[4] delta_a = parameters[5] tau_a = parameters[6] v_offset = parameters[7] mem_tau = parameters[8] noise_strength = parameters[9] time = np.arange(0, total_time_s * 1000, step_size) length = len(time) output_voltage = np.zeros(length) adaption = np.zeros(length) stimulus_values = rectified_stimulus_array spiketimes = [] output_voltage[0] = v_zero adaption[0] = a_zero for i in range(len(time)-1): noise_value = np.random.normal() noise = noise_strength * noise_value / np.sqrt(step_size) output_voltage[i] = output_voltage[i-1] + ((v_base - output_voltage[i-1] + v_offset + stimulus_values[i] - adaption[i-1] + noise) / mem_tau) * step_size adaption[i] = adaption[i-1] + ((-adaption[i-1]) / tau_a) * step_size if output_voltage[i] > threshold: output_voltage[i] = v_base spiketimes.append(i*step_size) adaption[i] += delta_a / (tau_a / 1000) return output_voltage, adaption, spiketimes