From 037ebeb9b4a8d9a1eef9080ac16251df950a49f4 Mon Sep 17 00:00:00 2001 From: Jan Grewe Date: Mon, 7 Sep 2020 13:09:59 +0200 Subject: [PATCH] add model and p-unit respones --- model.py | 83 +++++++++++++++++++++++++++++++ punit_responses.py | 119 +++++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 202 insertions(+) create mode 100644 model.py create mode 100644 punit_responses.py diff --git a/model.py b/model.py new file mode 100644 index 0000000..b1155ef --- /dev/null +++ b/model.py @@ -0,0 +1,83 @@ +from IPython.terminal.embed import embed +from numba.core.types import float64 +import numpy as np + +try: + from numba import jit +except ImportError: + def jit(nopython): + def decorator_jit(func): + return func + return decorator_jit + + +def load_models(file): + """ Load model parameter from csv file. + + Parameters + ---------- + file: string + Name of file with model parameters. + + Returns + ------- + parameters: list of dict + For each cell a dictionary with model parameters. + """ + parameters = [] + with open(file, 'r') as file: + header_line = file.readline() + header_parts = header_line.strip().split(",") + keys = header_parts + for line in file: + line_parts = line.strip().split(",") + parameter = {} + for i in range(len(keys)): + parameter[keys[i]] = float(line_parts[i]) if i > 0 else line_parts[i] + parameters.append(parameter) + return parameters + + +@jit(nopython=True) +def simulate(stimulus, deltat=0.00005, v_zero=0.0, a_zero=2.0, threshold=1.0, v_base=0.0, + delta_a=0.08, tau_a=0.1, v_offset=-10.0, mem_tau=0.015, noise_strength=0.05, + input_scaling=60.0, dend_tau=0.001, ref_period=0.001, **kwargs): + """ Simulate a P-unit. + + Returns + ------- + spike_times: 1-D array + Simulated spike times in seconds. + """ + # initial conditions: + v_dend = stimulus[0] + v_mem = v_zero + adapt = a_zero + + # prepare noise: + noise = np.random.randn(len(stimulus)) + noise *= noise_strength / np.sqrt(deltat) + + # rectify stimulus array: + stimulus = stimulus.copy() + stimulus[stimulus < 0.0] = 0.0 + + # integrate: + spike_times = [] + for i in range(len(stimulus)): + v_dend += (-v_dend + stimulus[i]) / dend_tau * deltat + v_mem += (v_base - v_mem + v_offset + ( + v_dend * input_scaling) - adapt + noise[i]) / mem_tau * deltat + adapt += -adapt / tau_a * deltat + + # refractory period: + if len(spike_times) > 0 and (deltat * i) - spike_times[-1] < ref_period + deltat/2: + v_mem = v_base + + # threshold crossing: + if v_mem > threshold: + v_mem = v_base + spike_times.append(i * deltat) + adapt += delta_a / tau_a + + return np.array(spike_times) diff --git a/punit_responses.py b/punit_responses.py new file mode 100644 index 0000000..cb4b311 --- /dev/null +++ b/punit_responses.py @@ -0,0 +1,119 @@ +import numpy as np +import nixio as nix +import os +from numpy.core.fromnumeric import repeat +from traitlets.traitlets import Instance +from chirp_ams import get_signals +from model import simulate, load_models +from IPython import embed +import matplotlib.pyplot as plt + +def append_settings(section, sec_name, sec_type, settings): + section = section.create_section(sec_name, sec_type) + for k in settings.keys(): + if isinstance(settings[k], dict): + append_settings(section, k, "settings", settings[k]) + else: + if isinstance(settings[k], np.ndarray): + if len(settings[k].shape) == 1: + section[k] = list(settings[k]) + else: + section[k] = settings[k] + + +def save(filename, name, stimulus_settings, model_settings, stimulus, responses, overwrite=False): + print("saving! ", filename, name) + if os.path.exists(filename) and not overwrite: + nf = nix.File.open(filename, nix.FileMode.ReadWrite) + else: + nf = nix.File.open(filename, mode=nix.FileMode.Overwrite, + compression=nix.Compression.DeflateNormal) + + if name in nf.blocks: + print("Data with this name is already stored! ", name) + nf.close() + return + mdata = nf.create_section(name, "nix.simulation") + append_settings(mdata, "model parameter", "nix.model.settings", model_settings) + append_settings(mdata, "stimulus parameter", "nix.stimulus.settings", stimulus_settings) + b = nf.create_block(name, "nix.simulation") + b.metadata = mdata + + # save stimulus + stim_da = b.create_data_array("stimulus", "nix.timeseries.sampled", dtype=nix.DataType.Float, + data=stimulus) + stim_da.label = "voltage" + stim_da.label = "mv/cm" + dim = stim_da.append_sampled_dimension(model_settings["deltat"]) + dim.label = "time" + dim.unit = "s" + + # save responses + for i in range(len(responses)): + da = b.create_data_array("response_%i" %i, "nix.timeseries.events.spike_times", + dtype=nix.DataType.Float, data=responses[i]) + da.label = "time" + da.unit = "s" + dim = da.append_range_dimension() + dim.link_data_array(da, [-1]) + nf.close() + pass + + +def plot_responses(): + pass + + +def simulate_responses(stimulus_params, model_params, repeats=10): + cell_params = model_params.copy() + cell = cell_params["cell"] + del cell_params["cell"] + del cell_params["EODf"] + for c in stimulus_params["contrasts"]: + print("creating stimuli\n\tcontrast: ", str(c), "\t condition: ", + stimulus_params["condition"]) + params = stimulus_params.copy() + del params["contrasts"] + del params["chirp_frequency"] + params["contrast"] = c + time, self_signal, self_freq, other_signal, other_freq = get_signals(**params) + signal = (self_signal + other_signal) + signal /= np.max(signal) + print("create p-unit responses for cell: ", cell) + + spikes = [] + for r in range(repeats): + spikes.append(simulate(signal, **cell_params)) + save("test.nix", "contrast_%.3f_condition_%s" %(c, stimulus_params["condition"]), params, + cell_params, signal, spikes) + + +def main(): + cell_id = 20 + models = load_models("models.csv") + deltaf = 20. # Hz, difference frequency between self and other + model_params = models[cell_id] + stimulus_params = {"eodfs": {"self": model_params["EODf"], + "other": model_params["EODf"] + deltaf}, + "contrasts": [20, 10, 5, 2.5, 1.25, 0.625, 0.3125], + "chirp_size": 100, # Hz, frequency excursion + "chirp_duration": 0.015, # s, chirp duration + "chirp_amplitude_dip": 0.05, # %, amplitude drop during chirp + "chirp_frequency": 5, # Hz, how often does the fish chirp + "duration": 5., # s, total duration of simulation + "dt": model_params["deltat"], # s, stepsize of the simulation + } + chirp_times = np.arange(stimulus_params["chirp_duration"], + stimulus_params["duration"] - stimulus_params["chirp_duration"], + 1./stimulus_params["chirp_frequency"]) + stimulus_params["chirp_times"] = chirp_times + conditions = ["other", "self"] + + for c in conditions: + stimulus_params["condition"] = c + simulate_responses(stimulus_params, model_params, repeats=25) + pass + + +if __name__ == "__main__": + main() \ No newline at end of file