\documentclass[12pt,a4paper,pdftex]{article}
\usepackage[left=25mm, right=25mm, top=20mm, bottom=25mm]{geometry}
\usepackage{graphicx}
\usepackage{amsmath}
\usepackage{natbib}
\usepackage{xcolor}
\usepackage[breaklinks=true,colorlinks=true,citecolor=blue!30!black,urlcolor=blue!30!black,linkcolor=blue!30!black]{hyperref}


\usepackage[utf8x]{inputenc}
\usepackage[english]{babel}
%\usepackage{float}
\usepackage{floatrow}
\usepackage{wrapfig}

\usepackage{listings} % für den code am Ende



\newcommand{\todo}[1]{{\color{red}(TODO: #1)}}

\newcommand{\AptLepto}{{\textit{Apteronotus leptorhynchus\:}}}

\newcommand{\lepto}{{\textit{A. leptorhynchus\:}}}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Ab hier beginnt der eigentliche Text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{document}


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Titelseite
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{titlepage}
		\begin{center}
		{\Huge Modeling the Heterogeneity of Electrosensory Afferents in Electric Fish \par}
		\vspace{0.75cm}
		{\Large Masterthesis \par}
		\vspace{0.25cm}
		{der Mathematisch-Naturwissenschaftlichen Fakultät \par} {der Eberhard Karls Universität Tübingen \par}
		\vspace{0.75cm}
		{Erstkorrektor: Prof.~Dr.~Philipp Berens\\
		Zweitkorrektor: Prof.~Dr.~Jan Benda \par}
		\vspace{0.25cm}
		{Lehrbereich für Neuroethologie}
		
		\vfill
		\large vorgelegt von \par
		\large Alexander Mathias Ott \par
		Abgabedatum: 21.09.2020
	\end{center}
\end{titlepage}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Erklärung
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section*{Eigenständigkeitserklärung}

\vspace{0.5cm}

Hiermit erkläre ich, dass ich die vorgelegte Arbeit selbstständig verfasst habe und keine anderen als die angegebenen Quellen und Hilfsmittel benutzt habe. 

\vspace{2mm}
\noindent
Außerdem erkläre ich, dass die eingereichte Arbeit weder vollständig noch in wesentlichen Teilen Gegenstand eines anderen Prüfungsverfahrens gewesen ist.






\vfill

\begin{tabular}{ll}
	$\overline{\text{Unterschrift}\hspace{6cm}}$ & $\overline{\text{Ort, Datum}\hspace{4cm}}$ \\
\end{tabular}
\newpage\newpage



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Inhalsverzeichnis
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

{
 \hypersetup{linkcolor=black}
 \tableofcontents
}

\newpage

\section*{Not to forget: TODO}

\begin{itemize}

\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

\end{itemize}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Zusammenfassung
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Zusammenfassung}
% Abstract in deutsch




\section{Abstract}
%Einleitung + Ergebnisse der Diskussion in kurz


\newpage
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Einleitung
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Introduction}
\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

The environment of an organism holds important information that it needs to survive. Information about predators to avoid, food to find and potential mates. That means that 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).


%

\begin{figure}[H]
{\caption{\label{fig:p_unit_example} Example behavior of a p-unit with a high baseline firing rate and an EODf of 744\,Hz. \textbf{A}: A 100\,ms voltage trace of the baseline recording with spikes marked by the black lines. \textbf{B}: The histogram of the ISI with the x-axis in EOD periods, showing the phase locking of the firing. \textbf{C}: The serial correlation of the ISI showing a negative correlation for lags one and two. \textbf{D}: The response of the p-unit to a step increase in EOD amplitude. In \todo{color} the averaged frequency over 10 trials. The P-unit strongly reacts to the onset of the stimulus but very quickly adapts to the new stimulus and then shows a steady state response. \textbf{E}: The f-I curve visualizes the onset and steady-state response of the neuron for different step sizes (contrasts). In \todo{color} the detected onset responses and the fitted Boltzmann, in \todo{color} the detected steady-state response and the linear fit.}}
{\includegraphics[width=1\textwidth]{figures/p_unit_example.pdf}}
\end{figure}
\newpage



\begin{figure}[H]
{\caption{\label{fig:heterogeneity_isi_hist} Variability in spiking behavior between P units under baseline conditions. \textbf{A--C} 100\,ms of cell membrane voltage and \textbf{D--F} interspike interval histograms, each for three different cells. \textbf{A} and \textbf{D}: A non bursting cell with a baseline firing rate of 133\,Hz (EODf: 806\,Hz), \textbf{B} and \textbf{E}: A cell with some bursts and a baseline firing rate of 235\,Hz (EODf: 682\,Hz) and \textbf{C} and \textbf{F}: A strongly bursting cell with longer breaks between bursts. Baseline rate of  153\,Hz and EODf of 670\,Hz }}
{\includegraphics[width=\textwidth]{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. 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}. Even though P units were already modelled 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}). There is up to this point 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, 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.






%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Methoden
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Materials and Methods}

\todo{some transition from the introduction}

\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 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 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)
For the recordings a standard glass mircoelectrode (borosilicate; 1.5 mm outer diameter; GB150F-8P, Science Products, Hofheim, Germany) was used. They were pulled to a resistance of 50--100\,M$\Omega$ using Model P-97 from Sutter Instrument Co. (Novato, CA, USA) and filled with 1\,M KCl solution. The electrodes were controlled using microdrives (Luigs-Neumann; Ratingen, Germany) and the potentials recorded with the bridge mode of the SEC-05 amplifier (npi-electronics GmbH, Tamm, Germany) and lowpass filtered at 10 kHz.

During the recording spikes were detected online using the peak detection algorithm from \cite{todd1999identification}. It uses a dynamically adjusted threshold value above the previously detected trough. To detect spikes through changes in amplitude the threshold was set to 50\% of the amplitude of a detected spike while keeping the threshold above a minimum set to be higher than the noise level based on a histogram of all peak amplitudes. Trials with bad spike detection were removed from further analysis.
The fish's EOD was recorded using two vertical carbon rods (11\,cm long, 8\,mm diameter) positioned in front of the head and behind its tail. The signal was amplified 200 to 500 times and band-pass filtered (3 − 1500 Hz passband, DPA2-FX, npi-electronics, Tamm, Germany). The electrodes were placed on iso-potential lines of the stimulus field to reduce the interference of the stimulus in the recording. All signals were digitized using a data acquisition board (PCI-6229; National Instruments, Austin TX, USA) at a sampling rate of 20--100\,kHz (54 cells at 20\,kHz, 20 at 100\,kHz and 1 at 40\,kHz)

The recording and stimulation was done using the ephys, efield, and efish plugins of the software RELACS (\href{www.relacs.net}{www.relacs.net}). It allowed the online spike and EOD detection, pre-analysis and visualization and ran on a Debian computer.



\subsection{Stimulus Protocols}
% image of Baseline stimulus as baseline doesn't mean no stimulus here
% image of Fi curve stimulus sinusoidal step
% image of SAM stimulus

The stimuli used during the recordings were presented from two vertical carbon rods (30 cm long, 8 mm diameter) as stimulus electrodes. They were positioned at either side of the fish parallel to its longitudinal axis. The stimuli were computer generated, attenuated and isolated  (Attenuator: ATN-01M, Isolator: ISO-02V, npi-electronics, Tamm, Germany) and then send to the stimulus electrodes.
For this work two types of recordings were made with all cells: baseline recordings and amplitude step recordings for the frequency-Intensity curve (f-I curve). 
The 'stimulus' for the baseline recording is purely the EOD field the fish produces itself with no external stimulus.

The amplitude step stimulus here is a step in EOD amplitude. The amplitude modulation (AM) is measured as a contrast. The contrast is calculated by dividing the EOD amplitude during the step by the normal EOD amplitude. To be able to cause a given AM in the fish's EOD, the EOD was recorded and multiplied with the modulation (see fig. \ref{fig:stim_examples}). This modified EOD can then be presented at the right phase with the stimulus electrodes, causing constructive interference and adding the used amplitude modulation to the EOD (Fig. \ref{fig:stim_examples}). This stimuli construction as seen in equation~\ref{eq:am_generation} works for any AM as long as the EOD of the fish is stable. 

\begin{equation}
V_{Stim}(t) = EOD(t)(1 + AM(t)) 
\label{eq:am_generation}
\end{equation}


\begin{figure}[H]
\floatbox[{\capbeside\thisfloatsetup{capbesideposition={left, center}, capbesidewidth=0.45\textwidth}}]{figure}[\FBwidth]
{\caption{\label{fig:stim_examples} Example of the stimulus construction. At the top a recording of the fish's EOD. In the middle: EOD recording multiplied with the AM, with a step between 0 and 50\,ms to a contrast of 30\,\% (marked in \todo{color}). At the bottom the resulting stimulus trace when the AM is added to the EOD. \todo{Umformulieren add figure labels A, B, C}}}
{\includegraphics[width=0.45\textwidth]{figures/amGeneration.pdf}}
\end{figure}


All step stimuli consisted of a delay of 0.2\,s followed by a 0.4\,s (n=68) or 1\,s (n=7) long step and a 0.8\,s long recovery time. The contrast range measured  was for the most cells 80--120\% of EOD amplitude. Some cells were measured in a larger range up to 20--180\%. In the range at least 7 contrasts were measured with at least 7 trials, but again many cells were measured with more contrasts and trials. The additionally measured contrasts were used for the model if they had at least 3 trials.


\subsection{Cell Characteristics}

The cells were characterized by ten parameters: 6 for the baseline and 4 for the f-I curve.
For the baseline the mean firing rate was calculated by dividing the number of spikes in the recording by the recording time. Then the set of all interspike intervals (ISI) $T$ was computed and further parameters were calculated from it. 

The coefficient of variation
\begin{equation}
CV = \frac{STD(T)}{\langle T \rangle}
\label{eq:CV}
\end{equation}
is defined as the standard deviation (STD) of $T$ divided by the mean ISI, see equation \ref{eq:CV} with angled brackets as the averaging operator. 


The vector strength (VS) is a measure of how strong the cell locks to a phase of the EOD. It was calculated as seen in Eq. \ref{eq:VS}, by placing each spike on a unit circle depending on the relative spike time $t_i$ of how much time has passed since the start of the current EOD period in relation to the EOD period length. This set of vectors is then averaged and the absolute value of this average vector describes the VS. If the VS is zero the spikes happen equally in all phases of the EOD while if it is one all spikes happen at the exact same phase of the EOD. 

\begin{equation}
vs = |\frac{1}{n} \sum_n e^{iwt_i}|
\label{eq:VS}
\end{equation}


The serial correlation with lag k ($SC_k$) of $T$ is a measure how the ISI $T_i$ (the $i$-th ISI) influences the $T_{i+k}$ the ISI with a lag of k intervals. This is calculated as,  
 
\begin{equation}
SC_k = \frac{\langle (T_{i} - \langle T \rangle)(T_{i+k} - \langle T \rangle) \rangle}{\sqrt{\langle (T_i - \langle T \rangle)^2 \rangle}\sqrt{\langle (T_{i+k} - \langle T \rangle)^2 \rangle}}
\label{eq:SC}
\end{equation}
with the angled brackets again the averaging operator.


Finally the ISI-histogram was calculated within a range of 0--50\,ms and a bin size of 0.1\,ms. The burstiness was calculated as the percentage of ISI smaller than 2.5 EOD periods multiplied by the average ISI. This gives a rough measure of how how often a cell fires in the immediately following EOD periods compared to its average firing frequency.   


%burstiness: \todo{how to write as equation, ignore and don't show an equation?}
% \begin{equation}
% 	b = (T < 1/2EODf)/ N * \langle T \rangle
% \end{equation}

\begin{figure}[H]
\centering
% trim={<left> <lower> <right> <upper>}
%\parbox[c][0mm][t]{80mm}{\hspace{-10.5mm}\large\sffamily A\hspace{50.5mm} \large\sffamily B} 
%\raisebox{70mm}[10]{\large\sffamily A)}
\includegraphics[trim={10mm 5mm 10mm 5mm}, scale=0.8]{figures/f_point_detection.pdf}

\caption{\label{fig:f_point_detection} \textbf{A}: The averaged response of a cell to a step in EOD amplitude. The step of the stimulus is marked by the back bar. The detected values for the onset ($f_0$) and steady-state ($f_{\infty}$) response are marked in \todo{color}. $f_0$ is detected as the highest deviation from the mean frequency before the stimulus while $f_{\infty}$ is the average frequency in the 0.1\,s time window, 25\,ms before the end of the stimulus. \textbf{B}: The fi-curve visualizes the onset and steady-state response of the neuron for different stimuli contrasts. In \todo{color} the detected onset responses and the fitted Boltzmann, in \todo{color} the detected steady-state response and the linear fit.}
\end{figure}


As already mentioned in the introduction, p-units react to a step in EOD amplitude with a strong onset response decaying back to a steady state response (fig.~\ref{fig:f_point_detection}~A). This adaption behavior of the cell was characterized by the f-I curve measurements. First the ISI frequency trace for each stimulus was calculated. The ISI frequency of a time point t is defined as $1/T_i$ with $T_i$ the ISI the time point t falls into. This gives a frequency trace starting with the first spike and ending at the last spike. For further analysis all trials of a specific contrast were averaged over the trials with the resolution of the sampling rate. This results in a trial-averaged step response for each contrast as illustrated in figure \ref{fig:f_point_detection}~A. 
In this firing frequency trace the baseline frequency, the onset $f_0$ and steady-state $f_{\infty}$ response were detected. The baseline frequency was measured as the mean of the firing frequency 25\,ms after recording start up to 25\,ms before the stimulus start. $f_0$ was then defined as the largest deviation from the baseline frequency, within the first 25\,ms after stimulus onset. If there was no deviation farther than the minimum or maximum before the stimulus start, then the average frequency in that 25\,ms time window was used. This approximation made the detection of $f_0$ more stable for small contrasts and trials with high variation.   
The $f_{\infty}$ response was estimated as the average firing frequency in the 100\,ms time window ending 25\,ms before the end of the stimulus (fig. \ref{fig:f_point_detection} A). 
Afterwards a Boltzmann:
\begin{equation}
f_{0}(I) = (f_{max}-f_{min}) (1 / (1 + e^{-k * (I - I_0)})) + f_{min}
\end{equation}
was fitted to the onset response and a rectified line:
\begin{equation}
f_{\infty}(I) = \lfloor mI+c \rfloor_0
\end{equation} 
(with $\lfloor x \rfloor_0$ the rectify operator) was fitted to the steady-state responses (fig.~\ref{fig:f_point_detection}~B).




\subsection{Leaky Integrate and Fire Model}

% add info about simulation by euler integration and which time steps!
% show voltage dynamics with resistance : 

The above described cell characteristics need to be reproduced by a simple and efficient model to be able to simulate bigger populations in a reasonable time. The model used in this thesis follows these equations:

\begin{align}
	\tau_m \frac{dV}{dt} &= -V+I_{Bias} +\alpha V_{dend} - I_{A} + \sqrt{2D}\frac{\xi}{\sqrt{\Delta t}} \label{eq:full_model_dynamics_voltage} \\	
	\tau_A \frac{dI_A}{dt} &= -I_A + \Delta_A \sum \delta (t) \label{eq:full_model_dynamics_adaption} \\
	\tau_{dend} \frac{dV_{dend}}{dt} &= -V_{dend} + \lfloor V_{stim} \rfloor_0 \label{eq:full_model_dynamics_dendrite} 
\end{align}

Equation \ref{eq:full_model_dynamics_voltage} describes the leaky dynamics of the membrane voltage with $\tau_m$ the membrane time constant, $I_{Bias}$ a bias current, $\alpha$ the cell specific gain factor for $V_{dend}$ the input voltage coming from the dendrite. $\sqrt{2D}$ is the strength of the normal distributed noise $\xi$. $I_A$ is an adaption current with the dynamics of equation~\ref{eq:full_model_dynamics_adaption}. $\tau_A$ is the time constant of the adaption, $\Delta_A$ its strength and $\delta (t)$ is the spike train of the cell. Equation~\ref{eq:full_model_dynamics_dendrite} shows the dynamics of the synapse and dendrite with $\tau_{dend}$ the time constant of the dendrite and $\lfloor V_{stim} \rfloor_0$ the rectified stimulus given. Finally the model also includes a refractory period $t_{ref}$, not shown in above equations, that keeps the membrane voltage $V$ at zero for its duration.

To arrive at this model the simplest commonly used neuron model the perfect integrate-and-fire (PIF) model was stepwise extended. The PIF's voltage can be described in one equation: $\tau_m \frac{dV}{dt} = \frac{I}{R_m}$ with $I$ the stimulus current, $R_m$ the membrane resistance and a voltage threshold $V_\theta$. In this model $I$ is integrated and when this threshold $\theta$ is reached the voltage is reset to zero and a spike is recorded (see fig. \ref{fig:model_comparison} PIF). The model is useful for basic simulations but cannot reproduce the richer behavior of the p-units, as it has no memory of previous spikes so it cannot show any adaption behavior and it is also very strongly locked to its limit cycle producing very constant ISI, not allowing the firing flexibility of the p-units. 

The next slightly more complex model is the leaky integrate-and-fire (LIF) model:
\begin{equation}
	\tau_m \frac{dV}{dt} = -V + IR_m
	\label{eq:basic_voltage_dynamics}
\end{equation}
As the name suggests it adds a leakage current to the PIF (fig.~\ref{fig:model_comparison} LIF). The leakage current adds sub threshold behavior to the model and allows for some more flexibility in suprathresold firing but it is still not flexible enough and cannot reproduce the adaption.

To reproduce the adaption behavior the model needs some form of memory of previous spikes. There are two main ways this can be added to the model as an adaptive current or a dynamic threshold. The biophysical mechanism of the adaption in p-units is unknown because the cell bodies are not accessible for intra-cellular recordings. Following the results of \cite{benda2010linear} a negative adaptive current was chosen, because the dynamic threshold causes divisive adaption instead of the subtractive adaption of p-units seen in \cite{benda2005spike}. This results in an leaky integrate-and-fire model with adaption current (LIFAC) (fig.~\ref{fig:model_comparison} LIFAC). The added adaptive current follow the dynamics:

\begin{equation}
	\tau_A \frac{dI_A}{dt} = -I_A + \Delta_A \sum \delta (t)
	\label{eq:adaption_dynamics}
\end{equation}

and gets subtracted from the input current $I$ of of the voltage dynamics eq.~\ref{eq:basic_voltage_dynamics}. It is modeled as an exponential decay with the time constant $\tau_A$ and an adaption strength $\Delta_A$. $\Delta_A$ is multiplied with the sum of spikes $t_i$ in the spike train ($\delta (t_i)$) of the model cell. For the simulation using the Euler integration this results in an increase of $I_A$ by $\frac{\Delta_A}{\tau_A}$ at every time step where a spike is recorded. \todo{image of model simulation with voltage adaption and spikes using the toy model?} The input current $I$ from equation \ref{eq:basic_voltage_dynamics} is a sum of those two currents and an additional bias current $I_{Bias}$ that is needed to adjusts the cells spontaneous spiking:

\begin{equation}
	I = \alpha I_{Input} - I_A + I_{Bias}
	\label{eq:currents_lifac}
\end{equation}

Note that in this p-unit model all currents are measured in mV because, as mentioned above, the cell body is not accessible for intra-cellular recordings and as such the membrane resistance $R_m$ is unknown \todo{ref mem res p-units}. The input current $I_{Input}$ is the current of the stimulus, an amplitude modulated sine wave mimicking the frequency EOD. This stimulus is then rectified to model the receptor synapse and low-pass filtered with a time constant of $\tau_{dend}$ to simulate the low-pass filter properties of the dendrite (fig. \ref{fig:stim_development}) according to:
\begin{equation}
\tau_{dend} \frac{dV_{dend}}{dt} = -V_{dend} + \lfloor I_{Input} \rfloor_0
\end{equation}

Afterwards it is multiplied with $\alpha$ a cell specific gain factor. This gain factor has the unit of cm because the $I_{Input}$ stimulus represents the EOD with a unit of mV/cm.

Finally, noise and an absolute refractory period were added to the model. The noise $\xi$ is drawn from a Gaussian noise distribution and divided by $\sqrt{\Delta t}$ to get a noise which autocorrelation function is independent of the simulation step size $\Delta t$.  The implemented form of the absolute refractory period $t_{ref}$ keeps the model voltage at zero for the duration of $t_{ref}$ after a spike. This gives us the full model described in equations \ref{eq:full_model_dynamics_voltage}--\ref{eq:full_model_dynamics_dendrite}.


\begin{figure}[H]
\includegraphics[scale=0.6]{figures/model_comparison.pdf}
\caption{\label{fig:model_comparison} Comparison of different simple models normed to a spontaneous firing rate of ~10 Hz stimulated with a step stimulus. In the left column y-axis in mV in the right column the y-axis shows the frequency in Hz.  PIF: Shows a continuously increasing membrane voltage with a fixed slope and as such constant frequency for a given stimulus strength. LIF: Approaches a stimulus dependent membrane voltage steady state exponentially Also has constant frequency for a fixed stimulus value. LIFAC: Exponentially approaches its new membrane voltage value but also shows adaption after changes in the stimulus the frequency takes some time to adapt and arrive at the new stable value. }
% LIFAC + ref: Very similar to LIFAC the added absolute refractory period keeps the voltage constant for a short time after the spike and limits high fire rates. \todo{how to deal with the parameters}

\end{figure}

Together this results in the dynamics seen in equations \ref{eq:full_model_dynamics_voltage}--\ref{eq:full_model_dynamics_dendrite}.



\begin{figure}[H]
\includegraphics[scale=0.6]{figures/stimulus_development.pdf}
\caption{\label{fig:stim_development} The stimulus modification in the model. The fish's EOD is simulated with a sin wave. It is rectified at the synapse and then low-pass filtered in the dendrite.}
\end{figure}


\begin{table}[H]
\begin{tabular}{c|l|c}

 parameter & explanation & unit \\
 \hline
  $\alpha$   	&	stimulus scaling factor & [cm] \\
 $\tau_m$    	&	membrane time constant & [ms]\\
 $I_{Bias}$ 	& 	bias current & [mV] \\
 $\sqrt{2D}$ 	&	noise strength & [mV$\sqrt{\text{s}}$]\\
 $\tau_A$ 		&	adaption time constant & [ms] \\
 $\Delta_A$ 	&	adaption strength & [mVms]\\
 $\tau_{dend}$ 	& 	time constant of dendritic low-pass filter & [ms] \\
 $t_{ref}$ 		&	absolute refractory period & [ms]
\end{tabular}
\caption{\label{tab:parameter_explanation} Overview about all parameters of the model that are fitted.}
\end{table}





\subsection{Fitting of the Model}
%überleitung!
The full model has, as described above, eight parameters that need to be fitted so it can reproduce the behavior of the cell. During the fitting and the analysis all models were integrated with at time step of 0.05\,ms.
The stimuli described in the stimulus protocols section above were recreated for the stimulation of the model during the fitting process. The pure fish EOD was approximated by a simple sine wave of the appropriate frequency, but it was decided to keep the amplitude of the sine wave at one to make the models more comparable. Changes in the amplitude can be compensated for by changing the input scaling factor so there is no qualitative difference.

During the fitting the baseline stimulus was simulated 3 times with 30\,s each and the step stimuli were simulated with a delay, step duration and recovery time of each 0.5\,s. The contrasts were the same as in the cell recordings. The step stimuli for the different contrasts were each repeated 8 times. The simulated data was analyzed in the same way as the cells (see above).

The error function was constructed from both the baseline characteristics: VS, CV, SC, ISI-histogram and burstiness and the f-I curve: the detections of $f_{inf}$ and $f_0$ responses for each contrast, the slope of the linear fit into the $f_{inf}$ and the frequency trace of one step response. 

The error of the VS, CV, SC, and burstiness was calculated as the scaled absolute difference:

\begin{equation}
err_i = c_i |x^M_i - x^C_i|
\end{equation}
with $x^M_i$ the model value for the characteristic $i$, $x^C_i$ the corresponding cell value and $c_i$ a scaling factor that is the same for all cells but different between characteristics. The scaling factor was used to make all errors a similar size. They are listed in table \ref{tab:scaling_factors}.

The error for the slope of the $f_{inf}$ fit was the scaled relative difference:

\begin{equation}
err_i = c_i|1 - ((x^M_i - x^C_i) / x^C_i)| 
\end{equation}

For the $f_{inf}$ and $f_0$ responses the average scaled difference off all contrasts was taken and finally the error for the ISI-histogram and the step-response was calculated with a mean-square error. For the histogram over all bins but for the step response only the first 50\,ms after stimulus onset as an error for the adaption time constant.  

\begin{equation}
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.    


\begin{table}[H]
\begin{tabular}{c|c}

 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  \\
 $f_0$ step response		& 0.001
\end{tabular}
\caption{\label{tab:scaling_factors} Scaling factors for fitting errors.}
\end{table}


\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?}

\section{Results}


\begin{figure}[H]
\centering
\includegraphics[width=\textwidth]{figures/dend_ref_effect.pdf}
\caption{\label{fig:dend_ref_effect} 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.}
\end{figure}


\begin{figure}[H]
\includegraphics[scale=0.6]{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. \textbf{A--C} 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.}
\end{figure}


\begin{figure}[H]
\includegraphics[scale=0.6]{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} }
\end{figure}

\begin{figure}[H]
\includegraphics[scale=0.6]{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.}
\end{figure}

\begin{figure}[H]
\includegraphics[scale=0.5]{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.}
\end{figure}

\begin{figure}[H]
\includegraphics[scale=0.5]{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. Excluded \todo{how many} value pairs from $f_0$ slope as they had slopes higher than \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.}
\end{figure}

\begin{figure}[H]
\includegraphics[scale=0.5]{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.}
\end{figure}


\begin{figure}[H]
\includegraphics[scale=0.6]{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.}
\end{figure}

\begin{figure}[H]
\includegraphics[scale=0.6]{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}$}
\end{figure}

\todo{image with rescaled time parameters to 800\,Hz, add to above figure?}

\begin{figure}[H]
\includegraphics[scale=0.6]{figures/parameter_correlations.pdf}
\caption{\label{fig:parameter_correlations} Significant correlations between model parameters $p < 0.05$ (Bonferroni corrected).}
\end{figure}


\begin{figure}[H]
\includegraphics[scale=0.6]{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}$}
\end{figure}


\begin{figure}[H]
\includegraphics[scale=0.6]{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}$}
\end{figure}

\begin{figure}[H]
\includegraphics[scale=0.6]{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).}
\end{figure}

\begin{figure}[H]
\includegraphics[scale=0.6]{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.}
\end{figure}


\section{Discussion}


\subsection*{Fitting quality}

\begin{itemize}
\item many bursty cells not well fitted 
\item bursties might need more data to be well defined (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 but some with weird structure in ISI hist are not possible
\end{itemize}

\subsection*{Behavior correlations}

\begin{itemize}
\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

\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}

\newpage
\bibliography{citations} 
\bibliographystyle{apalike}


\end{document}