correct mat&met 1

Figures_Baseline.py

@@ -16,6 +16,7 @@ DATA_COLOR = "blue"
 
 DATA_SAVE_PATH = "data/figure_data/"
 
+EXAMPLE_CELL = "data/final/2012-12-20-ac-invivo-1"
 
 def main():
     # data_isi_histogram()
@@ -28,8 +29,7 @@ def main():
 
 
 def p_unit_example():
-    # cell = "data/final/2013-04-17-ac-invivo-1"
-    cell = "data/final/2012-04-20-af-invivo-1"
+    cell = EXAMPLE_CELL
     cell_data = CellData(cell)
 
     base = BaselineCellData(cell_data)
@@ -39,11 +39,9 @@ def p_unit_example():
     step = cell_data.get_sampling_interval()
 
     # Overview figure for p-unit behaviour
-
     fig = plt.figure(tight_layout=True, figsize=consts.FIG_SIZE_MEDIUM)
     gs = gridspec.GridSpec(3, 2)
 
-
     # a bit of trace with detected spikes
     ax = fig.add_subplot(gs[0, :])
 
@@ -71,7 +69,7 @@ def p_unit_example():
     ax.hist(isi, bins=bins)
     ax.set_ylabel("Count")
     ax.set_xlabel("ISI in EOD periods")
-    ax.set_title("ISI-histogram")
+    ax.set_title("ISI Histogram")
 
 
     # Serial correlation
@@ -80,8 +78,7 @@ def p_unit_example():
     sc = base.get_serial_correlation(10)
     ax.plot(range(11), [0 for _ in range(11)], color="darkgrey", alpha=0.8)
     ax.plot(range(11), [1] + list(sc))
-
-
+    ax.plot(range(11), [1] + list(sc), '+')
 
     ax.set_xlabel("Lag")
     ax.set_ylabel("SC")
@@ -97,13 +94,13 @@ def p_unit_example():
 
     f_trace_times, f_traces = fi.get_mean_time_and_freq_traces()
 
-    part = 1 + 0.2 + 0.2  # stim duration + delay up front and a part of the "delay" at the back
+    part = 0.4 + 0.2 + 0.2  # stim duration + delay up front and a part of the "delay" at the back
     idx = int(part/step)
 
     ax.plot(f_trace_times[-1][:idx], f_traces[-1][:idx])
-    strength = 200
-    smoothed = np.convolve(f_traces[-1][:idx], np.ones(strength)/strength)
-    ax.plot(f_trace_times[-1][:idx], smoothed[int(strength/2):idx + int(strength/2)])
+    # strength = 200
+    # smoothed = np.convolve(f_traces[-1][:idx], np.ones(strength)/strength)
+    # ax.plot(f_trace_times[-1][:idx], smoothed[int(strength/2):idx + int(strength/2)])
     ax.set_xlim((-0.2, part-0.2))
     ylim = ax.get_ylim()
     ax.set_ylim((0, ylim[1]))
@@ -132,16 +129,15 @@ def p_unit_example():
     # ax.set_ylim((-1, 1))
     ax.set_xlabel("Contrast")
     ax.set_ylabel("Frequency in Hz")
-    ax.set_title("FI-Curve")
+    ax.set_title("f-I Curve")
 
     plt.tight_layout()
-    plt.savefig("thesis/figures/p_unit_example.png")
+    plt.savefig("thesis/figures/p_unit_example.png", transparent=True)
     plt.close()
 
 
 def fi_point_detection():
-    # cell = "data/final/2013-04-17-ac-invivo-1"
-    cell = "data/final/2012-04-20-af-invivo-1"
+    cell = EXAMPLE_CELL
     cell_data = CellData(cell)
 
     fi = FICurveCellData(cell_data, cell_data.get_fi_contrasts())
@@ -151,7 +147,7 @@ def fi_point_detection():
 
     f_trace_times, f_traces = fi.get_mean_time_and_freq_traces()
 
-    part = 1 + 0.2 + 0.2  # stim duration + delay up front and a part of the "delay" at the back
+    part = 0.4 + 0.2 + 0.2  # stim duration + delay up front and a part of the "delay" at the back
     idx = int(part / step)
     f_zero = fi.get_f_zero_frequencies()[-1]
     f_zero_idx = fi.indices_f_zero[-1]
@@ -166,8 +162,9 @@ def fi_point_detection():
     # mark stim start and end:
     stim_start = cell_data.get_stimulus_start()
     stim_end = cell_data.get_stimulus_end()
-    axes[0].plot([stim_start]*2, (0, fi.get_f_baseline_frequencies()[0]), color="darkgrey")
-    axes[0].plot([stim_end]*2, (0, fi.get_f_baseline_frequencies()[0]), color="darkgrey")
+    axes[0].plot([stim_start, stim_end], (100, 100), color="black", linewidth=3)
+    # axes[0].plot([stim_start]*2, (0, fi.get_f_baseline_frequencies()[0]), color="darkgrey")
+    # axes[0].plot([stim_end]*2, (0, fi.get_f_baseline_frequencies()[0]), color="darkgrey")
     axes[0].set_xlim((-0.2, part - 0.2))
     ylimits = axes[0].get_ylim()
 
@@ -191,12 +188,12 @@ def fi_point_detection():
 
     axes[1].set_xlabel("Contrast in %")
     # axes[1].set_ylabel("Frequency in Hz")
-    axes[1].set_title("FI-Curve")
+    axes[1].set_title("f-I Curve")
     axes[1].set_ylim((0, ylimits[1]))
 
 
     plt.tight_layout()
-    plt.savefig("thesis/figures/f_point_detection.png")
+    plt.savefig("thesis/figures/f_point_detection.png", transparent=True)
     plt.close()
 
 def data_fi_curve():

thesis/Masterthesis.tex

@@ -178,7 +178,7 @@ The stimuli used during the recordings were presented from two vertical carbon r
 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. To be able to cause an amplitude modulation (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. 
+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)) 
@@ -244,7 +244,7 @@ vs = |\frac{1}{n} \sum_n e^{iwt_i}|
 \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 x intervals. This is calculated as,  
+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}}
@@ -264,16 +264,26 @@ Finally the ISI-histogram was calculated within a range of 0--50\,ms and a bin s
 \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)}\hspace{50mm}\raisebox[10]{70mm}{\large\sffamily B)}
 \includegraphics[trim={10mm 5mm 10mm 5mm}, scale=0.8]{figures/f_point_detection.png}
 
-\caption{\label{fig:f_point_detection} \todo{place right in text}On the left: The averaged response of a cell to a step in EOD amplitude. The beginning (at 0\,s) and end (at 1\,s) of the stimulus are marked by the gray lines. The detected values for the onset ($f_0$) and steady-state ($f_{inf}$) response are marked in \todo{color}. $f_0$ is detected as the highest deviation from the mean frequency before the stimulus while $f_{inf}$ is the average frequency in the 0.1\,s time window, 25\,ms before the end of the stimulus. On the right: 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.}
+\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}
 
 
-The adaption behavior of the cell was characterized by the fi-curve consisting of the onset ($f_0$) and steady-state ($f_{inf}$) response. First the ISI frequency trace for each stimulus was calculated. The ISI frequency of a time point t is defined as the $1/T_i$ with $T_i$ the ISI the time point t falls into. This gives a frequency trace starting by the first spike and ending at the last spike. For the further analysis all trials done for a specific contrast were pointwise averaged after cutting them to the same length. This gives an averaged step response for each contrast as seen in figure \ref{fig:f_point_detection} on the left. 
-In this frequency trace the onset $f_0$ and steady-state $f_{inf}$ response were detected (fig. \ref{fig:f_point_detection}). $f_0$ was defined as the farthest deviation from the mean frequency before the stimulus in the range of 25\,ms after stimulus onset. If there was no deviation farther than the variations before the stimulus, then the average frequency in that time window was used. This approximation made the detection of $f_0$ more stable for small contrasts and traces with high variation.   
-The $f_{inf}$ response was defined as the average frequency of the trace in the 0.1\,s time window, 25\,ms before the end of the stimulus. 
-Afterwards a Boltzmann was fitted to the onset response and a rectified line was fitted to the steady-state responses (FI-Curve fig. \ref{fig:f_point_detection}).
+The adaption behavior \todo{explain adaption behavior/ ref the introduction} of the cell was characterized by the fi-curve measurements result in the onset ($f_0$) and steady-state ($f_{\infty}$) response. 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 \todo{with the time resolution of the sampling rate} after cutting them to the same length. 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 onset $f_0$ and steady-state $f_{\infty}$ response were detected. \todo{explain baseline frequency} $f_0$ was defined as the largest deviation from the baseline frequency before the stimulus, within the first 25\,ms after stimulus onset. If there was no deviation farther than the minimum or maximum before the stimulus, then the average frequency in that 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 (FI-Curve fig. \ref{fig:f_point_detection} B).
 
 
 
@@ -283,7 +293,7 @@ Afterwards a Boltzmann was fitted to the onset response and a rectified line was
 % 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 simplest commonly used neuron model is the perfect integrate-and-fire (PIF) model. It'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 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 neither a memory of previous spikes that could cause the negative serial correlation between successive spikes nor can it show any adaption behavior. 
+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 simplest commonly used neuron model is the perfect integrate-and-fire (PIF) model. It'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 neither a memory of previous spikes that could cause the negative serial correlation between successive spikes nor can it show any adaption behavior. 
 
 The next slightly more complex model is the leaky integrate-and-fire (LIF) model. As the name suggests it adds a leakage current to the PIF and as follows the equation \ref{eq:basic_voltage_dynamics} (fig. \ref{fig:model_comparison} LIF). The leakage current adds sub threshold behavior to the model but still cannot reproduce the adaption or serial correlation.
 
@@ -299,16 +309,16 @@ To reproduce this behavior the model needs some form of memory of previous spike
 	\label{eq:adaption_dynamics}
 \end{equation}
 
-It is modeled as an exponential decay with the time constant $\tau_A$ and a strength called $\Delta_A$. $\Delta_A$ is multiplied with the sum of events in the spike train ($\delta (t)$) 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 current of the from equation \ref{eq:basic_voltage_dynamics}  can thus be split into three currents for the modeling of the neuron:
+It is modeled as an exponential decay with the time constant $\tau_A$ and a strength called $\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 current of the from equation \ref{eq:basic_voltage_dynamics}  can thus be split into three currents for the modeling of the neuron:
 
 \begin{equation}
 	I = \alpha I_{Input} - I_A + I_{Bias}
 	\label{eq:currents_lifac}
 \end{equation}
 
-The stimulus current $I_{Input}$, the bias current $I_{Bias}$ and the already discussed adaption current $I_A$. 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}. $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}). 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. $I_{Bias}$ is the bias current that causes the cells spontaneous spiking.   
+The stimulus current $I_{Input}$, the bias current $I_{Bias}$ and the already discussed adaption current $I_A$. 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}. $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}). 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. $I_{Bias}$ is the bias current that adjusts the cells spontaneous spiking.   
 
-Finally noise and an absolute refractory period were added to the model. The noise $\xi$ is drawn in 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.
+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.
 
 
 \begin{figure}[H]
@@ -318,15 +328,14 @@ Finally noise and an absolute refractory period were added to the model. The noi
 
 \end{figure}
 
-Together this results in the dynamics seen in equation \ref{eq:full_voltage_dynamics}. Not shown in the equation is the refractory period $t_{ref}$ and the $\tau_{dend}$ \todo{extra function describing $I_{input}$? rectify and filtering}
+Together this results in the dynamics seen in equation \ref{eq:full_voltage_dynamics}. Not shown in the equation is the refractory period $t_{ref}$ and the $\tau_{dend}$ \todo{move to the start.}.
+
+\begin{align}
+	\tau_m \frac{dV}{dt} &= -V+I_{Bias} +\alpha V_d - 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_d \frac{dV_d}{dt} &= -V_d + \lfloor V_{stim} \rfloor_0 \label{eq:full_model_dynamics_dendrite} 
+\end{align}
 
-\begin{equation}
-\begin{split}
-	\tau_m \frac{dV}{dt} &= -V+I_{Bias} +\alpha I_{Input} - I_{A} + \sqrt{2D}\frac{\xi}{\sqrt{\Delta t}} \\
-	\tau_A \frac{dI_A}{dt} &= -I_A + \Delta_A \sum \delta (t)
-	\label{eq:full_voltage_dynamics}
-\end{split}
-\end{equation}
 
 
 \begin{figure}[H]
@@ -349,7 +358,7 @@ Together this results in the dynamics seen in equation \ref{eq:full_voltage_dyna
  $\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 variables of the model that are fitted.}
+\caption{\label{tab:parameter_explanation} Overview about all parameters of the model that are fitted.}
 \end{table}
 
 
@@ -360,7 +369,7 @@ Together this results in the dynamics seen in equation \ref{eq:full_voltage_dyna
 
 % describe used errors 
 % describe Nelder Mead 
-This leaves the eight variables to be fitted to the cell. During the fitting and the analysis all models were simulated with at time step of 0.05\,ms.
+This leaves the eight parameters to be fitted to the cell. During the fitting and the analysis all models were simulated with at time step of 0.05\,ms.
 The stimuli as 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 sin function of the appropriate frequency, but it was decided to keep the amplitude of the sin at one to make the models more comparable. Changes in the amplitude can be compensated in the model by changing the input scaling factor and the time constant of the dendritic low-pass filter, 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 0.5\,s delay, 0.5\,s duration and 0.5\,s recovery time. The contrasts were the same as in the cell recording. The step stimuli were repeated 8 times. 
@@ -382,6 +391,10 @@ err_i = |1 - ((x^M_i - x^C_i) / x^C_i)| * c_i
 
 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 = (\langle (x^M_i - x^C_i)²\rangle) * c_i
+\end{equation}
+
 All errors were then weighted and 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.