Alexander/P-unit_model
2
0
Fork 0
Code Issues Pull Requests Releases Wiki Activity

version 1 results

Browse Source
This commit is contained in:
a.ott 2020-09-20 13:35:28 +02:00
parent 64d696935e
commit 44174e608f
40 changed files with 351 additions and 237 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

12
Figure_constants.py
View File

@ -19,13 +19,13 @@ FIG_SIZE_LARGE_HIGH = (6, 8)
colors_muted = ptc.colors_muted
COLOR_DATA = ptc.colors_vivid["blue"]
COLOR_DATA_f0 = colors_muted["red"]
COLOR_DATA_finf = ptc.colors_vivid["blue"]
COLOR_DATA = "C0" # '#0020C075' # ptc.colors_vivid["blue"] with alpha=50
COLOR_DATA_f0 = '#C02717' # colors_muted["red"]
COLOR_DATA_finf = '#0020C0' # ptc.colors_vivid["blue"]
COLOR_MODEL = colors_muted["orange"]
COLOR_MODEL_f0 = colors_muted["orange"]
COLOR_MODEL_finf = ptc.colors_vivid["green"]
COLOR_MODEL = "C1" # "#F7801775" # colors_muted["orange"] with alpha=75
COLOR_MODEL_f0 = "C1" # "#F7801775" # colors_muted["orange"]
COLOR_MODEL_finf = "#30D700" # ptc.colors_vivid["green"]
f0_marker = (3, 0, 0) # "^"
finf_marker = (4, 0, 0) # "s"

24
Figures_Baseline.py
View File

@ -63,10 +63,10 @@ def p_unit_heterogeneity():
event_tick_length = (y_lims[1] - y_lims[0]) / 10
ax.eventplot([s * 1000 for s in spikes if time_offset <= s < time_offset+duration],
colors="black", lineoffsets=max(v1[idx_start:idx_end])+1.5, linelengths=event_tick_length)
ax.set_ylabel("Voltage in mV")
ax.set_ylabel("Voltage [mV]")
ax.set_xlim((0, duration*1000))
if i == 2:
ax.set_xlabel("Time in ms")
ax.set_xlabel("Time [ms]")
ax.set_yticks([-5, 5, 15])
for i, cell in enumerate(cells):
@ -81,7 +81,7 @@ def p_unit_heterogeneity():
ax.set_ylabel("Density")
ax.set_yticklabels(["{:.1f}".format(t) for t in ax.get_yticks()])
if i == 2:
ax.set_xlabel("ISI in EOD periods")
ax.set_xlabel("ISI [EOD periods]")
plt.tight_layout()
fig.align_ylabels()
@ -119,8 +119,8 @@ def p_unit_example():
ax.plot((np.array(time[:int(duration/step)]) - start) * 1000, v1[:int(duration/step)], consts.COLOR_DATA)
ax.eventplot([s * 1000 for s in spiketimes if start < s < start + duration], lineoffsets=max(v1[:int(duration/step)])+1.25,
color="black", linelengths=2)
ax.set_ylabel('Voltage in mV')
ax.set_xlabel('Time in ms')
ax.set_ylabel('Voltage [mV]')
ax.set_xlabel('Time [ms]')
ax.set_title("Baseline Firing")
ax.set_xlim((0, duration*1000))
@ -134,7 +134,7 @@ def p_unit_example():
bins = np.arange(0, maximum * 1.01, 0.1)
ax.hist(isi, bins=bins, color=consts.COLOR_DATA, density=True)
ax.set_ylabel("Density")
ax.set_xlabel("ISI in EOD periods")
ax.set_xlabel("ISI [EOD periods]")
ax.set_title("ISI Histogram")
@ -170,8 +170,8 @@ def p_unit_example():
ax.set_xlim((-0.2, part-0.2))
ylim = ax.get_ylim()
ax.set_ylim((0, ylim[1]))
ax.set_xlabel("Time in s")
ax.set_ylabel("Frequency in Hz")
ax.set_xlabel("Time [s]")
ax.set_ylabel("Frequency [Hz]")
ax.set_title("Step Response")
# FI-Curve
@ -193,7 +193,7 @@ def p_unit_example():
# ax.set_xlim((0, 10))
# ax.set_ylim((-1, 1))
ax.set_xlabel("Contrast")
ax.set_ylabel("Frequency in Hz")
ax.set_ylabel("Frequency [Hz]")
ax.set_xticks([-0.2, -0.1, 0, 0.1, 0.2])
ax.set_xlim((-0.21, 0.2))
ylim = ax.get_ylim()
@ -242,8 +242,8 @@ def fi_point_detection():
axes[0].set_xlim((-0.2, part - 0.2))
ylimits = axes[0].get_ylim()
axes[0].set_xlabel("Time in s")
axes[0].set_ylabel("Frequency in Hz")
axes[0].set_xlabel("Time [s]")
axes[0].set_ylabel("Frequency [Hz]")
axes[0].set_title("Step Response")
contrasts = fi.stimulus_values
@ -260,7 +260,7 @@ def fi_point_detection():
axes[1].plot(x_values, f_zero_fit, color=consts.COLOR_DATA_f0)
axes[1].plot(x_values, f_inf_fit, color=consts.COLOR_DATA_finf)
axes[1].set_xlabel("Contrast in %")
axes[1].set_xlabel("Contrast")
# axes[1].set_ylabel("Frequency in Hz")
axes[1].set_title("f-I Curve")
axes[1].set_ylim((0, ylimits[1]))

18
Figures_Model.py
View File

@ -34,7 +34,7 @@ def model_comparison():
axes[0, 0].plot(np.arange(-0.5, duration, step)[:len(stimulus)], stimulus)
axes[0, 0].set_title("Stimulus")
axes[0, 0].set_ylabel("V in mV")
axes[0, 0].set_ylabel("V [mV]")
axes[0, 1].set_frame_on(False)
axes[0, 1].set_axis_off()
axes[0, 0].set_ylim((-0.1, 1.5))
@ -48,8 +48,8 @@ def model_comparison():
time, freq = hF.calculate_time_and_frequency_trace(spikes, step)
axes[1, 1].plot(time, freq)
axes[1, 0].set_title("PIF")
axes[1, 0].set_ylabel("V in mV")
axes[1, 1].set_ylabel("Frequency in Hz")
axes[1, 0].set_ylabel("V [mV]")
axes[1, 1].set_ylabel("Frequency [Hz]")
v1, spikes = sM.lif_simulation(stimulus, step)
spikes = np.array(spikes)-0.5
@ -58,8 +58,8 @@ def model_comparison():
time, freq = hF.calculate_time_and_frequency_trace(spikes, step)
axes[2, 1].plot(time, freq)
axes[2, 0].set_title("LIF")
axes[2, 0].set_ylabel("V in mV")
axes[2, 1].set_ylabel("Frequency in Hz")
axes[2, 0].set_ylabel("V [mV]")
axes[2, 1].set_ylabel("Frequency [Hz]")
v1, spikes = sM.lifac_simulation(stimulus, step)
spikes = np.array(spikes) - 0.5
@ -68,11 +68,11 @@ def model_comparison():
time, freq = hF.calculate_time_and_frequency_trace(spikes, step)
axes[3, 1].plot(time, freq)
axes[3, 0].set_title("LIFAC")
axes[3, 0].set_ylabel("V in mV")
axes[3, 1].set_ylabel("Frequency in Hz")
axes[3, 0].set_ylabel("V [mV]")
axes[3, 1].set_ylabel("Frequency [Hz]")
axes[3, 0].set_xlabel("Time in s")
axes[3, 1].set_xlabel("Time in s")
axes[3, 0].set_xlabel("Time [s]")
axes[3, 1].set_xlabel("Time [s]")
# v1, spikes = sM.lifac_ref_simulation(stimulus, step)
# spikes = np.array(spikes) - 0.5
# axes[4, 0].plot(np.arange(-0.5, duration, step)[:len(v1)], v1)

6
Figures_Stimuli.py
View File

@ -16,7 +16,7 @@ def main():
def am_generation():
cell = "data/final/2013-04-17-ac-invivo-1"
cell_data = CellData(cell)
fig, axes = plt.subplots(3, 1, sharey=True, sharex=True, figsize=consts.FIG_SIZE_SMALL)
fig, axes = plt.subplots(3, 1, sharey=True, sharex=True, figsize=consts.FIG_SIZE_MEDIUM)
start = 0
end = 0.05
@ -47,11 +47,11 @@ def am_generation():
axes[1].set_title("Amplitude modulation")
axes[1].plot(time, values*am)
axes[1].plot(time, am)
axes[1].set_ylabel("Voltage in mV")
axes[1].set_ylabel("Voltage [mV]")
axes[2].set_title("Resulting stimulus")
axes[2].plot(time, (values*am) + values)
axes[2].set_xlabel("Time in ms")
axes[2].set_xlabel("Time [ms]")
axes[2].set_xlim((time[0], time[-1]))
axes[2].set_xticks([-25, 0, 25, 50, 75])
plt.tight_layout()

72
Figures_results.py
View File

@ -2,7 +2,7 @@ import numpy as np
import matplotlib.pyplot as plt
import matplotlib.gridspec as gridspec
import matplotlib as mpl
from analysis import get_filtered_fit_info, get_behaviour_values, get_parameter_values, behaviour_correlations, parameter_correlations
from analysis import get_filtered_fit_info, get_behaviour_values, get_parameter_values, behaviour_correlations, parameter_correlations, calculate_percent_errors
from ModelFit import get_best_fit
from Baseline import BaselineModel, BaselineCellData
from FiCurve import FICurveModel, FICurveCellData
@ -33,26 +33,27 @@ def main():
# run_all_images()
# quit()
dir_path = "results/final_2/"
# dend_tau_and_ref_effect()
# quit()
dir_path = "results/final_2/"
fits_info = get_filtered_fit_info(dir_path, filter=True)
print("Cells left:", len(fits_info))
cell_behaviour, model_behaviour = get_behaviour_values(fits_info)
plot_cell_model_comp_baseline(cell_behaviour, model_behaviour)
# cell_behaviour, model_behaviour = get_behaviour_values(fits_info)
# plot_cell_model_comp_baseline(cell_behaviour, model_behaviour)
# plot_cell_model_comp_burstiness(cell_behaviour, model_behaviour)
# plot_cell_model_comp_adaption(cell_behaviour, model_behaviour)
# behaviour_correlations_plot(fits_info)
# parameter_correlation_plot(fits_info)
create_parameter_distributions(get_parameter_values(fits_info))
parameter_correlation_plot(fits_info)
#
# create_parameter_distributions(get_parameter_values(fits_info))
create_parameter_distributions(get_parameter_values(fits_info, scaled=True, goal_eodf=800), "scaled_to_800_")
# errors = calculate_percent_errors(fits_info)
# create_boxplots(errors)
# example_good_hist_fits(dir_path)
# example_bad_hist_fits(dir_path)
# example_good_fi_fits(dir_path)
# example_bad_fi_fits(dir_path)
@ -82,7 +83,7 @@ def run_all_images():
def dend_tau_and_ref_effect():
cells = ["2012-12-21-am-invivo-1", "2014-03-19-ad-invivo-1", "2018-05-08-ac-invivo-1"]
cells = ["2012-12-21-am-invivo-1", "2014-03-19-ad-invivo-1", "2014-03-25-aa-invivo-1"]
cell_type = ["no burster", "burster", "strong burster"]
folders = ["results/ref_and_tau/no_dend_tau/", "results/ref_and_tau/no_ref_period/", "results/final_2/"]
title = [r"without $\tau_{dend}$", r"without $t_{ref}$", "with both"]
@ -94,6 +95,11 @@ def dend_tau_and_ref_effect():
cell_baseline = BaselineCellData(cell_data)
cell_baseline.load_values(cell_data.get_data_path())
eodf = cell_data.get_eod_frequency()
print(cell)
print("EODf:", eodf)
print("base rate:", cell_baseline.get_baseline_frequency())
print("bursty:", cell_baseline.get_burstiness())
print()
for j, folder in enumerate(folders):
fit = get_best_fit(folder + cell)
model_baseline = BaselineModel(fit.get_model(), eodf)
@ -101,12 +107,14 @@ def dend_tau_and_ref_effect():
model_isis = model_baseline.get_interspike_intervals() * eodf
bins = np.arange(0, 0.025, 0.0001) * eodf
if i == 0 and j == 2:
axes[i, j].hist(model_isis, density=True, bins=bins, color=consts.COLOR_MODEL, alpha=0.75,
label="model")
axes[i, j].hist(cell_isis, density=True, bins=bins, color=consts.COLOR_DATA, alpha=0.5, label="data")
axes[i, j].hist(model_isis, density=True, bins=bins, color=consts.COLOR_MODEL, alpha=0.75, label="model")
axes[i, j].legend(loc="upper right", frameon=False)
else:
axes[i, j].hist(cell_isis, density=True, bins=bins, color=consts.COLOR_DATA, alpha=0.5)
axes[i, j].hist(model_isis, density=True, bins=bins, color=consts.COLOR_MODEL, alpha=0.75)
axes[i, j].hist(cell_isis, density=True, bins=bins, color=consts.COLOR_DATA, alpha=0.5)
if j == 0:
axes[i, j].set_ylabel(cell_type[i])
axes[i, j].set_yticklabels([])
@ -114,7 +122,7 @@ def dend_tau_and_ref_effect():
axes[0, j].set_title(title[j])
plt.xlim(0, 17.5)
fig.text(0.5, 0.04, 'Time in EOD periods', ha='center', va='center') # shared x label
fig.text(0.5, 0.04, 'Time [EOD periods]', ha='center', va='center') # shared x label
fig.text(0.06, 0.5, 'ISI Density', ha='center', va='center', rotation='vertical') # shared y label
fig.text(0.11, 0.9, 'A', ha='center', va='center', rotation='horizontal', size=16, family='serif')
@ -206,7 +214,7 @@ def create_correlation_plot(ax, labels, correlations, p_values, title, y_label=T
else:
cleaned_cors[i, j] = np.NAN
if j >= i:
if j > i:
cleaned_cors[i, j] = np.NAN
im = ax.imshow(cleaned_cors, vmin=-1, vmax=1)
@ -231,9 +239,6 @@ def create_correlation_plot(ax, labels, correlations, p_values, title, y_label=T
# Loop over data dimensions and create text annotations.
for i in range(len(labels)):
for j in range(len(labels)):
if j >= i:
continue
if not np.isnan(cleaned_cors[i, j]):
if cleaned_cors[i, j] > 0:
text = ax.text(j, i, "{:.2f}".format(cleaned_cors[i, j]), ha="center", va="center", color="black", size=6)
@ -290,6 +295,8 @@ def example_bad_hist_fits(dir_path):
fig, axes = plt.subplots(1, 3, sharex="all", figsize=consts.FIG_SIZE_SMALL_EXTRA_WIDE) # , gridspec_kw={"top": 0.95})
for i, cell in enumerate([bursty_cell, strong_bursty_cell, extra_structure_cell]):
fit_dir = dir_path + cell + "/"
fit = get_best_fit(fit_dir)
@ -298,9 +305,16 @@ def example_bad_hist_fits(dir_path):
model = fit.get_model()
baseline_model = BaselineModel(model, eodf, trials=5)
cell_baseline = BaselineCellData(cell_data)
print(cell)
print("EODf:", eodf)
print("base rate:", cell_baseline.get_baseline_frequency())
print("bursty:", cell_baseline.get_burstiness())
print()
model_isi = np.array(baseline_model.get_interspike_intervals()) * eodf
cell_isi = BaselineCellData(cell_data).get_interspike_intervals() * eodf
cell_isi = cell_baseline.get_interspike_intervals() * eodf
bins = np.arange(0, 0.025, 0.0001) * eodf
if i == 0:
@ -311,7 +325,7 @@ def example_bad_hist_fits(dir_path):
axes[i].hist(model_isi, bins=bins, density=True, alpha=0.75, color=consts.COLOR_MODEL)
axes[i].hist(cell_isi, bins=bins, density=True, alpha=0.5, color=consts.COLOR_DATA)
axes[i].set_xlabel("ISI in EOD periods")
axes[i].set_xlabel("ISI [EOD periods]")
axes[0].set_ylabel("Density")
plt.tight_layout()
consts.set_figure_labels(xoffset=-2.5, yoffset=1.25)
@ -322,14 +336,24 @@ def example_bad_hist_fits(dir_path):
def example_good_fi_fits(dir_path):
fig, axes = plt.subplots(1, 3, figsize=consts.FIG_SIZE_SMALL_EXTRA_WIDE, sharey="all")
for i, cell in enumerate(["2012-12-21-am-invivo-1", "2013-02-21-ag-invivo-1", "2014-03-19-ae-invivo-1"]):
for i, cell in enumerate(["2012-12-21-am-invivo-1", "2014-03-19-ae-invivo-1", "2014-03-25-aa-invivo-1" ]):
fit_dir = dir_path + cell + "/"
fit = get_best_fit(fit_dir)
cell_data = fit.get_cell_data()
eodf = cell_data.get_eod_frequency()
cell_baseline = BaselineCellData(cell_data)
print(cell)
print("EODf:", eodf)
print("base rate:", cell_baseline.get_baseline_frequency())
print("bursty:", cell_baseline.get_burstiness())
print()
model = fit.get_model()
contrasts = cell_data.get_fi_contrasts()
fi_curve_data = FICurveCellData(cell_data, contrasts, save_dir=cell_data.get_data_path())
@ -357,6 +381,7 @@ def example_good_fi_fits(dir_path):
marker=consts.finf_marker, alpha=0.75, color=consts.COLOR_MODEL_finf, label=r"model $f_{\infty}$")
axes[i].set_xlabel("Contrast")
axes[i].set_xlim((-0.22, 0.22))
axes[0].legend(loc="upper left", frameon=False)
axes[0].set_ylabel("Frequency [Hz]")
@ -378,6 +403,14 @@ def example_bad_fi_fits(dir_path):
cell_data = fit.get_cell_data()
eodf = cell_data.get_eod_frequency()
cell_baseline = BaselineCellData(cell_data)
print(cell)
print("EODf:", eodf)
print("base rate:", cell_baseline.get_baseline_frequency())
print("bursty:", cell_baseline.get_burstiness())
print()
model = fit.get_model()
contrasts = cell_data.get_fi_contrasts()
fi_curve_data = FICurveCellData(cell_data, contrasts, save_dir=cell_data.get_data_path())
@ -405,6 +438,7 @@ def example_bad_fi_fits(dir_path):
marker=consts.finf_marker, alpha=0.75, color=consts.COLOR_MODEL_finf, label=r"model $f_{\infty}$")
axes[i].set_xlabel("Contrast")
axes[i].set_xlim((-0.22, 0.2))
axes[0].set_ylabel("Frequency [Hz]")
axes[0].legend(loc="upper left", frameon=False)

7
analysis.py
View File

@ -60,8 +60,11 @@ def get_filtered_fit_info(folder, filter=True):
# filter:
if filter:
if model_behaviour["f_zero_slope"] > 50000 or cell_behaviour["f_zero_slope"] > 50000:
print("f_zero_slope used to filter a fit, ", cell_folder)
if model_behaviour["f_zero_slope"] > 50000:
print("model bad: f_zero_slope used to filter a fit, ", cell_folder)
continue
if cell_behaviour["f_zero_slope"] > 50000:
print("BECAUSE OF CELL: f_zero_slope used to filter a fit, ", cell_folder)
continue
# if (abs(model_behaviour["f_inf_slope"] - cell_behaviour["f_inf_slope"]) / cell_behaviour["f_inf_slope"]) > 0.25:
# print("f_inf_slope used to filter a fit")

2
models/LIFACnoise.py
View File

@ -184,7 +184,7 @@ class LifacNoiseModel(AbstractModel):
def get_eodf_scaled_parameters(self, factor):
scaled_parameters = self.parameters.copy()
time_param_keys = ["refractory_period", "tau_a", "mem_tau", "dend_tau"]
time_param_keys = ["refractory_period", "tau_a", "mem_tau", "dend_tau", "delta_a"]
for key in time_param_keys:
scaled_parameters[key] = self.parameters[key] / factor

33
random_models.py
View File

@ -106,10 +106,11 @@ def compare_distribution_random_vs_fitted_params(par_list, scaled_param_values):
# if max_v > limit:
# print("For {} the max value was limited to {}, {} values were excluded!".format(l, limit, np.sum(np.array(cell_b_values[l]) > limit)))
# max_v = limit
step = (max_v - min_v) / 30
bins = np.arange(min_v, max_v + step, step)
axes_flat[i].hist(fitted_model_values, bins=bins, alpha=0.5, density=True)
axes_flat[i].hist(rand_model_values, bins=bins, alpha=0.5, density=True)
values = list(rand_model_values)
values.extend(fitted_model_values)
bins = calculate_bins(values, 30)
axes_flat[i].hist(fitted_model_values, bins=bins, color=consts.COLOR_MODEL, alpha=0.75, density=True)
axes_flat[i].hist(rand_model_values, bins=bins, color="black", alpha=0.5, density=True)
axes_flat[i].set_xlabel(parameter_titles[l] + " " + x_labels[i])
axes_flat[i].set_yticks([])
axes_flat[i].set_yticklabels([])
@ -127,7 +128,7 @@ def compare_distribution_random_vs_fitted_params(par_list, scaled_param_values):
def parameter_correlation_plot(par_list, fits_info):
fig = plt.figure(tight_layout=True, figsize=consts.FIG_SIZE_MEDIUM_WIDE)
gs = gridspec.GridSpec(2, 2, width_ratios=(1, 1), height_ratios=(5, 1), hspace=0.05, wspace=0.05)
gs = gridspec.GridSpec(2, 2, width_ratios=(1, 1), height_ratios=(5, 1), hspace=0.3, wspace=0.05)
# fig, axes = plt.subplots(1, 2, figsize=consts.FIG_SIZE_MEDIUM_WIDE)
labels, corr_values, corrected_p_values = parameter_correlations(fits_info)
@ -150,8 +151,6 @@ def parameter_correlation_plot(par_list, fits_info):
ax_col.imshow(data)
ax_col.set_xlabel("Correlation Coefficients")
plt.savefig(consts.SAVE_FOLDER + "rand_parameter_correlations_comparison.pdf")
plt.close()
@ -258,21 +257,14 @@ def load_behavior():
def create_behaviour_distributions(drawn_model_behaviour, fits_info):
fig, axes = plt.subplots(4, 2, gridspec_kw={"left": 0.1, "hspace":0.5}, figsize=consts.FIG_SIZE_LARGE_HIGH)
cell_behaviour, fitted_model_behaviour = get_behaviour_values(fits_info)
labels = ['Burstiness', 'baseline_frequency', 'coefficient_of_variation', 'f_inf_slope', 'f_zero_slope', 'serial_correlation', 'vector_strength']
unit = ["[%ms]", "[Hz]", "", "[Hz]", "[Hz]", "", ""]
labels = ['baseline_frequency', 'serial_correlation', 'vector_strength', 'Burstiness', 'coefficient_of_variation', 'f_inf_slope', 'f_zero_slope']
unit = ["[Hz]", "", "", "[%ms]", "", "[Hz]", "[Hz]"]
axes_flat = axes.flatten()
for i, l in enumerate(labels):
min_v = min(drawn_model_behaviour[l]) * 0.95
max_v = max(drawn_model_behaviour[l]) * 1.05
# limit = 50000
# if max_v > limit:
# print("For {} the max value was limited to {}, {} values were excluded!".format(l, limit, np.sum(np.array(cell_b_values[l]) > limit)))
# max_v = limit
step = (max_v - min_v) / 20
bins = np.arange(min_v, max_v + step, step)
axes_flat[i].hist(drawn_model_behaviour[l], bins=bins, alpha=0.75, density=True, color=consts.COLOR_MODEL)
axes_flat[i].hist(cell_behaviour[l], bins=bins, alpha=0.5, density=True, color=consts.COLOR_DATA)
bins = calculate_bins(drawn_model_behaviour[l], 20)
axes_flat[i].hist(drawn_model_behaviour[l], bins=bins, density=True, color=consts.COLOR_MODEL, alpha=0.75)
axes_flat[i].hist(cell_behaviour[l], bins=bins, density=True, color=consts.COLOR_DATA, alpha=0.5)
axes_flat[i].set_xlabel(behaviour_titles[l] + " " + unit[i])
axes_flat[i].set_yticks([])
axes_flat[i].set_yticklabels([])
@ -362,7 +354,7 @@ def plot_distributions_with_set_fits(param_values):
fig, axes = plt.subplots(4, 2, gridspec_kw={"left": 0.1, "hspace":0.5}, figsize=consts.FIG_SIZE_LARGE_HIGH)
gauss_fits = get_gauss_fits()
bin_number = 30
bin_number = 20
labels = ["input_scaling", "v_offset", "mem_tau", "noise_strength",
"tau_a", "delta_a", "dend_tau", "refractory_period"]
@ -417,7 +409,6 @@ def plot_distributions(param_values):
# normal hist:
values = param_values[key]
bins = calculate_bins(values, bin_number)
normal, n_bins, patches = axes[i, 0].hist(values, bins=calculate_bins(values, bin_number), density=True)
axes[i, 0].set_title(key)

12
read_data_infos.py
View File

@ -12,11 +12,9 @@ def main():
def eod_info():
cells = []
for item in sorted(os.listdir("data/invivo/")):
cells.append(os.path.join("data/invivo/", item))
for item in sorted(os.listdir("data/final/")):
cells.append(os.path.join("data/final/", item))
for item in sorted(os.listdir("data/invivo_bursty/")):
cells.append(os.path.join("data/invivo_bursty/", item))
eod_freq = []
@ -33,11 +31,9 @@ def eod_info():
def fish_info():
cells = []
for item in sorted(os.listdir("data/invivo/")):
cells.append(os.path.join("data/invivo/", item))
for item in sorted(os.listdir("data/final/")):
cells.append(os.path.join("data/final/", item))
for item in sorted(os.listdir("data/invivo_bursty/")):
cells.append(os.path.join("data/invivo_bursty/", item))
cell_type = []
weight = []
size = []

6
run_Fitter.py
View File

@ -14,8 +14,8 @@ from helperFunctions import plot_errors
import multiprocessing as mp
SAVE_DIRECTORY = "./results/final_3"
SAVE_DIRECTORY_BEST = "./results/final_3_best/"
SAVE_DIRECTORY = "./results/ref_and_tau/no_ref_period/"
SAVE_DIRECTORY_BEST = "./results/ref_and_tau/nrp_best/"
# SAVE_DIRECTORY = "./results/ref_and_tau/no_dend_tau/"
# SAVE_DIRECTORY_BEST = "./results/ref_and_tau/ndt_best/"
@ -68,7 +68,7 @@ def fit_cell_base(parameters):
time1 = time.time()
fitter = Fitter()
fitter.set_data_reference_values(parameters[0])
fmin, res_par = fitter.fit_routine(parameters[2], ERROR_WEIGHTS)
fmin, res_par = fitter.fit_routine_no_ref_period(parameters[2], ERROR_WEIGHTS)
cell_data = parameters[0]
cell_name = os.path.split(cell_data.get_data_path())[-1]

2
run_cellwise_fitting.sh
View File

@ -1,5 +1,5 @@
for file in data/final/*; do
for file in data/ref_and_dend_set/*; do
if [ -d "$file" ]; then
nice python3 run_Fitter.py --cell $file
fi

21
thesis/Masterthesis.bbl
View File

@ -33,6 +33,11 @@ Chacron, M.~J., Longtin, A., and Maler, L. (2001).
\newblock Simple models of bursting and non-bursting p-type electroreceptors.
\newblock {\em Neurocomputing}, 38:129--139.
\bibitem[Chacron et~al., 2004]{chacron2004burst}
Chacron, M.~J., Longtin, A., and Maler, L. (2004).
\newblock To burst or not to burst?
\newblock {\em Journal of computational neuroscience}, 17(2):127--136.
\bibitem[Gao and Han, 2012]{gao2012implementing}
Gao, F. and Han, L. (2012).
\newblock Implementing the nelder-mead simplex algorithm with adaptive
@ -56,6 +61,11 @@ Maciver, M.~A., Sharabash, N.~M., and Nelson, M.~E. (2001).
effects of water conductivity.
\newblock {\em Journal of experimental biology}, 204(3):543--557.
\bibitem[Nelder and Mead, 1965]{nelder1965simplex}
Nelder, J.~A. and Mead, R. (1965).
\newblock A simplex method for function minimization.
\newblock {\em The computer journal}, 7(4):308--313.
\bibitem[Olypher and Calabrese, 2007]{olypher2007using}
Olypher, A.~V. and Calabrese, R.~L. (2007).
\newblock Using constraints on neuronal activity to reveal compensatory changes
@ -103,6 +113,17 @@ Walz, H., Grewe, J., and Benda, J. (2014).
evoked by transient communication signals.
\newblock {\em Journal of Neurophysiology}, 112(4):752--765.
\bibitem[Xu et~al., 1996]{xu1996logarithmic}
Xu, Z., Payne, J.~R., and Nelson, M.~E. (1996).
\newblock Logarithmic time course of sensory adaptation in electrosensory
afferent nerve fibers in a weakly electric fish.
\newblock {\em Journal of neurophysiology}, 76(3):2020--2032.
\bibitem[Zeldenrust et~al., 2018]{zeldenrust2018neural}
Zeldenrust, F., Wadman, W.~J., and Englitz, B. (2018).
\newblock Neural coding with bursts—current state and future perspectives.
\newblock {\em Frontiers in computational neuroscience}, 12:48.
\bibitem[Zupanc et~al., 2006]{zupanc2006electric}
Zupanc, G., S{\^\i}rbulescu, R., Nichols, A., and Ilies, I. (2006).
\newblock Electric interactions through chirping behavior in the weakly

60
thesis/Masterthesis.blg
View File

@ -3,44 +3,44 @@ Capacity: max_strings=35307, hash_size=35307, hash_prime=30011
The top-level auxiliary file: Masterthesis.aux
The style file: apalike.bst
Database file #1: citations.bib
You've used 19 entries,
You've used 23 entries,
1935 wiz_defined-function locations,
606 strings with 7343 characters,
and the built_in function-call counts, 7914 in all, are:
= -- 800
> -- 300
< -- 8
+ -- 100
- -- 98
* -- 731
:= -- 1365
add.period$ -- 57
call.type$ -- 19
change.case$ -- 142
chr.to.int$ -- 19
cite$ -- 19
duplicate$ -- 264
empty$ -- 566
format.name$ -- 128
if$ -- 1522
632 strings with 7966 characters,
and the built_in function-call counts, 9564 in all, are:
= -- 964
> -- 367
< -- 9
+ -- 124
- -- 120
* -- 884
:= -- 1651
add.period$ -- 69
call.type$ -- 23
change.case$ -- 173
chr.to.int$ -- 23
cite$ -- 23
duplicate$ -- 320
empty$ -- 682
format.name$ -- 156
if$ -- 1838
int.to.chr$ -- 1
int.to.str$ -- 0
missing$ -- 18
newline$ -- 98
num.names$ -- 57
pop$ -- 98
missing$ -- 22
newline$ -- 118
num.names$ -- 69
pop$ -- 120
preamble$ -- 1
purify$ -- 143
purify$ -- 174
quote$ -- 0
skip$ -- 199
skip$ -- 240
stack$ -- 0
substring$ -- 706
swap$ -- 19
substring$ -- 843
swap$ -- 23
text.length$ -- 0
text.prefix$ -- 0
top$ -- 0
type$ -- 114
type$ -- 138
warning$ -- 0
while$ -- 73
while$ -- 88
width$ -- 0
write$ -- 249
write$ -- 301

BIN
thesis/Masterthesis.pdf
View File

Binary file not shown.

246
thesis/Masterthesis.tex
View File

@ -97,10 +97,10 @@ Außerdem erkläre ich, dass die eingereichte Arbeit weder vollständig noch in
\item update the colors in all plots to be consistent.
\item make plot labels consistent (Units: in mV vs [mV])
\item update number of cells / fish etc
\item check time form in text results present!
\end{itemize}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
@ -121,40 +121,17 @@ Außerdem erkläre ich, dass die eingereichte Arbeit weder vollständig noch in
% Einleitung
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Introduction}
\begin{comment}
\begin{enumerate}
\item sensory input important for all life etc.
\item electric fish
\begin{enumerate}
\item general: habitat,
\item as model animal for ethology
\item electric organ + eod
\item sensory neurons p- and t(?)-type
\end{enumerate}
\item sensory perception
\begin{enumerate}
\item receptor $\rightarrow$ heterogenic population
\item further analysis limited by what receptors code for - P-Units encoding
\item p-type neurons code AMs
\end{enumerate}
\item goal be able to simulate heterogenic population to analyze full coding properties $\rightarrow$ many cells at the same time needed $\rightarrow$ only possible in vitro/ with model simulations
\item Possible to draw representative values for model parameters to generate a population ?
\end{enumerate}
\newpage
\end{comment}
\todo{A complete characterisation of their activity has been the sub-
ject of previous studies (Gabbiani et al., 1996; Gussin et al., 2007; Scheich et al., 1973;
Wessel et al., 1996; Xu et al., 1996; Zakon, 1986) ???}
The environment of an organism holds important information that it needs to survive. Information about predators to avoid, food to find and potential mates. The ability to sense and process this information is of vital importance for any organism. At the same time the environment also contains a lot of information that is irrelevant to an organism. \cite{barlow1961possible} suggested already that the sensory systems of an organism should be specialized to extract the information it needs while filtering out the noise and irrelevant information, to efficiently use the limited coding capacity of the sensory systems.
One interesting model system for questions adaptive signal processing is the electric fish \AptLepto (Brown ghost knife fish).
\lepto generate a sinusoidal electric field with the electric organ in their tail enabling them to use active electroreception which they use to find prey and communicate with each other (\cite{maciver2001prey}, \cite{zupanc2006electric}). The different use cases of this electric organ discharge (EOD) come with the necessity to detect a wide range of different amplitude modulations (AMs). Electrolocation of object in the surrounding water like small prey or rocks cause small low frequency AMs \citep{babineau2007spatial}. At the same time other electric fish can cause stronger and higher frequency AMs through interference between the electric fields and their communication signals like chirps, short increases in their EOD frequency \citep{zupanc2006electric}. This means that the electroreceptors need to be able to encode a wide range of changes in EOD amplitude, in speed as well as strength.
The EOD and its AMs are encoded by electroreceptor organs in the skin. \lepto have two kinds of tuberous electrosensory organs: the T and P type units \citep{scheich1973coding}. The T units (time coder) are strongly phase locked to the EOD and fire regularly once every EOD period. They encode the phase of the EOD in their spike timing. The P units (probability coders) on the other hand do not fire every EOD period. Instead they fire irregularly with a certain probability that depends on the EOD amplitude. That way they encode information about the EOD amplitude in their firing probability \citep{scheich1973coding}. An example of the firing behavior of a P unit is shown in figure~\ref{fig:p_unit_example}.
When the fish's EOD is unperturbed P units fire every few EOD periods but they have a certain variability in their firing (fig. \ref{fig:p_unit_example} B) and show negative correlation between successive interspike intervals (ISIs)(fig. \ref{fig:p_unit_example} C). When presented with a step increase in EOD amplitude P units show strong adaption behavior. After a strong increase in firing rate reacting to the onset of the step, the firing rate quickly decays back to a steady state (fig. \ref{fig:p_unit_example} D). When using different sizes of steps both the onset and the steady state response scale with its size and direction of the step (fig. \ref{fig:p_unit_example} E).
When the fish's EOD is unperturbed P units fire every few EOD periods but they have a certain variability in their firing (fig. \ref{fig:p_unit_example} B) and show negative correlation between successive interspike intervals (ISIs)(fig. \ref{fig:p_unit_example} C). When presented with a step increase in EOD amplitude P units show strong adaption behavior. After a strong increase in firing rate reacting to the onset of the step, the firing rate quickly decays back to a steady state (fig. \ref{fig:p_unit_example} D). When using different sizes of steps both the onset and the steady state response scale with its size and direction of the step (fig. \ref{fig:p_unit_example} E).
\begin{figure}[H]
@ -170,38 +147,26 @@ When the fish's EOD is unperturbed P units fire every few EOD periods but they h
{\includegraphics{figures/isi_hist_heterogeneity.pdf}}
\end{figure}
\todo{heterogeneity more, bursts important for coding in other systems}
Furthermore show P-units a pronounced heterogeneity in their spiking behavior (fig.~\ref{fig:heterogeneity_isi_hist}, \cite{gussin2007limits}). This is an important aspect one needs to consider when trying to understand what and how information is encoded in the spike trains of the neuron (\cite{padmanabhan2010intrinsic}, \cite{tripathy2013intermediate}). A single neuron might be an independent unit from all other neurons but through different tuning curves a full picture of the stimulus can be encoded in the population even when a single neuron only encodes a small feature space. This type of encoding is ubiquitous in the nervous system and is used in the visual sense for color vision, PLUS MORE... \todo{refs labeled line vs summation code}. Even though P units were already modeled based on a simple leaky integrate-and-fire neuron \citep{chacron2001simple} and conductance based \citep{kashimori1996model} and well studied (\cite{bastian1981electrolocation}, \cite{ratnam2000nonrenewal} \cite{benda2005spike}). Up to this point there is no model that tries to cover the full breadth of heterogeneity of the P unit population. Having such a model could help shed light into the population code used in the electric sense and allow researchers gain a better picture how higher brain areas might process the information and get one step closer to the full path between sensory input and behavioral output.
\todo{viel wichtiger: das ist das perfekte modellsystem um die Kodierung in einer heterogenen population veon neuronen zu untersuchen. Relativ unabhaengig vom efish.}
Furthermore show P-units a pronounced heterogeneity in their spiking behavior (fig.~\ref{fig:heterogeneity_isi_hist}, \cite{gussin2007limits}). Currently the spiking behavior of P-units is often split into two distinct categories of bursting and non-bursting cells (\cite{xu1996logarithmic}, \cite{chacron2004burst}) or the bursting behavior not considered at all \citep{walz2013Phd}, but when one is trying to understand how information is encoded in the spike trains and populations of neurons the bursts (review: \cite{zeldenrust2018neural}) and general heterogeneity (\cite{padmanabhan2010intrinsic}, \cite{tripathy2013intermediate}) are important aspects to consider.
A single neuron might be an independent unit from all other neurons but through different tuning curves a full picture of the stimulus can be encoded in the population even when a single neuron only encodes a small feature space. This type of encoding is ubiquitous in the nervous system and is,for example in form of a labeled line code, used in the visual sense for color vision. Even though P-units were already modeled based on a simple leaky integrate-and-fire neuron \citep{chacron2001simple} and conductance based \citep{kashimori1996model} and well studied (\cite{bastian1981electrolocation}, \cite{ratnam2000nonrenewal} \cite{benda2005spike}). Up to this point there is no model that tries to cover the full breadth of heterogeneity of the P unit population. Having such a model could help shed light into the population code used in the electric sense, but also of heterogeneous neuron populations in general as there are currently few model systems that have well defined heterogeneous populations.
Further it could allow researchers gain a better picture how higher brain areas might process the information and get closer to the full path between sensory input and behavioral output.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Methoden
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Materials and Methods}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\todo{some transition from the introduction}
\section{Materials and Methods}
\subsection{Cell recordings}
%sampling: 100000.0 count: 20
%sampling: 20000.0 count: 54
%sampling: 40000.0 count: 1
% each fish only once in calculation (as determined with max one fish per day)
% # of fish = 32
% EOD-freq: min 601.09, mean 753.09, max 928.45, std 82.30
% Sizes: min 11.00, mean 15.78, max 25.00, std 3.48
The cell recordings for this master thesis were collected as part of other previous studies (\cite{walz2013Phd}, \citep{walz2014static})\todo{ref other studies} and the recording procedure is described there but will also be repeated below. The recordings of altogether 457 p-units were inspected. Of those 88 fulfilled basic necessary requirements: including a measurement of at least 30 seconds of baseline behavior and containing at least 7 different contrasts with each at least 7 trials for the f-I curve (see below fig. \ref{fig:f_point_detection} B). After pre-analysis of those cells an additional 15 cells were excluded because of spike detection difficulties.
The cell recordings for this master thesis were collected as part of other previous studies (\cite{walz2013Phd}, \cite{walz2014static})\todo{ref other studies} and the recording procedure is described there but will also be repeated below. The recordings of altogether 457 p-units were inspected. Of those 88 fulfilled basic necessary requirements: including a measurement of at least 30 seconds of baseline behavior and containing at least 7 different contrasts with each at least 7 trials for the f-I curve (see below fig. \ref{fig:f_point_detection} B). After pre-analysis of those cells an additional 15 cells were excluded because of spike detection difficulties.
The 73 used cells came from 32 \AptLepto (brown ghost knifefish). The fish were between 11--25\,cm long (15.8 $\pm$ 3.5\,cm) and their electric organ discharge (EOD) frequencies ranged between 601 and 928\,Hz (753 $\pm$ 82\,Hz). The sex of the fish was not determined.
The 67 used cells came from 32 \AptLepto (brown ghost knifefish). The fish were between 11--25\,cm long (15.8 $\pm$ 3.5\,cm) and their electric organ discharge (EOD) frequencies ranged between 601 and 928\,Hz (753 $\pm$ 82\,Hz). The sex of the fish was not determined.
The in vivo intracellular recordings of P-unit electroreceptors were done in the lateral line nerve. The fish were anesthetized with MS-222 (100-130 mg/l; PharmaQ; Fordingbridge, UK) and the part of the skin covering the lateral line just behind the skull was removed, while the area was anesthetized with Lidocaine (2\%; bela-pharm; Vechta, Germany). The fish were immobilized for the recordings with Tubocurarine (Sigma-Aldrich; Steinheim, Germany, 25--50\,$\mu l$ of 5\,mg/ml solution) and placed in the experimental tank (47 $\times$ 42 $\times$ 12\,cm) filled with water from the fish's home tank with a conductivity of about 300$\mu$\,S/cm and the temperature was around 28°C.
All experimental protocols were approved and complied with national and regional laws (files: no. 55.2-1-54-2531-135-09 and Regierungspräsidium Tübingen no. ZP 1/13 and no. ZP 1/16)
@ -308,7 +273,7 @@ f_{\infty}(I) = \lfloor mI+c \rfloor_0
\subsection{Leaky Integrate and Fire Model}
\subsection{Leaky Integrate-and-Fire Model}
% add info about simulation by euler integration and which time steps!
% show voltage dynamics with resistance :
@ -420,7 +385,7 @@ For the $f_{inf}$ and $f_0$ responses the average scaled difference off all cont
err_i = c_i (\langle (x^M_i - x^C_i)²\rangle)
\end{equation}
All errors were then summed up for the full error. The fits were done with the Nelder-Mead algorithm of scipy minimize \citep{gao2012implementing}. All model variables listed above in table \ref{tab:parameter_explanation} were fit at the same time except for $I_{Bias}$. $I_{Bias}$ was determined before each fitting iteration and set to a value giving the correct baseline frequency within 2\,Hz.
All errors were then summed up for the full error. The fits were done with the simplex algorithm from \cite{nelder1965simplex} implemented in the python package Scipy according to \cite{gao2012implementing}. All model variables listed above in table \ref{tab:parameter_explanation} were fit at the same time except for $I_{Bias}$. $I_{Bias}$ was determined before each fitting iteration and set to a value giving the correct baseline frequency within 2\,Hz.
\begin{table}[H]
@ -428,13 +393,13 @@ All errors were then summed up for the full error. The fits were done with the N
behavior & scaling factor \\
\hline
vector strength & 100 \\
coefficient of variation & 20 \\
serial correlation & 10 \\
ISI-histogram & 1/600\\
$f_0$ detections & 0.1 \\
$f_{\infty}$ detections & 1 \\
$f_\infty$ slope & 20 \\
vector strength & 100 \\
coefficient of variation & 20 \\
serial correlation & 10 \\
ISI-histogram & 1/600 \\
$f_0$ detections & 0.1 \\
$f_{\infty}$ detections & 1 \\
$f_\infty$ slope & 20 \\
$f_0$ step response & 0.001
\end{tabular}
\caption{\label{tab:scaling_factors} Scaling factors for fitting errors.}
@ -442,97 +407,162 @@ All errors were then summed up for the full error. The fits were done with the N
\todo{Fitting more in detail number of start parameters the start parameters themselves}
\todo{explain removal of bad models! Filter criteria how many filtered etc.}
\todo{explain how to draw random models from fitted model parameters! here probably?}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% RESULTS
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\newpage
\section{Results}
The here proposed model was, as described above, fit to each in vivo recording of a single P-unit with 12 different start parameters, choosing the best of the 12 resulting fitted models. The resulting 67 fitted models were then filtered for $f_0$ slope above 50000 and a CV more than 33\% different from the cell value. This filtering removed 13 cell-model pairs. 11 were filtered by the CV requirement and 2 by the $f_0$ slope requirement. This left 54 fitted models for further analysis.
\subsection{Model Examples}
First the effect of the dendritic low pass filter and the refractory period was investigated. This was done by fitting the model to three cells representing the continuously firing, weakly and strongly bursting cells, while removing the parameter to be investigated e.g. once without the low pass filter and once without a refractory period (fig. \ref{fig:dend_ref_effect}).
The model without the dendritic low pass filter was able to fit the shape of the ISI histogram even for the strongly bursting cell but was not able to correctly match the height of the distribution. It has very thin peaks which show that this model locks very strongly to the phase of the EOD and cannot match the weaker locking of the cells (fig. \ref{fig:dend_ref_effect} A).
When the refractory period of the model is disabled it can still match ISI histogram of the continuously firing cell perfectly but is not able to match the ISI histogram shape for the bursting cells. These cell have more than one local maximum in their histogram. One at the first EOD period caused by the bursts and a second maximum at a latter EOD period showing the times between bursts. The model without $t_{ref}$ is only able to produce a single maximum in their ISI histogram and cannot match the high firing probability at the first EOD period (fig. \ref{fig:dend_ref_effect} B).
The full model with both parameters can match both the shape and the height of the full ISI histogram for all three cells shown in figure \ref{fig:dend_ref_effect} C, but there are also cases in which the model fails to reproduce the ISI histogram of the cell (fig. \ref{fig:example_bad_isi_fits}). The cell A in figure \ref{fig:example_bad_isi_fits} shows a very strongly bursting cell. This cell has a very high peak at the first EOD period and there is some distance to the rest of the ISI distribution. This means that the cell has a few EOD periods in which it doesn't fire after a burst (compare fig. \ref{fig:heterogeneity_isi_hist} C). The fitted model fails at reproducing the long pauses between the bursts and in this case also shows a too low phase locking. In B the cell shows a high firing probability for the first two EOD periods after a spike and only after that shows the normal distribution \todo{explain} of firing probabilities. The fitted model does not show the high probability for the first two EOD periods. Instead it has only a high firing probability at the first EOD and a very low probability at the second EOD.
The last cell shown in figure \ref{fig:example_bad_isi_fits} has a higher order structure in its ISI histogram. It has high firing probability only at every second EOD period starting at the fourth EOD. This higher order structure is not matched by the model instead it shows a continuously increasing firing probability for each EOD period up to the maximum and then decreases again without being reduced every second EOD.
Comparing the f-I curves of the same three cells as in figure \ref{fig:dend_ref_effect} with the f-I curves of their fits shows good agreement for all three (fig. \ref{fig:example_good_fi_fits}), except the steady state response of cell C where the model shows a steeper slope. Note the miss-detections of the $f_0$ response in the negative contrasts of cell C influencing the fit of the Boltzmann function. The cells also demonstrate the variability of the cells in the strength of their response to the steps stimuli.
All fitted models were again inspected and representative failure cases are shown in figure \ref{fig:example_bad_fi_fits}. The fitted model in A overestimates the slope of the steady state response and underestimates the frequency in the lower half of the onset response.
The second example of the problematic cases (fig. \ref{fig:example_bad_fi_fits} B) is a cell with a very high baseline firing rate where the model does not manage to reproduce the plateau the cell quickly reaches for positive contrasts. Instead its the frequency of its onset response continues to increase above the possible range of one spike per EOD period.
\begin{figure}[H]
\centering
\includegraphics[trim={5mm 0mm 0mm 0mm}]{figures/dend_ref_effect.pdf}
\caption{\label{fig:dend_ref_effect} Effect of the dendritic filter and the refractory period on baseline firing. In each row data (blue) and model fits (orange) to three example cells are shown that differ in their burstiness as indicated on the left. \todo{Top: cell 2012-xxx, r=z Hz, b=y; center: cell 2014-xxx, r= b=...; bottom ...}
A: Without dendritic filter ($\tau_{dend}$) the spikes are too strongly locked to the EOD, resulting in very high vector strength and too narrow peaks in the baseline ISIH. B: Without refactory period ($t_{ref}$) the model cannot capture burstiness. While this is no problem for the non-bursting cell (top), the peak in the ISIH at one EOD period cannot be reproduced without refacory period. C: with both the dendritic filter and the refracory period ISIhs can be faithfully reproduced for all three cells.
%Effect of the addition of $\tau_{dend}$ and $t_{ref}$ to the model. Rows \textbf{1--3} different cells: not bursting, bursting or strongly bursting in \todo{color} the cell and in \todo{color} the model. The fits in each column (\textbf{A--C}) were done with different parameters to show their effect on the model. Column \textbf{A}: The cells were fit without $\tau_{dend}$. This causes the model to be unable to fit the vector strength correctly and the models are too strongly locked to the EOD phase. \textbf{B}: The models were fit without $t_{ref}$, because of that the model cannot match the burstiness of the cell. Visible in the missing high peak at the first EOD period. In column \textbf{C} the model all parameters. It can match the full spiking behavior of the cells for the different strengths of bursting.
}
\caption{\label{fig:dend_ref_effect} Effect of the dendritic filter and the refractory period on baseline firing. In each row data (blue) and model fits (orange) to three example cells are shown that differ in their burstiness as indicated on the left. Top: cell 2012-12-21-am, r=135\,Hz, b=0.02\,\%ms; center: cell 2014-03-19-ad-invivo-1, r=237\,Hz b=1.69\,\%ms; bottom: cell 2014-03-25-aa r=204\,Hz , b=1.9\,\%ms
A: Without dendritic filter ($\tau_{dend}$) the spikes are too strongly locked to the EOD, resulting in very high vector strength and too narrow peaks in the baseline ISI histogram. B: Without refractory period ($t_{ref}$) the model cannot capture the burstiness. While this is no problem for the non-bursting cell (top), the peak in the ISI histogram at one EOD period cannot be reproduced without refractory period. C: with both the dendritic filter and the refractory period ISI histograms can be faithfully reproduced for all three cells.}
\end{figure}
\begin{figure}[H]
\includegraphics{figures/example_bad_isi_hist_fits.pdf}
\caption{\label{fig:example_bad_isi_fits} \todo{Add pointer arrows in plot?} Problem cases in which the model ISI histogram wasn't fit correctly to the cell. ISI histograms of different cells (\todo{color}) and their corresponding model (\todo{color}). \textbf{A}: Strongly bursting cell with large pauses between bursts, where the model doesn't manage to reproduce the long pauses. \textbf{B}: Bursting cell with a high probability of firing in the first and second following EOD period. Here the model can't reproduce the high probability on the second following EOD period. \textbf{C}: Cell with a higher order structure \todo{??} in its ISI histogram. It only has a high firing probability every second EOD period which is also not represented in the model.}
\caption{\label{fig:example_bad_isi_fits} Problem cases in which the model ISI histogram wasn't fit correctly. \textbf{A}: (cell 2014-06-06-ag r=117\,Hz , b=3.9) Strongly bursting cell with large pauses between bursts, where the model doesn't manage to reproduce the long pauses. \textbf{B}: (cell 2018-05-08-ab r=112\,Hz , b=2.8) Bursting cell with a high probability of firing in the first and second following EOD period. Here the model can't reproduce the high probability on the second following EOD period. \textbf{C}: (cell 2014-12-11-ad rate=50\,Hz , b=0) Cell with a higher order structure in its ISI histogram. It only has a high firing probability every second EOD period which is also not represented in the model.}
\end{figure}
\begin{figure}[H]
\includegraphics{figures/example_good_fi_fits.pdf}
\caption{\label{fig:example_good_fi_fits} Good fit examples of the f-I curve. \textbf{A--C}: Three cells with different response patterns which are all well matched by their models. \todo{Color explanation} }
\caption{\label{fig:example_good_fi_fits} Good fit examples of the f-I curves. The red line fitted Boltzmann function and blue line fitted line for the cell's $f_0$ and $f_\infty$ response respectively. A: cell 2012-12-21-am r=135\,Hz, EODf=806\,Hz; B: cell 2014-03-19-ad-invivo-1 r=237\,Hz, EODf=658\,Hz; C: cell 2014-03-25-aa r=204\,Hz, EODf=870\,Hz. The cells show different response strengths to the contrasts. Which are all well matched by their models.}
\end{figure}
\begin{figure}[H]
\includegraphics{figures/example_bad_fi_fits.pdf}
\caption{\label{fig:example_bad_fi_fits} Examples of bad fits of the f-I curve. \textbf{A--C}: Different cells. \todo{Color explanation}. \textbf{A}: model that did not fit the negative contrast responses of the $f_0$ response well but was successful in the positive half. It also was not successful in the $f_\infty$ response and shows a wrong slope. \textbf{B}: A fit that was successful for the lower $f_0$ response but overshoots the limit of the cell and fires too fast for high positive contrasts. It also has a slightly wrong $f_\infty$ response slope.}
\caption{\label{fig:example_bad_fi_fits} Examples of bad fits of the f-I curves. The red line fitted Boltzmann function and blue line fitted line for the cell's $f_0$ and $f_\infty$ response respectively. \textbf{A}: (cell 2012-12-13-ao r=146\,Hz, EODf=657\,Hz) Model that did not fit the negative contrast responses of the $f_0$ response well but was successful in the positive half. It also was not successful in the $f_\infty$ response and shows a too steep slope. \textbf{B}: (cell 2014-01-23-ab r=431\,Hz, EODf=775\,Hz) A fit that was successful for the lower $f_0$ response but overshoots the limit of one spike per EOD period. It also has a slightly too steep $f_\infty$ response slope.}
\end{figure}
\begin{figure}[H]
\subsection{Population Comparison}
The general fitting quality was inspected by comparing the distributions in firing behavior between cells and fitted models as well as directly comparing each cell and its respective model. The baseline rate was matched perfectly (fig. \ref{fig:comp_baseline} A) because it was set to be equal within 2\,Hz during the fitting procedure by adjusting the bias current $I_{Bias}$ appropriately. Its approximately log-normal distribution was in the expected range of around 50--400\,Hz of the literature \citep{gussin2007limits}.
The vector strength (VS) was matched quiet well. Many cells have a VS of around 0.85 and with a few cells having a VS as low as 0.5. The fitted models show the same range of VSs but for cells with a VS above 0.8, the models often underestimating the true VS.
The models failed too fit the full breadth of the serial correlation shown by the cells, that have serial correlations between -0.8-- -0.1. The fits fail to match the strong negative SCs and reducing the range on the lower end to -0.7. The fits of the SC also show a high variability.
\begin{figure}[H]
\makebox[\textwidth][c]{\includegraphics{figures/fit_baseline_comparison.pdf}}
\caption{\label{fig:comp_baseline} Comparison of the cell behavior and the behavior of the corresponding fit: \textbf{A} baseline firing rate, \textbf{B} vector strength (VS) and \textbf{C} serial correlation (SC). The histograms compare the distributions of the cell (\todo{color}) and the model (\todo{color}). Below \todo{what is this plot called} In grey the line on which cell and model values are equal. \textbf{A}: The baseline firing rate of the cell and the model. The base rate agrees near perfectly as it is set to be equal within a margin of 2\,Hz during the fitting process. \textbf{B}: The vector strength agrees well for most cells but if the cells have a vector strength above 0.8 the probability increases for the models to show often a weaker VS than the cell. \textbf{C}: Comparison of the SC with lag 1. Here the models cluster more strongly and don't show quite the same range like the cells do. Models of cells with a strongly negative SC often have a weaker negative SC while the models in the opposite case show too strong negative correlations. In general is the fitting of the SC a lot more variable than the precise fitting of the VS.}
\caption{\label{fig:comp_baseline} Comparison of baseline firing properties between cells and their corresponding fits. The histograms on top compare the distributions of the n= 54 cells in blue and their respective models in orange. The scatter plot at the bottom directly compares them. Points on the identity line (grey) indicate perfect model predictions. \textbf{A}: The baseline firing rate of the cell and the model. The base rate agrees near perfectly as it is set to be equal within a margin of 2\,Hz during the fitting process. \textbf{B}: The vector strength agrees well for most cells but for some cells with a VS above 0.8 the models to show a weaker VS than the cell. \textbf{C}: The models fail to show the same strong negative SC at lag 1. this effect gets stronger the more the SC deviates from -0.4.}
\end{figure}
The last two baseline firing behaviors are the burstiness and the coefficient of variation (CV) and both were fit quite similar, because both are correlated as shown below in figure \ref{fig:behavior_correlations} and by the color coding by the computed cell burstiness of the scatter plot (fig. \ref{fig:comp_burstiness}). The model fits lower half of them well but doesn't manage to match the high values, where it consistently underestimates the corresponding value.
The cells' burstiness distribution has two peaks: the continuously firing cells around 0 and the bursting cells around 2. These two main peaks are still fit quite well but the rare very strongly bursting cell is not matched by the model, but this may be an artifact of how burstiness was defined here, as the ISI histograms seem to contain a full continuum between regular firing as seen in fig. \ref{fig:heterogeneity_isi_hist} A and the strongly bursting cells as in C. The covariance is distributed more evenly but still shows the two peaks seen in the burstiness measure.
\begin{figure}[H]
\includegraphics{figures/fit_burstiness_comparison.pdf}
\caption{\label{fig:comp_burstiness} Comparison of the cell behavior and the behavior of the corresponding fit: \textbf{A} burstiness, \textbf{B} coefficient of variation (CV). \textbf{A}: The model values for the burstiness agree well with the values of the model but again show a tendency that the higher the value of the cell the more the model value is below it. \textbf{B}: The CV also shows the problem of the burstiness but the values drift apart more slowly starting around 0.6.}
\caption{\label{fig:comp_burstiness} Comparison of baseline bursting properties between cells and their corresponding fits. The histograms on top compare the distributions of the n= 54 cells in blue and their respective models in orange. The scatter plot at the bottom directly compares them. Points on the identity line (grey) indicate perfect model predictions. \textbf{A}: The model values for the burstiness agree well with the values of the model but again show a tendency that the higher the value of the cell the more the model value is below it. \textbf{B}: The CV also shows the problem of the burstiness but the values drift apart more slowly starting around 0.6.}
\end{figure}
The models were able to fit the steady state response of the cells very well (\ref{fig:comp_adaption} A). There are fluctuations but they are even so there is no strong bias in one direction. Most cells react to a 10\% increase in EOD amplitude with around an increase of their firing frequency of around 40\,Hz.
The fit of the onset response characterized by the slope of the Boltzmann function shows very strong fluctuations which make an accurate judgment difficult, as even with these large differences the quality of the model is still often decent.
\begin{figure}[H]
\includegraphics{figures/fit_adaption_comparison.pdf}
\caption{\label{fig:comp_adaption} Comparison of the cell behavior and the behavior of the corresponding fit: \textbf{A} steady state $f_\infty$ and \textbf{B} onset $f_0$ response slope. In grey the line on which cell and model values are equal. \todo{how many} value pairs from \textbf{B} lie outside of the shown area. They had slopes between \todo{}. In \textbf{A} the $f_\infty$ slope pairs. Cell and models show good agreement with a low scattering in both direction. \textbf{B} The $f_0$ values show a higher spread and for steeper slopes the models have more often too flat slopes.}
\caption{\label{fig:comp_adaption} Comparison of adaption properties between cells and their corresponding fits. The histograms on top compare the distributions of the n=54 cells in blue and their respective models in orange. The scatter plot at the bottom directly compares them. Points on the identity line (grey) indicate perfect model predictions. \textbf{A}: The $f_\infty$ slope pairs show good agreement with mostly low scattering in both direction. \textbf{B} The $f_0$ values show a higher spread and for steeper slopes the models have more often too flat slopes.}
\end{figure}
Given the differences in the between the cell firing properties and the ones of the model the correlations between the different firing properties were calculated and show differences (fig. \ref{fig:behavior_correlations}). Of the seven correlations found in the data set the fitted models show all except the correlation between the VS and baseline firing rate, but the models also show four additional correlations. These are between the baseline firing rate and the $f_0$ slope, base rate and burstiness again base rate and SC and between SC and $f_\infty$ slope.
Before the parameter distributions (fig. \ref{fig:parameter_distributions}) and correlations (fig. \ref{fig:parameter_correlations}) of the model parameters were closer investigated, the potential influence of the different EOD frequencies was removed by scaling the time dependent parameters for all models. This was done by calculating the factor between the fish's EOD frequency and the chosen EOD frequency of 800\,Hz and then multiplying all time parameters appropriately to their dependence with the factor. These scaled parameter distributions are shown in figure \ref{fig:parameter_distributions}. With these scaled distributions the correlations between the parameters were computed giving the matrix in figure \ref{fig:parameter_correlations}. it shows extensive correlations between most parameters. The correlations indicate that the parameter can compensate for each other and that the model can produce similar firing properties for different parameter sets. A notable exception is the refractory period $t_{ref}$ which is independent of all other parameters and could as such the only variable influencing the burstiness in this model.
\begin{figure}[H]
\includegraphics{figures/behaviour_correlations.pdf}
\caption{\label{fig:behavior_correlations} Significant correlations between the behavior variables in the data and the fitted models $p < 0.05$ (Bonferroni corrected). The models contain all the same correlations as the data except for the correlation between the baseline firing rate and the VS, but they also show four additional correlations not seen within the cells: bursting - base rate, SC - $f_\infty$ slope, $f_0$ slope - base rate, SC - base rate.}
\caption{\label{fig:behavior_correlations} Significant correlations between the firing properties in the data and the fitted models (Significance $p < 0.05$ Bonferroni corrected). The models contain all the same correlations as the data except for the correlation between the baseline firing rate and the VS, but they also show four additional correlations not seen within the cells: bursting - base rate, SC - $f_\infty$ slope, $f_0$ slope - base rate, SC - base rate.}
\end{figure}
\begin{figure}[H]
\includegraphics{figures/parameter_distributions.pdf}
\caption{\label{fig:parameter_distributions} Distributions of all eight model parameters. \textbf{A}: input scaling $\alpha$, \textbf{B}: Bias current $I_{Bias}$, \textbf{C}: membrane time constant $\tau_m$, \textbf{D}: noise strength $\sqrt{2D}$, \textbf{E}: adaption time constant $\tau_A$, \textbf{F}: adaption strength $\Delta_A$, \textbf{G}: time constant of the dendritic low pass filter $\tau_{dend}$, \textbf{H}: refractory period $t_{ref}$}
\caption{\label{fig:parameter_distributions} Distributions of all eight model parameters with the time scaled for all models so their driving EOD frequency has 800\,Hz. \textbf{A}: input scaling $\alpha$, \textbf{B}: Bias current $I_{Bias}$, \textbf{C}: membrane time constant $\tau_m$, \textbf{D}: noise strength $\sqrt{2D}$, \textbf{E}: adaption time constant $\tau_A$, \textbf{F}: adaption strength $\Delta_A$, \textbf{G}: time constant of the dendritic low pass filter $\tau_{dend}$, \textbf{H}: refractory period $t_{ref}$}
\end{figure}
\todo{image with rescaled time parameters to 800\,Hz, add to above figure?}
\begin{figure}[H]
\includegraphics{figures/parameter_correlations.pdf}
\caption{\label{fig:parameter_correlations} Significant correlations between model parameters $p < 0.05$ (Bonferroni corrected).}
\caption{\label{fig:parameter_correlations} Correlations between model parameters (Significance $p < 0.05$ Bonferroni corrected). The model parameters show many correlation between each other indicating strong compensatory effects where multiple parameter sets can lead to the same firing properties. The only parameter without correlations to any other is $t_{ref}$ showing a certain independence compared to the other parameters.}
\end{figure}
\subsection{Random Model Population}
Finally the scaled parameter distributions of the fitted models were used to compute a representative multivariante normal space out of which random parameter sets for the models could be drawn.
Each parameter distribution was fit by Gaussian function and then tweaked by hand to remove overly wide Gaussian fits that would produce many parameters outside of the observed distributions (fig. \ref{fig:parameter_dist_with_gauss_fits}). The standard deviation of the fits and the calculated correlations between the parameters were used to compute the covariances. The covariances and the means of the Gaussian fits were then used to compute the multivariante normal space.
\begin{figure}[H]
\includegraphics{figures/parameter_distribution_with_gauss_fits.pdf}
\caption{\label{fig:parameter_dist_with_gauss_fits} Gauss fits used as approximations for the parameter distribution. In black the Gaussian fit used. All parameters except for $t_{ref}$ and $I_{Bias}$ were log transformed to get a more Gaussian distribution. \textbf{A}: Log input scaling $\alpha$, \textbf{B}: bias current $I_{Bias}$, \textbf{C}: Log membrane time constant $\tau_m$, \textbf{D}: Log noise strength $\sqrt{2D}$, \textbf{E}: Log adaption time constant $\tau_A$, \textbf{F}: Log adaption strength $\Delta_A$, \textbf{G}: Log time constant of the dendritic low pass filter $\tau_{dend}$, \textbf{H}: refractory period $t_{ref}$}
\caption{\label{fig:parameter_dist_with_gauss_fits} Gaussian fits (black) used as approximations for the parameter distributions. All parameters except for $t_{ref}$ and $I_{Bias}$ were log transformed to get a more Gaussian-like distribution. \textbf{A}: Logarithmic input scaling $\alpha$, \textbf{B}: bias current $I_{Bias}$, \textbf{C}: Logarithmic membrane time constant $\tau_m$, \textbf{D}: Logarithmic noise strength $\sqrt{2D}$, \textbf{E}: Logarithmic adaption time constant $\tau_A$, \textbf{F}: Logarithmic adaption strength $\Delta_A$, \textbf{G}: Logarithmic time constant of the dendritic low pass filter $\tau_{dend}$, \textbf{H}: refractory period $t_{ref}$}
\end{figure}
\newpage
100 parameter sets were drawn their distribution and correlation compared to the parameters of fitted models to verify the drawing. The distributions shown in figure \ref{fig:drawn_parameter_dist} for most parameters good agreement. the distributions of the input scaling $\alpha$ and the adaption strength $\Delta_A$ are a bit shifted from the reference distribution and the distribution of the bias current $I_{Bias}$ doesn't show the tail towards negative values. The correlations also showed differences even though they were part of the computation of the multivariante normal space. The drawn parameter sets have 3 missing and one additional correlation. The missing correlations are between $\tau_{dend}$ and $\sqrt{2D}$, $\tau_{dend}$ and $\Delta_A$ and $\tau_A$ and $\alpha$, while the correlation between $\alpha$ and $\tau_m$ was only found in the drawn parameter sets. This calculation was repeated a few times during the analysis and these four correlations were inconsistent in the drawn models indicating they may not be well defined in the calculated covariances.
Even with these differences the firing property distributions (fig. \ref{fig:drawn_behavior_dist}) shown by the random models show large overlaps. The random models show in the baseline firing rate distribution, that there are too few models with low firing rates and also has a long tail up too 800\,Hz which does not match the cells. The serial correlation distribution is shifted towards weaker negative correlations. The worst fit distribution is the VS where the random models have a nearly distinct distribution compared to the cells with a lot lower VS.
\begin{figure}[H]
\includegraphics{figures/compare_parameter_dist_random_models.pdf}
\caption{\label{fig:drawn_parameter_dist} Parameter distribution between randomly drawn models \todo{color}orange and the fitted ones blue\todo{color}. \textbf{A}: input scaling $\alpha$, \textbf{B}: Bias current $I_{Bias}$, \textbf{C}: membrane time constant $\tau_m$, \textbf{D}: noise strength $\sqrt{2D}$, \textbf{E}: adaption time constant $\tau_A$, \textbf{F}: adaption strength $\Delta_A$, \textbf{G}: time constant of the dendritic low pass filter $\tau_{dend}$, \textbf{H}: refractory period $t_{ref}$}
\includegraphics{figures/rand_parameter_correlations_comparison.pdf}
\caption{\label{fig:drawn_parameter_corr} Parameter correlation comparison between the fitted parameters and the ones drawn from the multivariante normal distribution. There are four correlations that do not agree between the two, but those are inconsistent in the drawn models (see discussion).}
\end{figure}
\begin{figure}[H]
\includegraphics{figures/rand_parameter_correlations_comparison.pdf}
\caption{\label{fig:drawn_parameter_corr} Parameter correlation comparison between the fitted parameters and the ones drawn from the multivariant normal distribution. There are four correlations that do not agree between the two, but those are inconsistent in the drawn models (see discussion).}
\includegraphics{figures/compare_parameter_dist_random_models.pdf}
\caption{\label{fig:drawn_parameter_dist} Parameter distribution between randomly drawn models (gray) and the fitted ones (orange). \textbf{A}: input scaling $\alpha$, randomly drawn parameters are shifted to higher values. \textbf{B}: Bias current $I_{Bias}$ doesn't match the long tail into the negative. \textbf{C}: membrane time constant $\tau_m$, \textbf{D}: noise strength $\sqrt{2D}$, \textbf{E}: adaption time constant $\tau_A$, \textbf{F}: adaption strength $\Delta_A$, \textbf{G}: time constant of the dendritic low pass filter $\tau_{dend}$, \textbf{H}: refractory period $t_{ref}$}
\end{figure}
\begin{figure}[H]
\includegraphics{figures/random_models_behaviour_dist.pdf}
\caption{\label{fig:drawn_behavior_dist} Behavior distribution of the randomly drawn models \todo{color}(orange) and the original cells \todo{color}(blue). The distribution of the seven behavior characteristics agree well for the most part, but especially the vector strength (VS) in \textbf{G} is offset to the distribution seen in the cells.}
\caption{\label{fig:drawn_behavior_dist} Distribution of firing properties from randomly drawn models (orange) and the original cells (blue). The distribution of the seven firing properties agree well, but especially the vector strength (VS) in \textbf{C} is offset to the distribution seen in the cells and shows manly lower values and the burstiness and SC are also slightly offset.}
\end{figure}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%DISCUSSION
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\newpage
\section{Discussion}
@ -540,9 +570,15 @@ In this thesis a simple model based on the leaky integrate-and-fire (LIF) model
\todo{comparison to the other already existing models}
%Previous studies limited themselves in the P-unit behavior that was modeled.
% \cite{kashimori1996model} introduced a conductances based model but the estimation of variables is difficult because of the little experimental data that is available.
%In \cite{chacron2001simple} two P-units were modeled considering only the baseline behaviour with one bursty cell and one regularly firing one as representatives.
\subsection{Model fit}
The dendtritic low pass filter and the refractory period were necessary for the model to match the firing behavior of the P-units (fig. \ref{fig:dend_ref_effect}). As \cite{walz2013Phd} demonstrated a model without the low pass filter is not able to match the VS and locks too strongly to the EOD. A refractory period $t_{ref}$ is necessary for the model to deviate from the Gaussian firing probability \todo{explain what that is} and show bursting behavior and is flexible enough to match different strengths of burstiness.
With these additions the behavior of the cells was generally matched well by the models with very similar final distributions of the firing properties but there were some limitations. The model failed to reproduce cells with a very high burstiness (long bursts with long pauses between) or a high coefficient of variation could not fully be matched by the model \todo{burst or CV not really correct because of the correlation!}. As seen in figure \ref{fig:comp_burstiness}. The example of fig. \ref{fig:example_bad_isi_fits} \textbf{A} is a case where the model can show this type firing behavior (long bursts and pauses) but it seems difficult to reach the parameter configuration needed with the fitting approach used. In contrast to that the firing behavior of the cells in fig. \ref{fig:example_bad_isi_fits} \textbf{B} and \textbf{C} are not possible for the model in its current form. The addition of the refractory period $t_{ref}$ does not also allow for an increased firing probability at the 2nd EOD period and the cell \textbf{C} shows a higher order structure in its ISI histogram on a comparatively long timescale which the proposed simple model cannot reproduce. These kind of cells showing higher order structure in their ISI histogram are rare but might provide interesting insights in the physiological properties of P-units when further studied.
With these additions the behavior of the cells was generally matched well by the models with very similar final distributions of the firing properties but there were some limitations. The model failed to reproduce cells with a very high burstiness (long bursts with long pauses between) or a high coefficient of variation could not fully be matched by the model \todo{burst or CV not really correct because of the correlation!} (fig.~\ref{fig:comp_burstiness}). The example of fig. \ref{fig:example_bad_isi_fits} \textbf{A} is a case where the model can show this type firing behavior (long bursts and pauses) but it seems difficult to reach the parameter configuration needed with the fitting approach used. In contrast to that the firing behavior of the cells in fig. \ref{fig:example_bad_isi_fits} \textbf{B} and \textbf{C} are not possible for the model in its current form. The addition of the refractory period $t_{ref}$ does not also allow for an increased firing probability at the 2nd EOD period and the cell \textbf{C} shows a higher order structure in its ISI histogram on a comparatively long timescale which the proposed simple model cannot reproduce. These kind of cells showing higher order structure in their ISI histogram are rare but might provide interesting insights in the physiological properties of P-units when further studied.
Two firing properties had a high spread in the fitted models. In the serial correlation the models had some tendency to underestimate the cell's SC. The second property was the slope of the $f_0$ response. Here one possible source is that the fitted Boltzmann function and its slope are quite sensitive to miss-detections of spikes. A wrong estimate of the firing frequency for a single contrast can strongly influence the slope of the fitted Boltzmann function. Unlike the baseline firing properties there don't seem to be cases in which the model cannot fit the f-I curves. The problematic cases shown in figure \ref{fig:example_bad_fi_fits} are both generally possible (fig. \ref{fig:example_good_fi_fits}) so improvements in the cost function and fitting routine should also further improve the model consistency for the adaption responses.
@ -552,54 +588,22 @@ This was also looked at in \cite{walz2013Phd} but only 23 cells were used and th
The parameters of the fitted models also showed extensive correlations between each other. This is an indication of strong compensation effects between them \citep{olypher2007using}. Which is especially clear for the input gain $\alpha$ and the bias current $I_{Bias}$ that have a nearly perfect correlation and control the models baseline firing rate together. Note that the refractory period $t_{ref}$ is the only completely independent variable. This might show a certain independence between the strength of the burstiness and the other firing characteristics, which could be more closely investigated by looking at the sensitivity of models firing properties to changes in $t_{ref}$.
\subsection{Heterogeneous Population}
The correlations and the estimated parameter distributions were used form of their covariances to draw random parameter sets from a multivariante normal distribution. The drawn parameters show the expected distributions but different correlations. That could mean that the \todo{number} models used to calculate them were to few to give enough statistical power for the correct estimation of all correlations. Drawing more models and compensating for the increase in power showed that the involved correlations stay inconsistent, which points to an uncertainty already in the measured covariance matrix of the data. This could be further investigated with a robustness analysis estimating the reliability of the computed covariances.
The correlations and the estimated parameter distributions were used form of their covariances to draw random parameter sets from a multivariante normal distribution. The drawn parameters show the expected distributions but different correlations. That could mean that the 54 models used to calculate them were to few to give enough statistical power for the correct estimation of all correlations. Drawing more models and compensating for the increase in power showed that the involved correlations stay inconsistent, which points to an uncertainty already in the measured covariance matrix of the data. This could be further investigated with a robustness analysis estimating the reliability of the computed covariances.
The firing behavior shown by the drawn models on the other hand fits the ones of the data quite well except for the VS, where it is consistently underestimating the VS of the data.
%Previous studies limited themselves in the P-unit behavior that was modeled.
% \cite{kashimori1996model} introduced a conductances based model but the estimation of variables is difficult because of the little experimental data that is available.
%In \cite{chacron2001simple} two P-units were modeled considering only the baseline behaviour with one bursty cell and one regularly firing one as representatives.
\subsection{Conclusion}
In general the model is the first that takes the burstiness as a continuum into account and seems to be able to accurately describe the firing behavior in a large part of the behavior space of the P-units. But further testing is required to get a clearer picture where and why discrepancies exist. An important next step is the verification of the models with a different type of stimulus. For this a stimulus with random or sinusoidal amplitude modulations could be used. The correlations also need further investigation. As a first step a robustness test could be done to estimate if there are correlations that are not well characterized in both the cells and the models.
\todo{Doesn't cover long timescales described in Gussin 2007??}
\todo{Why do we want such models - analysis of coding in heterogeneic neuron populations, possibility to "measure" responses from whole population at the same time to a single stimulus, separating the different types and analyzing their specific coding properties - finding out why the heterogeneity is necessary!}
\todo{comparison to existing models Chacron, Waltz, Kashimori what does this model add which the others "missed" don't deliver on.}
\newpage
\begin{comment}
\subsection*{Fitting quality}
\begin{itemize}
\item strongly bursty cells not well fitted with a gap between first EOD and the gauss distribution, but this is possible in the model just seems to be a "hard to reach" parameter combination. more start parameters changes in cost function.
\item additionally bursties might need more data to be "well defined" because of the more difficult pre-analysis of the cell itself (high variance in rate traces even with a mean of 7-10 trials)
\item different fitting routine / weights might also improve consistency of the fitting
\item Model CAN fit most types of cells except for double burst spikes in ISI hist (1st and 2nd EOD have high probability) and some with higher level structure in ISI hist are probably not possible with the current model
\item f-I curve works but is still not fully consistent slope of $f_0$ strongly affected by miss detections so not the best way of validating the $f_0$ response.
\item b-correlations:
\item Data correlation between base rate and VS unexpected and missing in the model.
\item data correlation between base rate and $f_0$ expected as $f_0$ ~ $f_\infty$ and $f_\infty$ ~ base rate (appeared in the model)
\item Model correlation between base rate and burstiness probably error because of the problems to "fit the gap" between burst and other ISIs. (appeared in model)
\item above probably causes base rate - SC correlation and this causes the additional SC - $f_\infty$ correlation over $f_\infty$ ~ base rate
\item In total it can deliver good models over a large space of the heterogeneity especially with the addition burstiness, but it is not yet verified with a different stimulus type like RAM or SAM and also need further investigation to it's quality for example because of the mismatched correlations $\rightarrow$ robustness analysis should be done
\end{itemize}
\subsection*{random models}
\begin{itemize}
\item Gauss fits are in some cases questionable at best
\item resulting parameter distributions as such also not really that good
\item Parameter correlations have 4--5 correlations whose significance is strongly inconsistent/random even when using 1000 drawn models (while compensating for higher power): thus acceptable result??
\item behavior distribution not perfect by any means but quite alright except for the VS. Which definitely needs improvement! Maybe possible with more tweaking of the gauss fits.
\end{itemize}
\end{comment}
\newpage
\bibliography{citations}
\bibliographystyle{apalike}

35
thesis/citations.bib
View File

@ -306,6 +306,41 @@
publisher={National Acad Sciences}
}
@article{chacron2004burst,
title={To burst or not to burst?},
author={Chacron, Maurice J and Longtin, Andr{\'e} and Maler, Leonard},
journal={Journal of computational neuroscience},
volume={17},
number={2},
pages={127--136},
year={2004},
publisher={Springer}
}
@article{nelder1965simplex,
title={A simplex method for function minimization},
author={Nelder, John A and Mead, Roger},
journal={The computer journal},
volume={7},
number={4},
pages={308--313},
year={1965},
publisher={Oxford University Press}
}
@article{zeldenrust2018neural,
title={Neural coding with bursts—current state and future perspectives},
author={Zeldenrust, Fleur and Wadman, Wytse J and Englitz, Bernhard},
journal={Frontiers in computational neuroscience},
volume={12},
pages={48},
year={2018},
publisher={Frontiers}
}

BIN
thesis/figures/amGeneration.pdf
View File

Binary file not shown.

BIN
thesis/figures/behaviour_correlations.pdf
View File

Binary file not shown.

BIN
thesis/figures/compare_parameter_dist_random_models.pdf
View File

Binary file not shown.

BIN
thesis/figures/dend_ref_effect.pdf
View File

Binary file not shown.

BIN
thesis/figures/example_bad_fi_fits.pdf
View File

Binary file not shown.

BIN
thesis/figures/example_bad_isi_hist_fits.pdf
View File

Binary file not shown.

BIN
thesis/figures/example_good_fi_fits.pdf
View File

Binary file not shown.

BIN
thesis/figures/example_good_isi_hist_fits.pdf
View File

Binary file not shown.

BIN
thesis/figures/f_point_detection.pdf
View File

Binary file not shown.

BIN
thesis/figures/fit_adaption_comparison.pdf
View File

Binary file not shown.

BIN
thesis/figures/fit_baseline_comparison.pdf
View File

Binary file not shown.

BIN
thesis/figures/fit_burstiness_comparison.pdf
View File

Binary file not shown.

BIN
thesis/figures/isi_hist_heterogeneity.pdf
View File

Binary file not shown.

BIN
thesis/figures/model_comparison.pdf
View File

Binary file not shown.

BIN
thesis/figures/p_unit_example.pdf
View File

Binary file not shown.

BIN
thesis/figures/parameter_correlations.pdf
View File

Binary file not shown.

BIN
thesis/figures/parameter_distribution_with_gauss_fits.pdf
View File

Binary file not shown.

BIN
thesis/figures/parameter_distributions.pdf
View File

Binary file not shown.

BIN
thesis/figures/rand_parameter_correlations_comparison.pdf
View File

Binary file not shown.

BIN
thesis/figures/random_models_behaviour_dist.pdf
View File

Binary file not shown.

BIN
thesis/figures/scaled_to_800_parameter_distributions.pdf
View File

Binary file not shown.

BIN
thesis/figures/stimulus_development.pdf
View File

Binary file not shown.

BIN
thesis/sources/notes/notes.pdf
View File

Binary file not shown.

32
thesis/sources/notes/notes.tex
View File

@ -38,6 +38,8 @@ Bastian 1981a Electrolocation I. How the electroreceptors of Apteronotus albifro
\section{Introduction}
\begin{enumerate}
\item sensory input important for all life etc.
\item electric fish
\begin{enumerate}
\item general: habitat,
@ -59,6 +61,7 @@ Bastian 1981a Electrolocation I. How the electroreceptors of Apteronotus albifro
\end{enumerate}
\subsection{Apteronotus leptorynchus}
\begin{itemize}
@ -227,11 +230,38 @@ FI-Curve:
\section{Discussion}
\subsection{Fitting quality}
\begin{itemize}
\item todo
\item strongly bursty cells not well fitted with a gap between first EOD and the gauss distribution, but this is possible in the model just seems to be a "hard to reach" parameter combination. more start parameters changes in cost function.
\item additionally bursties might need more data to be "well defined" because of the more difficult pre-analysis of the cell itself (high variance in rate traces even with a mean of 7-10 trials)
\item different fitting routine / weights might also improve consistency of the fitting
\item Model CAN fit most types of cells except for double burst spikes in ISI hist (1st and 2nd EOD have high probability) and some with higher level structure in ISI hist are probably not possible with the current model
\item f-I curve works but is still not fully consistent slope of $f_0$ strongly affected by miss detections so not the best way of validating the $f_0$ response.
\item b-correlations:
\item Data correlation between base rate and VS unexpected and missing in the model.
\item data correlation between base rate and $f_0$ expected as $f_0$ ~ $f_\infty$ and $f_\infty$ ~ base rate (appeared in the model)
\item Model correlation between base rate and burstiness probably error because of the problems to "fit the gap" between burst and other ISIs. (appeared in model)
\item above probably causes base rate - SC correlation and this causes the additional SC - $f_\infty$ correlation over $f_\infty$ ~ base rate
\item In total it can deliver good models over a large space of the heterogeneity especially with the addition burstiness, but it is not yet verified with a different stimulus type like RAM or SAM and also need further investigation to it's quality for example because of the mismatched correlations $\rightarrow$ robustness analysis should be done
\end{itemize}
\subsection{random models}
\begin{itemize}
\item Gauss fits are in some cases questionable at best
\item resulting parameter distributions as such also not really that good
\item Parameter correlations have 4--5 correlations whose significance is strongly inconsistent/random even when using 1000 drawn models (while compensating for higher power): thus acceptable result??
\item behavior distribution not perfect by any means but quite alright except for the VS. Which definitely needs improvement! Maybe possible with more tweaking of the gauss fits.
\end{itemize}