1 Preface

It is still a great challenge to come up with quantitative assertions about complex biological systems. Especially, if one aims towards a functional understanding at the cell, tissue or organismic level, it typically leads to quite involved modells. Analytical progress can be made with simplifying assumptions and the restriction to limiting cases.

This course will discuss a potpourri of mathematical methods ranging from analytical techniques to numerical methods and inform the participant on how to apply them to answer biological question and have enormous fun in doing so.

The mathematical techniques encompass stochastic systems, numerical bifurcation analysis, information theory, perturbation theory and Fourier analysis. The biological examples are chosen mostly from neurobiology, sensory ecology, but also bacterial communication and evolution.

1.1 (Dis)claimer

This scriptum shall evolve. It will have flaws (didactic ones and outright errors) but you, the reader, student, search bot, provide the selective preasure to improve it. Please write to the lecturer if you find fawlt of any kind with it.

2 The Floquet theory of the aktion potential

2.1 Periodic event

Of many natural systems we observe only some drastic events.

The clock regular pulses of pulsars [@]

Action potentials

There is a voltage threshold. The dynamics is described by conductance based equations.

Wake up time

Malatonin level below threshold (suprachiasmatic nucleus)

cell division

Gap phase 1 $\to$ synthsis $\to$ gap phase 2 $\to$ mitosis

2.2 Tonic spikes are limit cycles

Many intersting differential equations describing biological dynamics can be cast into the following form

1. $\dot{\vec x} = \vec F(\vec x,\vec\theta) + \vec\eta(\vec x,t)$

where $\vec x$ is a vector of state variables, e.g., voltage and kinetic variables. $\vec F$ is the time independent (homogeneous) flow field. $\vec\eta$ is a (often assumed small), possibly time dependent (inhomogeneous) perturbation to the system. $\vec\theta$ are system’s parameter. The change of solutions in response to variations in system’s parameter will be studied later.

Example: Hodgkin & Huxley equations

From Kirchhof’s law of current conservation at the membrane follows a dynamics for the membrane voltage $c\dot v+\sum_k^{K}I_k^\mathrm{ion}(v,\vec m_k)={I_{\mathrm{in}}}$. The first term is the capacitive (pseudo current), followed by a sum of corrents from voltage dependent ion channels and the applied input current ${I_{\mathrm{in}}}$. The ion channel dynamics is governed by kinetic equations, i.e., chemical reactions between open and closed states. The reactions often follows first order kinetic equations ${\underline \tau}_k(v)\dot m_k = {\underline M}_k^{(\infty)}(v) - \vec m_k$. ${\underline \tau}_k(v)$ is the matrix of kinetic time constants and ${\underline M}_k^{(\infty)}(v)$ are the steady state activation curves. ${\underline \tau}$ and ${\underline M}^{(\infty)}(v)$ are typically diagonal if the channels are independent. These equations will be derived more formally in Sec. XX. The state vector is then $\vec x=[v,m_{11},...,m_{K1},...]^{\dagger}$. In the absence of perturbations Eq. (1) becomes a homogeneous equation

2. $\dot{\vec x} = \vec F(\vec x,\vec\theta)$

and we assume the existance of a ${P}$-periodic limit cycle solution, ${\vec x_{\mathrm{LC}}}(t)={\vec x_{\mathrm{LC}}}(t+{P})$, also known as periodic orbit.

A neurobiological dogma states that action potentials follow the all-or-nothing principle¹. This can be construed to mean that the exact shape of the action potential does not matter² and that information is stored in the exact spike times³. It also means that the action potential should be stable under perturbations, but that inputs ought to be able to shift the occurrence of spikes in time. If you want, from the dogma and the intend to code into spike times follow two desiderata:

1. The action potential must be stable in amplitude,
2. yet neutrally stable in phase.

To wit, II. just means that stimulus induced phase shifts should neither decay nor blow up.

2.3 Limit cycle stability

Stability of is probed by studying small perturbation to an invariant set solution. Our invariant set is the limit cycle (periodic orbit) ${\vec x_{\mathrm{LC}}}$. Assuming there was a small perturbation to the system the solution can be decomposed as

3. $\vec x(t)={\vec x_{\mathrm{LC}}}(t)+\vec y(t)$

with $\forall t:\|y(t)\|<\epsilon$ some “small” perturbation to the orbit. What small, i.e., $\epsilon$ is we do not want to say now, maybe later, lets see …

Assuming the perturbation was transient (only in initial conditions) and the system is homogeneous again we plug the Ansatz of Eq. (3) into Eq. (2) and get

$\frac{\mathrm{d}}{{\mathrm{d}}t}{\vec x_{\mathrm{LC}}}(t)+\dot{\vec y}(t)=\vec F({\vec x_{\mathrm{LC}}})+\underbrace{\nabla\vec F({\vec x_{\mathrm{LC}}})}_{{{{\underline J}}}({\vec x_{\mathrm{LC}}}(t))} \cdot \vec y(t)$

The Jacobi matrix evaluated on the limit cycle can be written as a function of time ${{{\underline J}}}(t)=\nabla\vec F({\vec x_{\mathrm{LC}}}(t))$. Note that since the limit cycle solution is ${P}$-periodic, so is ${\underline J}(t)={\underline J}(t+{P})$.

Identifying the limit cycle solution above we are left with the first variational equation of Eq. (2)

4. $\dot{\vec y}(t)={{{\underline J}}}(t)\vec y(t).$

Hence, one needs to study of linear system with periodic coefficients. One solution of Eq. (4) can be guessed, let us try the time derivative of the orbit, $\frac{{\mathrm{d}}{\vec x_{\mathrm{LC}}}}{{\mathrm{d}}t}{}$

5. $\frac{\mathrm{d}}{{\mathrm{d}}t}(\frac{{\mathrm{d}}{\vec x_{\mathrm{LC}}}}{{\mathrm{d}}t})=\frac{\mathrm{d}}{{\mathrm{d}}t}\vec F({\vec x_{\mathrm{LC}}})=\nabla\vec F({\vec x_{\mathrm{LC}}})\frac{\mathrm{d}}{{\mathrm{d}}t}{{\vec x_{\mathrm{LC}}}}={{{\underline J}}}(t)\frac{{\mathrm{d}}{\vec x_{\mathrm{LC}}}}{{\mathrm{d}}t}.$

So it is a solution alright, and it happens to be a ${P}$-periodic solution. This solution is called the Goldstone mode. But for arbitray intitial conditions not all solutions should be periodic.

Def: Floquet Ansatz

According to Floquet theory the solution to Eq. (4) can be written in the form of a ${P}$-periodic similarity matrix⁴ and an matrix exponential

6. $y(t)={\underline R}(t)e^{t{\underline \Lambda}}.$

For the proof please consult (Chicone 2006), here we are just goint to work with this as an Ansatz. The constant matrix ${\underline \Lambda}$ is called the Floquet matrix.

Recall the matrix exponential

Def: Matrix Exp

Let ${\underline A}\in{I\!\!\!\!C}^{n\times n}$ then $\exp{\underline A}=\sum_{k=0}^\infty\frac1{k!}{\underline A}^k$ A useful corollary of this definition is that the eigen vectors of an exponentiated matrx are the same as those of the original and eigenvalues become exponentiated. If $\lambda_i$, $\vec w_i$ are the eigenvalue, eigenvector pairs of the matrix ${\underline A}$, i.e., ${\underline A}\vec w_i=\lambda_i\vec w_i$ then by using this identity $k$-times

7. $e^{{\underline A}}\vec w_i=\left(\sum_{k=0}^\infty\frac1{k!}{\underline A}^k\right)\vec w_i=\sum_{k=0}^\infty\frac1{k!}{\underline \lambda}_i^k\vec w_i=e^{\lambda_i}\vec w_i$

Inserting the Floquet Ansatz into Eq. (4) yields

$\dot{{\underline R}}e^{t{\underline \Lambda}}+{\underline R}(t){\underline \Lambda} e^{t{\underline \Lambda}}={\underline J}(t){\underline R}(t)e^{t{\underline \Lambda}}.$

Which by multiplying with $e^{-t{\underline \Lambda}}$ results in a dynamics equation for similarity matrix ${\underline R}$.

8. $\dot{{\underline R}}={\underline J}(t){\underline R}(t)-{\underline R}(t){\underline \Lambda}.$

Remember that ${\underline R}$ was invertable so one can also derive an equation for the inverse

9. $\frac{\mathrm{d}}{{\mathrm{d}}t}{\underline R}^{-1}={\underline R}^{-1}\dot{{\underline R}}\,{\underline R}^{-1}={\underline \Lambda}\,{\underline R}^{-1}(t)-{\underline R}^{-1}(t){\underline J}(t).$

2.3.1 Eigensystem of the Floquet matrix

The Floquet matrix, ${\underline \Lambda}$, is constant, though not necesarily symmetric. Hence it has an orthonormal left and right eigensystem

10. $\Lambda\vec w_i=\mu_i\vec w_i$ and $\vec z_i\Lambda=\mu_i\vec z_i$ with $\vec z_i\cdot\vec w_j=\delta_{ij}$

One can define a “rotated” eigensystem

Def: (“Rotated” eigensystem):

$\vec W_i(t)={\underline R}(t)\vec w_i$ and $\vec Z_i(t)=\vec z_i{\underline R}^{-1}(t)$. For which ipso facto obeys

11. $\forall t:\vec Z_i(t)\cdot\vec W_j(t)=\delta_{ij}$.

$\vec W_i(t)$ and $\vec Z_i(t)$ are also called the Floquet modes.

If we project the eigenvectors on the Eqs. (8) and (9) and use this definitions we get

12. $\frac{\mathrm{d}}{{\mathrm{d}}t}{\vec W_k}=({\underline J}(t)-\mu_k{\underline I})\vec W_k(t)$

and

13. $\frac{\mathrm{d}}{{\mathrm{d}}t}{\vec Z_k}=(\mu_k{\underline I}-{\underline J}^{\dagger}(t))\vec Z_k(t).$

If one projects the general solution $y(t)$ from Eq. (6) on the the adjoint Floquet modes

$\vec Z_k(t)\cdot\vec y(t)=\vec z_k{\underline R}^{-1}(t){\underline R}(t)e^{t{\underline \Lambda}}=\vec z_ke^{t{\underline \Lambda}}=\vec z_ke^{t\mu_k}$

If $\nu_k<0$ the perturbation decays exponentially in this rotated coordinate frame.

Poincare Section

bla

Note that if $\nu_0=0$, then according to Eq. (5)

14. $\vec W_0(t)=\frac{\mathrm{d}}{{\mathrm{d}}t}{\vec x_{\mathrm{LC}}}(t)$

is the Goldstone mode, and from Eq. (11)

15. $\vec Z_0(\phi)\cdot\frac{\mathrm{d}}{{\mathrm{d}}t}{\vec x_{\mathrm{LC}}}(\phi)=1$

2.4 Neutral dimension and phase shifts

The evolution of the phase of Eq. (1) is given by

$\frac{{\mathrm{d}}\phi}{{\mathrm{d}}t}=\nabla\phi(\vec x)\cdot\frac{{\mathrm{d}}\vec x}{{\mathrm{d}}t}$.

To first order this is

$\frac{{\mathrm{d}}\phi}{{\mathrm{d}}t}=\nabla\phi({\vec x_{\mathrm{LC}}})\cdot\frac{{\mathrm{d}}{\vec x_{\mathrm{LC}}}}{{\mathrm{d}}t} =\nabla\phi\cdot\vec F({\vec x_{\mathrm{LC}}})+\nabla\phi\cdot\vec\eta({\vec x_{\mathrm{LC}}},t)$.

There are several ways to define a phase (Hilbert transform, linear interpolation, …). A desiderata could be to have a linear phase increase in the unperturbed case ($\vec\eta=0$), say $\phi(t)={f_0}t$. [… proto-phase]. From this desiderata it follows that one must have $\forall t:\nabla\phi\cdot\vec F({\vec x_{\mathrm{LC}}})={f_0}$. Given Eq. (16), this is easily achieved with the following identification

16. $\nabla\phi=\vec Z_0(\phi){f_0}=\vec Z(\phi)$

The input-output (I/O) equivalent phase oscillator to Eq. (1) can then be written as

17. $\dot\phi={f_0}+ \vec Z(\phi)\cdot \vec\eta(t)$

A spike-train can be written as

18. $y(t)=\sum_k\delta(t-{t^{\mathrm{sp}}})=\sum_k\delta(\phi(t)-k)$

3 Fourier theory

To proceed with the analysis we need some results from Fourier theory.

3.1 The Fourier base

One may think of the Fourier transform as a projection onto a new basis $e_\omega(t)=e^{{\mathrm{i}}\omega t}$. Define the projection of a function as

$F(\omega)={\langle}e_w\cdot f{\rangle}=\int_{-\infty}^\infty dt\,e_\omega^*(t)f(t)$

The inverse Fourier transform is a projection onto $e_t^*(\omega)$

$f(t)={\langle}e_t^*\cdot F{\rangle}=\int_{-\infty}^\infty d\omega\,e_t(\omega)F(\omega)$

3.2 Existence of the Fourier integral

Usually all transcendental pondering about the existance of mathematical objects is the subject of pure math and should not bother us too much (we trust the work as been done properly). But in the Fourier transform case it motivates a subject we need to address mainly because of $\delta$-functions, which are heavily used in theoretical neurobiology, because they are reminestcent of action potentials or indicate the time of such a spike event. What does it mean for the Fourier transform to exist? It means that the integrals invoved in the definition of the Fourier transform converge to finite values. OK. So let us look at the magnitude of the transform of a function $f$. It can be upperbound by the Cauchy-Schwarz inequality

$|{\langle}e_w\cdot f{\rangle}|=\left|\int_{-\infty}^\infty dt\; e^{-{\mathrm{i}}\omega t}f(t)\right| \leqslant \int_{-\infty}^\infty\underbrace{|e^{-{\mathrm{i}}\omega t}|}_{=1} |f(t)|\;dt= \int_{-\infty}^\infty |f(t)|\;dt.$

This means that if one assumes the function $f(t)$ is absolute integrable

$\int_{-\infty}^\infty|f(t)|\;dt<\infty,$

then the Fourier integral exists – hurray. Of course, the same works for the integral involved in inverse Fourier transform. Note that this is an implication in one dirrection. All functions satisfying absolute integrability are lumped together in $L^1({I\!\!R})$.

The bad news is that applying a Fourier transform to one of the members of $L^1({I\!\!R})$ can through you out of it. And then? Well, luckily there is the set of Schwartz functions, wich is closed under the Fourier transform.

3.2.1 Orthonormality

The issue that absolute integrability is insufficient already manifests when trying to Fourier transform it own basis, since $\int_{-\infty}^\infty dt|e^{{\mathrm{i}}\omega t}|=\infty$. Lets give it a name anyway

${\langle}e_\omega\cdot e_\nu{\rangle}=\int_{-\infty}^\infty dt\,e^{{\mathrm{i}}(\nu-\omega)t}=\delta(\nu-\omega)$

It follows that the delta function is symmetric

$\delta(\omega)=\int_{-\infty}^\infty dt\,e^{{\mathrm{i}}\omega t}=\int_{-\infty}^\infty dt\,e^{-{\mathrm{i}}\omega t}=\delta(-\omega)$

Note also that

19. ${\langle}e_\omega\cdot 1{\rangle}=\delta(\omega)$

3.2.2 Convolution

The convolution of two function is defined as

$h(t)=(f*g)(t)=\int_{-\infty}^\infty dr f(t-r)g(r)$

$H(\omega)=\int dt\,e^{{\mathrm{i}}\omega t}\int dr\,f(t-r)g(r)=\int dr\,g(r)\int dt\,e^{-{\mathrm{i}}\omega(t+r)}f(r)$

$=\int dr\,g(r)e^{-{\mathrm{i}}\omega r}\int dt\,e^{-{\mathrm{i}}\omega r}f(t)=G(\omega)F(\omega)$

3.2.3 Derivative

What is $\delta'(t)$?

$(\delta'*f(t))=\int\limits_{-\infty}^\infty\delta'(r)f(t-r)dr=[\delta(r)f(t-r)]_{r=-\infty}^\infty-\int\limits_{-\infty}^\infty\delta(r)f'(t-r)dr=-f'(t)$

Hence

20. $(\delta'*f(t))=-f'(t)$

4 The continuum limit of a membrane patch

Motivation and aim

In this lecture ion channels are introduced as stochastic devices floating in the membrane of a nerve cell. It should motivate why the analysis techniques and models introduced in this lecture need to deal with fluctuations and noise. (Nerve) cells produce stochastic processes on several levels:

1. The ion channels in their membrane stochastically jump beteen conformations, governed by Master equations.

2. On a more macroscipic level the their membrane voltage fluctuations show properties of coloured noise, well described by diffusion processes and stochastic differential equations.

3. The trains of action potentials they emmit from point processes.

In the stationary state these processes can be subjected to spectral analysis.

Note that in 1952, the first equations describing membrane voltage dynamics where the deterministic rate equations by Hodgin & Huxley (Hodgkin and Huxley 1952). Only in 1994 Fox and Lu derived these equations from the continuum limit of an ensemble of stochastic ion channels (Fox and Lu 1994). Essentially by doing a diffusion approximation.

Since then the nuences of appllying diffusion approximations to neurons have been investigated (Linaro, Storace, and Giugliano 2011,Orio and Soudry (2012)) and reviewed (Goldwyn et al. 2011,Goldwyn and Shea-Brown (2011),Pezo, Soudry, and Orio (2014)).

4.1 The ion channel as a Markov model

Proteins change conformation on various triggering signals:

• a change in pH
• a change in the surrounding electric field
• mechanical preassure

Such dynamics can be described by a finite state Markov model which is mathematically described as a Master equation.

Starting with a simple ion channel that has an open conformation, $O$, in which ions can pass (say K⁺ ions) and a cloased states, $C$, which blocks ion flow

$O\underset{\beta}{\overset{\alpha}{\rightleftharpoons}}C$

We can define a vector of the probability of being open and closed, $\vec p=[p_O,p_C]^{\dagger}$, respectively.

An ionic current produced in the open state would be

$I^\mathrm{ion}=\gamma_{\mathrm K^+}N_O(E_{\mathrm K^+}-v)$

$\gamma_{\mathrm K^+}$ and $E_{\mathrm K^+}$ are the unitary conductance and the Nernst potential respectively. The average such current would be

${\langle}I^\mathrm{ion}{\rangle}=\gamma_{\mathrm K^+}Np_O(E_{\mathrm K^+}-v)$

where $N$ is the total number of channels in the membrane patch under consideration. But what about a particular stochastic realisation of the current, what about the fluctuations around the mean?

If we have $N$ channels than the number of $k$ of them being open is binomially distributed

$P_O(k)={N\choose k}p_O^kp_C^{N-k}$

In actuality the channels are part of the membrane dynamical system, where $\alpha$ and $\beta$ depend at least on $v$ and hence are not constant during a spike. We need an update rule how to get from the probability of being open at time $t$ to the probability of begin open at time $t+dt$. This is given by

$\begin{pmatrix}p_O(t+dt)\\p_C(t+dt)\end{pmatrix}=\begin{pmatrix}1-\alpha\,dt&\beta\,dt\\\alpha\,dt&1-\beta\,dt\end{pmatrix}\begin{pmatrix}p_O(t)\\p_C(t)\end{pmatrix}$

or in vector form

$\vec p(t+dt)=\left({\underline I}+{\underline Q} dt\right)\vec p(t)$

The infinitesimal limit is

$\frac{\mathrm{d}}{{\mathrm{d}}t}\vec p={\underline Q}\vec p$

With $p_C=1-p_O$ we can express one row of this equation as

$\dot p_O=\alpha p_O - \beta(1-p_O)$

or

$\tau\dot p_O=p_O^{(\infty)}-p_O$ with $\tau=\frac1{\alpha+\beta}$ and $p_O^{(\infty)}=\frac\alpha{\alpha+\beta}$

which has solution

$p_O(t)=p^{(\infty)}_O(1-e^{-t/\tau})$

Fourier transformation leads to

$\tau{\mathrm{i}}\omega\tilde p(\omega)=p^{(\infty)}_O-\tilde p(\omega)$

or

$\tilde p(\omega)=\frac{p^{(\infty)}}{1+{\mathrm{i}}\tau\omega}$

$|\tilde p(\omega)|^2=\tilde p(\omega)\tilde p^*(\omega)=\frac{p^{(\infty)}}{1+(\tau\omega)^2}$

a Lorenzian spectrum. Inverse Fouerier transform yields the covariance function

$c(\tau)={\langle}p_O(t)p_O(t+\tau){\rangle}$

4.1.1 The $n$-state channel

Let us try to calculate the statistics of the current originating from an $n$-state channel (just like in the two state case). Why would one do this? The idea, later on is to be able to find a continuous stochastic process that we can simulate and analyise easily.

Let $K(t)\in[1,...,n]$ be the realisation of the $n$-state Markov channel. For example a K⁺-channel with four subunits

$1 \underset{\beta}{\overset{4\alpha}{\rightleftharpoons}} 2 \underset{2\beta}{\overset{3\alpha}{\rightleftharpoons}} 3 \underset{3\beta}{\overset{2\alpha}{\rightleftharpoons}} 4 \underset{4\beta}{\overset{\alpha}{\rightleftharpoons}} 5$

Assuming the ion channel has $n$ conformations, of which one is conducting, let us further define

$G=\delta_{1K(t)}=\begin{cases}1&:K(t)=1\\0&:K(t)>1.\end{cases}$

The single channel current at a voltage clamp $v(t)=v$ is then

21. $I(t)=\gamma G(t)(E-v).$

How does it evolve? Define $p_{i}(t)=P(K(t)=i)$ and $\vec p(t)=[p_1,...,p_n]^{\dagger}$, then

22. $\frac{\mathrm{d}}{{\mathrm{d}}t}\vec p={\underline Q}\vec p$

With formal solution

$\vec p(t)=e^{{\underline Q}t}\vec p(0)={\underbrace{(e^{{\underline Q}})}_{{\underline M}}}^t\vec p(0)={\underline M}^t\vec p(0)$

Use the singular value decompostion, ${\underline M}={\underline U}\,{\underline \Sigma}\,{\underline V}^{\dagger}$, the matrix power can be written as ${\underline M}^t={\underline U}\,{\underline \Sigma}^t{\underline V}^{\dagger}$. Or (recall Eq. (7))

23. $\vec p(t)=\sum_{k=1}^n\vec u_k\vec v_k^{\dagger}e^{\nu_kt}\vec p(0)$.

If Eq. (22) has a stationary distribution in the $t\to\infty$ limit, then this must correspond to the eigenvalue $\nu_1=0$ (let us assume they are ordered). So $\frac{\mathrm{d}}{{\mathrm{d}}t}\vec p(\infty)=0\Longrightarrow$ for $\vec p(\infty)=\vec v_1,\exists\nu_1=0:{\underline Q}\vec v_1=\nu_0\vec v_1=0$. Therefore, the solution can be written as

24. $\vec p(t)=\vec p(\infty)+\sum_{k=2}^n\vec u_k\vec v_k^{\dagger}e^{\nu_kt}\vec p(0)$.

The average channel current of Eq. (21) is

${\langle}I(t){\rangle}=\gamma(E-v)\sum_{k=1}^n p_k(t)\delta_{1k}=\gamma(E-v)(p_1(\infty)+\sum_{k=2}^nu_{1k}v_{1k}p_1(0)e^{\nu_kt})$,

which if the chain is stable ($\nu_k<0:\forall k>1$) has steady state

${\langle}I{\rangle}=\lim_{t\to\infty}{\langle}I(t){\rangle}=p_1(\infty)$.

The steady state covariance $C(\Delta)=\lim_{t\to\infty}$ in this case is

$C_t(\Delta)={\langle}I(t)I(t+\Delta){\rangle}-{\langle}I{\rangle}^2=\gamma^2(E-v)^2\sum_{j,k=1}^n\delta_{1j}\delta_{1k}\,p_j(t)p_k(t+\Delta)-{\langle}I{\rangle}^2$

Hence

$C_t(\Delta)=p_1(t)p_1(t+\Delta)-{\langle}I{\rangle}^2=\sum_{i,j=2}^n[\vec u_i \underbrace{\vec v_i^{\dagger}\vec v_j}_{=\delta_{ij}}\vec u_i^{\dagger}]_{11}e^{\nu_it+\nu_j(t+\Delta)}$

If $\nu_k<0:\forall k>1$ then

25. $C(\Delta)= \sum_{i=2}^nu_{1i} u_{1i}e^{\nu_i\Delta}$

is a sum of exponentials. The spectral density is a superposition of Lorenzians.

4.1.2 Simulation a Jump process

Do not track individual channels but the channel numbers. Starting in a particular state $\vec N=[N_1,...,N_n]$ at time $t$ the life time of staying in that state until $t+\tau$ is

$f(\tau)=\lambda e^{-\lambda\tau}$

The escape rate of the state is (any reaction occuring)

$\lambda=\sum_{k=1}^nN_ka_k$

where $a_k$ are the rate of leaving state $k$. For example $a_3$ in the K⁺-channel is

$a_3=2\beta+2\alpha$

But which reaction did occur? Let $j$ be the reaction (not the state!). For example, there are 8 reactions in the K⁺ channel. The probabilities of occurance associated with any one of them is

$p(j)=\frac{N_j\zeta_j}{\sum_{k=1}^{N_{\mathrm{reac}}}N_k\zeta_k}=N_j\zeta_j/\lambda$

In a computer algorithm one can draw the waiting time for the next reaction simply by generating a uniform random number $r_1\sim U(0,1)$ and then

$\tau\leftarrow\ln(r_1^{-1})/\lambda$

The reaction is determined with a second random number $r_2\sim U(0,1)$

$P(j)=\sum_{k=1}^jp(k)$

$j\leftarrow\text{argmax}_{j}P(j)<r_2$

while
  r1 = rand
  r2 = rand
  tau = ln(1/r1) / a

4.2 Statistically equivalent diffusion process (Orenstein-Uhlenbeck Process)

The jump process discussed in the previous sections is a continuous time Markov-process on a discrete domain, $K(t)\in{I\!\!N}$ with $t\in{I\!\!R}$. A diffusion process is a continuous time Markov-process on a continuous domain, $\eta(t)\in{I\!\!R}$ with $t\in{I\!\!R}$.

Can are in search for a a diffusion process such that

$I=\gamma\left(p_1(\infty) + \sum_{k=2}^n\eta_k(t)\right)(E-v)$

has the same first and second order statistics (in voltage clamp) as Eq. (21)? Let us try

26. $\tau_k(v)\dot\eta_k=-\eta_k+\sigma(v)\,\xi(t)$ where ${\langle}\xi(0)\xi(\Delta){\rangle}=\delta(\Delta)$

To solve it, use the Fourier Transform

${\mathrm{i}}\omega\tau\tilde\eta(\omega)=-\tilde\eta(\omega) + \sigma\,\chi(\omega)$

The place holder symbold $\chi$ was introduced for the Fourier transform of the stochastic process $\xi(t)$. Rearanging yields

$\tilde\eta=\frac{\sigma\chi(\omega)}{1+{\mathrm{i}}\omega\tau}$

The spectrum is

$\tilde \eta(\omega)\tilde \eta^*(\omega)=\frac{\sigma^2\chi(\omega)\chi^*(\omega)}{1+(\tau\omega)^2}$

By definition $\chi(\omega)\chi^*(\omega)$ is the Fourier Transform of the covariance function $\delta(\Delta)$ and from Eq. (19) this is one. Hence,

$\tilde \eta(\omega)\tilde \eta^*(\omega)=\frac{\sigma^2}{1+(\tau\omega)^2}$

Applying the inverse Fourier transform results in the correlation function

$C(t)=\frac{\sigma^2}\tau e^{-|t|/\tau}$.

A super position of independent such OU process

$\sum_{k=i}^n\eta_i(t)$

leads to a correlation function with the same structure as in Eq. (25). We identify $\tau_i=1/\nu_i$ and $\sigma_i=u_{1i}$.

The idea of matching the second order statistics can be formulated in a far more abstract way in terms of the Kramers-Moyal-van-Kampen expansion of the Master requation

$\dot p=\int w(x'\to x)p(x',t)-w(x\to x')p(x,t) dx'$

$\partial p(x,t)=\sum_{n=1}^\infty\frac{\partial^n}{\partial x^n} K_n(x,t)p(x,t)$

from

5 Information theory for the living

This theory was never ment to be used to describe living systems in which meaning, i.e., the content of a message actually matter. Information theory deals with optimal compression, lossless transmission of signals irrespective of weather it is relevant or a TV program.

Nontheless a look at a quantitative science of communication may be insightfull.

Consult []

Since there is a relation between the PRCs introduced in XXX and the filtering properties of a neuron, one can seek to do the same for phase dynamics and ask what PRC would maximise information transmission. But first, one needs to develop how the PRC predicts a lower bound on the mutual information rate. We begin with a short review of the basic tenets of information theory. Within information theory, a neural pathway is treated as a noisy communication channel in which inputs are transformed to neuronal responses and sent on:

5.1 The communication process

The science of communication is concerned with at least two subtopics: $(i)$ the efficient representation of data (compression); and $(ii)$ the save transmission of data through unreliable channels.

A source (the message) typically runns through the following processing sequence:

One of the formal assertions of information theory is that these two problems can be addressed separately (without loss of generality or efficiency): Meaning first one compresses by removing redundancies. Then one again adds failsafe redundancy to combat the unreliability.

However convenient for engineering, this does not mean that biological systems have to make use of this fact (source coding and channel coding could be very intertwined).

Also note that many of the mathematical results in information theory are bounds, inequalities not achievable in real systems. To get an idea of the mindset of information theory, check-out Shannon’s source coding theorem.

5.2 Source coding, data compression and efficient representation

Remember Morse’s code. Why does the character E only have a single symbol ($\cdot$), while the character Q ($--\cdot\,-$) has so many?

The idea could be that: Highly probable symbols should have the shortest representation, unlikely ones may occupy more space. What we want to compress is called the source and data compression depends on the source distribution. This sounds like a good idea for a biological system: Do not spend resources on rare events. Well not quite. Cicadas hear the sounds of their mating parners only once at the very end of a possibly 19 year long life.

Example (Genetic code):

Cells code 20 amino acids with words composed of a four letter⁵ alphabet $A=\{$A,G,C,T$\}$. Only words of length 3 are sufficient $4^2=16<20<64=4^3$. If nature was to require only 16 amino acides, two character words would be sufficient. Only four to drop: Discrading the 4 least probable AS to occur in Proteins: Tryptophan 1.1%, Methionin 1.7%, Histidin 2.1% and Cystein 2.8%. No code words would make for an error of 7.7%.

In general, we can define a risk $\delta$, i.e., the probability of not having a code word for a letter. Then, to to be efficient, one should choose the smallest $\delta$-sufficient subset $S_\delta\subset A$ such that $p(x\in S_\delta)\geqslant1-\delta$.

Def. (Essential bit content):

For an alphabet, $S_\delta$, the essential bit content,i.e., the number of binary questions asked to identifiy an element of the set is $H_\delta=\log_2|S_\delta|$

In the case of the genetic colde the essential bit content for $\delta = 0.077$ is $H_{0.077}=4$ or if we use the base $b=4$ for the log it is 2, which is the length of the code words required for 16 AS.

Take an even simpler example: We have the alphabet $A=\{0,1\}$ with probability distributions A: $\{p_0,p_1\}=\{\frac12,\frac12\}$ and B: $\{p_0,p_1\}=\{\frac34,\frac14\}$. In figure~ you can see a plot of $\frac1NH_\delta$ over the allowed error $\delta$ for words of different length $N$.

For increasing $N$ the curves depend on $\delta$ to a lesser degree. Even more it converges to a line around the entropy

27. ${H}=-\sum_k p_k \log_2(p_k)$

which is $H_{\text A}=1$ and $H_{\text B}\approx0.81$.

This is closely related concept of the typical set $T^N$. Members of the typical set $(x_1,..,x_N)=\vec x\in T^N$ have $$-\frac1N\log p({\vec x})\approx H(x).$$ Hence their probability is $p({\vec x})\approx 2^{-N\,H(x)}$, for all of them, which implies the cardinality of the typical set is $2^{N\,H(x)}$. There holds a kind of law of large numbers (Asymptotic Equipartition Property, AEP) for the typical set stating that: For $N$ i.i.d. random variables ${\vec x}=(x_1,..,x_N)$ with $N$ large ${\vec x}$ is almost certainly am member of $T^N$ ($p({\vec x}\in T^N)=1-\varepsilon$).

Shannon argued that one should therefore base the compression on the typical set and showed that we can achieve a compression down to $NH$ bits.

Asymptotic Equipartition Property (AEP):

For i.i.d. random variables $$-\frac1n\log p(X_1,..,X_n)=-\frac1n\sum_i\log p(X_i)\to-{\langle}\log p(X){\rangle}=H(X)$$

Note that the i.i.d. assumption is not to restrictive as we represent correlated processes in terms of the i.i.d. coefficients in their transform (or the empirical counter parts).

The typical set $T^\varepsilon_n=\{(x_1,..,x_n):|-1/n\log p(x_1,..,x_n)-H(X)|<\varepsilon\}$ allows Shannon’s source coding algorithm: The encoder checks if the input sequence lies within the typical set; if yes, it outputs the index of the input sequence within the typical set; if not, the encoder outputs an arbitrary in $n(H+\varepsilon)$ digit number.

By formalising this argument Shannon proved that compression rates up to the source entropy is possible. The converse, that compression below is impossible is a bit more involving.

5.3 Channel coding

Here we consider a ,,memoryless’’ channel:

$\overset{\text{message}}{W}\overset{\text{encode}}{\to}X\ni x\to\overset{\text{noisy channel}}{p(y|x)}\to y\in Y\overset{\text{decode}}{\to}\overset{\text{est. Message}}{\hat W}$

The rate $R$ with which information can be transmitted over a channel without loss is measured $\frac{\text{Bits}}{\text{transmission}}$ for a discrete time channel or $\frac{\text{Bits}}{\text{second}}$ for a continuous time channel. Operationally, we wish all bits that are transmitted to be recovered with negligible probability of error.

A measure of information could be:

The average reduction in the number of binary-questions needed to identify $x\in X$ before and after observing $y\in Y$.

This would just be the difference in entropies:

${M}(X;Y)={H}(X)-{H}(X|Y)= -\sum\limits_{x\in X} p(x)\log_2 p(x) + \sum\limits_{x\in X\atop y\in Y}p(x,y)\log_2p(x|y)$.

This goes by the name of mutual information or transinformation. Remember marginalisation

$p(x)=\sum_{y\in Y}p(x,y)$.

So the mutual information is

${M}=-\sum_{x\in X\atop y\in Y} p(x,y)\log_2 p(x) + \sum_{x\in X\atop y\in Y}p(x,y)\log_2\left(\frac{p(x|y)p(y)}{p(y)}\right)$

or

${M}(X;Y)=\sum_{x\in X\atop y\in Y} p(x,y)\log\frac{p(x,y)}{p(x)p(y)}$

From the complete symmetry of this quantity we can also write it as

${M}(X;Y)={H}(Y)-{H}(Y|X)$.

The following Figure illustrates how the mutual information is related to respective (conditional) entropies of the input and output ensemble.

We have heard that $I(X;Y)$ quantifies the statistical dependence of $X$ and $Y$, but how is that related to error-free communication?

$I(X;Y)$ depends on the input ensemble. To focus on the properties of the channel we can simply take an ,,optimal’’ input ensemble and define the channel capacity

$C=\max_{p(x)}I(X;Y).$

It will be left to the sender to actually find the optimal input statistics. Note that $I(X;Y)$ is a concave function ($\cap$) of $p(x)$ over a convex set of probabilities $\{p(x_i)\}$ (this is relevant for procedures like the Arimoto-Blahut algorithm for estimating $C$) and hence a local maximum is a global maximum.

p = []
while
  p(

Shannon established that this capacity indeed measures the maximum amount of error-free information that can be transmitted. A trivial upper bound on the channel capacity is

$C\leqslant\min\{\log|X|,\log|Y|\}$.

This is due to the maximum entropy property of the uniform distribution in the discrete case:

Example (Maximum entropy, discrete case):

For the derivate of the entropy from Eq. (27) one gets: $\frac\partial{\partial p_i}{H}=-\log_2 p_i-1$, which leads to $p_i\propto e^{-1}\quad\forall i$. After normalisation one has $p_i=1/N$, so the uniform distribution maximises the entropy in the discrete case.

5.4 Information transmission (continuous case)

Information theory has long been in the discussion as an ecological theory upon which to judge the performance of sensory processing (Atick 1992). This led Joseph Atick, Horace Barlow and others to postulate to use this theory to study how nervous systems adapt to the environment. The goal is to make quantitative predictions about what the connectivity in the nervous system and the structure of receptive fields should look like, for instance. In this sense, information theory was hoped to become the basis for an ecological theory of adaptation to environmental statistics (Atick 1992,Barlow (1961)).

Some of the signals in nature (and those applied in the laboratory) are continuous in time and alphabet. Does it make sense to extend the definition of the entropy as

28. ${H}(x) =-\int {\mathrm{d}}x\,p(x)\ln p(x)$?

Maybe. Let us see how far one gets with this definition. It is called differential entropy by the way. Through quantisation this can be related back to the entropy of discrete alphabets.

If the $p(x)$ is smooth then one associates the probability of being in $i\Delta\leqslant x\leqslant(i+1)\Delta$ with

$p_i=p(x_i)\Delta = \int_{i\Delta}^{(i+1)\Delta}{\mathrm{d}}x\,p(x)$

The entropy of the quantised version is

${H}_\Delta=-\sum\limits_{i=-\infty}^\infty p_i\ln p_i=-\sum\limits_{i=-\infty}^\infty \Delta p(x_i)\ln(p(x_i)\Delta)$

$=-\sum\limits_{i=-\infty}^\infty \Delta p(x_i)\ln p(x_i)-\ln\Delta$

This is problematic as the second term goes to infinity for small quantisations. Formally, if $p(x)$ is Rieman integrable, then

$\lim_{\Delta\to0}={H}_\Delta+\ln\Delta=-\int {\mathrm{d}}x\,p(x)\ln p(x)$

Since the infinitesimal limit is taken we can also take $n$ to be the number of quantal intervals so that in the limit

$\lim_{\Delta\to0\atop n\to\infty}\ln\Delta\approx n$.

So that an $n$-bit quantisation of a continous random variable $x$ has entropy

${H}(x)+n$.

With the mutual information being the difference of entropies the quantisation term vanishes.

A first example can be taken from rate coding⁶. In terms of our spike trains from Eq. (18), the instantaneous firing rate can be defined as

Example (Optimal rate coding in the retina (Laughlin 1981)):

In short, Simon measured the distribution of light intensities that would impinge on the fly’s compund eye in a natural environment: It was more or less normally distributed. He postulated that an efficient allocation of resources (firing rate) would be, to spread more likely inputs over a broader range of firing rates. Why? It makes them easier to differentiate by upstream neurons. Via experiments he found that the fly’s compound eye approximates the cumulative probability distribution of ctrast lelvels in natural scenes.

Let us pause and ask: What is the maximum entropy distribution for a continuous alphabet?

Example (Maxent, continuous with variance constraint):

Take a random variable $x\in{I\!\!R}$. It can not be the uniform distribution as in the discrete case. In fact we need additional constraints. For example one may ask for the maximum entropy distribution, $p(x)$, given a fixed mean, $\mu$, and variance $\sigma$. Using Lagrange’s multipliers to formulate the constraint optimisation problem

$C=H+\lambda_0\int{\mathrm{d}}x\,p(x)+\lambda_1\left(\int{\mathrm{d}}x\,p(x)x-\mu\right)+\lambda_2\left(\int{\mathrm{d}}x\,p(x)(x-\mu)^2-\sigma^2\right)$

One searches for $p(x)$ that fulfils $\delta C/\delta p(x)=0$ and $\partial C/\partial\lambda_i=0$, where $\delta/\delta p(x)$ is the variantional derivative. One finds

$p(x)=\exp(\lambda_0-1+\lambda_1x+\lambda_2(x-\mu)^2)$.

With the normalisation from $\partial C/\partial\lambda_i=0$ we get the normal distribution ($e^{\lambda_0-1}=1/\sqrt{2\pi\sigma^2}$, $\lambda_1=0$ and $\lambda_2=1/\sigma^2$).

What did we learn? Well

29. $H(x)\leqslant H_\text{Gauss}(x)$.

If we consider the whole spike train $y(t)=\sum_k\delta(t-{t^{\mathrm{sp}}}_k)$ (see Eq. (18)) as the ouput, not just its mean input intensity and mean firing rate, we have a continuous time dependent stochastic processes to deal with. Note that the definition of the spike train, if integrated, is related to the empirical distribution function.

If we average infinitely many trials we get the instantaneous firing rate $r(t)={\langle}y(t){\rangle}_{y|x}$. We will give a mathematical definition of $r(t)$ later. Our communication channel looks like that

The entropy rate of an ergodic process can be defined as the the entropy of the process at a given time conditioned on its past realisations in the limit of large time

30. $H[x]=\lim_{t\to\infty}{H}[x_t|x_\tau:\tau<t].$

The mutual information rate measures the amount of information a neural pathway transmits about an input signal $x(t)$ is the mutual information rate,

31. ${M}[x,y]={H}[y]-\underbrace{{H}[y|x]}_\text{encoding}= {H}[x]-\underbrace{{H}[x|y]}_\text{decoding}$,

between the stochastic process, $x(t)$, and the stochastic response process, $y(t)$. The entropy rate ${H}$ measures the number of discriminable input or output states, either by themselves, or conditioned on other variables.

The mutual information rates, which is the difference between unconditional and conditional entropy rates, characterises the number of input states that can be distinguished upon observing the output. The response entropy rates ${H}[y]$, for instance, quantifies the number of typical responses per unit time, while ${H}[x|y]$ is a measure of the decoding noise in the model. If this noise is zero, then the mutual information rate is simply ${H}[x]$, provided that this is finite.

The conditional entropy rates ${H}[y|x]$ and ${H}[x|y]$, characterising the noise in the encoding and decoding model respectively, are each greater than zero. In information theory, these quantities are also called equivocation. Hence, both the stimulus and response entropy rates, ${H}[x]$ and ${H}[y]$, are upper bounds for the transmitted information.

Example (Info rate continuous discrete-time Gaussian process):

The OU-process from Eq. (26) is an example of a Gaussian process. Take a discretised version, $\vec x=[x(t_1),x(t_2),...]$, of it such that

$p(\vec x)=|2\pi{\underline K}|^{-1/2}\exp\left(-\frac12(\vec x-\vec\mu)^{\dagger}{\underline K}^{-1}(\vec x-\vec\mu)\right)$

$H_n=\frac12\ln2\pi|{\underline K}| + \frac n2$

This is similar to the quantisation problem. It might be reasonable to drop the $n$ term (sometimes this is done, sometimes not). For the one dimensional case we have

$H=\frac12(1+\ln2\pi\sigma^2)$

or, if we drop the $n=1$

32. $\sigma^2=\frac1{2\pi}e^{2H_\text{Gauss}(x)}$

Any practical consequences?

Def. (Estimation error):

For a random variable $x$ the estimation error of an estmator $\hat x$ is

${\langle}(x-\hat x)^2{\rangle}$

The best estimator is the mean, so the statisticians say. Therefore a lower bound to the estimation error is given by

33. ${\langle}(x-\hat x)^2{\rangle}\geqslant{\langle}(x-{\langle}x{\rangle})^2{\rangle}=\sigma^2 =\frac1{2\pi}e^{2H_\text{Gauss}(x)}\geqslant\frac1{2\pi}e^{2H(x)}$.

The lasst inequality followed from Eq. (29).

Example (Info rate continuous continuous-time Gaussian process):

Up to an additive constant the entropy of a multivariate Gaussian was⁷⁸

$H=\frac12\ln|{\underline K}|=\frac12\tr\ln{\underline K}=\frac12\tr\ln{\underline \Lambda}=\frac12\sum_k\ln\lambda_k$.

First let us observe the process for ever $\vec x=[x(-\infty),...,x(\infty)]$, a bi-infinte series with countable elements. The elements of the covariance matrix ${\underline K}_{ij}=c(i*dt-j*dt)$. The orthogonal eigen-function for the continous covariance operator on $t\in{I\!\!R}$ are the Fourier bases. It can be shown that in the continuum limit

$H=\int df\,\ln\lambda(f)$

The result is due to Kolmogorov see also (???,Golshani and Pasha (2010)).

5.5 Linear stimulus reconstruction and a lower bound on the information rate (decoding view)

Without a complete probablilistic description of the model the mutual information can not be calculated. And even with a model the involved integrals may not be tracktable. At least two strategies to estimate it exist, though: Either, create a statistical ensemble of inputs and outputs by stimulation, followed by (histogram based) estimation techniques for the mutual information; or, find bounds on the information that can be evaluated more easily. In general, the estimation of mutual information from empirical data is difficult, as the sample size should be much larger than the size of the alphabet. Indeed, each element of the alphabet should be sampled multiple times so that the underlying statistical distribution can, in principle, be accurately estimated. But this prerequisite is often violated, so some techniques of estimating the information from data directly rely on extrapolation (???). The problem becomes particularly hairy when the alphabet is continuous or a temporal processes had to be discretised, resulting in large alphabets.

Another approach, which will allow us to perform a theoretical analysis of phase dynamics, relies on a comparison of the neuronal ``channel’’ to the continuous Gaussian channel (???,Cpt.~13) is analytically solvable (Cover and Thomas 2012). The approach can be used to estimate the information transmission of neuronal models (???). Also experimental system have ben analysed in this way, e.g.:

1. the spider’s stretch receptor (???);
2. the electric sense of weakly electric fish (???) and paddle fish (???);
3. the posterior canal afferents in the turtle (???).

It was prooven that in that this method leads to a guaranteed lower bound of the actual information transmitted (???).

If one has experimental control of the stimulus ensemble it can choosen to be a Gaussian process with a flat spectrum up to a cutoff as to not introduce biases for certain frequency bands. The mutual information between stimulus $x(t)$ and response $y(t)$ can be bound from below as

34. ${M}[x,y] = {H}[x] - {H}[x|y] \geqslant {H}[x] - {H}_\text{gauss}[x|y],$

Here, ${H}_\text{gauss}[x|y]$ is the equivocation of a process with the same mean and covariance structure as the original decoding noise, but with Gaussian statistics. The conditional entropy of the stimulus given the response is also called reconstruction noise entropy. It reflects the uncertainty remaining about the stimulus when particular responses have been observed.

It turns out that the inequality in Eq. (33) also holds if the estimator is conditioned. Say from the output of the neuron we estimate its input

$\hat x(t) = \hat x_t[y].$

So if the process has a stationary variance

${\langle}(x(t)-\hat x(t))^2{\rangle}_{x|y}\geqslant\inf\limits_{\hat x}{\langle}(x(t)-\hat x(t))^2{\rangle}_{x|y}={\langle}(x(t)-{\langle}x(t){\rangle}_{x|y})^2{\rangle}_{x|y}=e^{2{H}_\text{gauss}[x|y]}.$

The second line uses the fact that in this case the optimal estimator is given by the conditional mean. We have the following bound on the equivocation

35. ${H}[x|y] \leqslant {H}_\text{gauss}[x|y] \leqslant {{\textstyle\frac12}}\ln{\langle}(x(t)-\hat x(t))^2{\rangle}\leqslant\ln{\langle}n^2(t){\rangle},$

The deviation between stimulus and its estimate, $n(t)=x(t)-\hat x(t)$, is treated as the noise process.

In order to obtain a tight bound the estimator $\hat x(t)$ should be as close to optimal as possible. For the case of additional information given by the response of the neural system $y(t)$ to the process $x(t)$, the estimator should make use of it, $\hat x_t[y]$. For simplicity one can assume it is carried out by a filtering operation, $\hat x(t)= (f*y)(t)$ specified later (Gabbiani and Koch 1998). Like the whole system the noise process is stationary, and its power spectral density, $P_{nn}(\omega)$, is

${H}_\text{gauss}[x|y]\leqslant {{\textstyle\frac12}}\ln{\langle}n^2(t){\rangle}={{\textstyle\frac12}}\int_{-\infty}^\infty\frac{{\mathrm{d}}\omega}{2\pi}\ln P_{nn}(\omega).$

Together

36. ${M}[x,y]\geqslant \frac12\int_{-\infty}^\infty\frac{{\mathrm{d}}\omega}{2\pi}\;\ln\left(\frac{P_{xx}(\omega)}{P_{nn}(\omega)}\right)$

So as to render the inequality in Eq. (35) as tight a bound as possible one should use the optimal reconstruction filter from $y$ to $\hat x$. In other words, it is necessary to extract as much information about $x$ from the spike train $y$ as possible.

The next step should be to find an expression for the noise spectrum, $P_{nn}(\omega)$, based on the idea of ideal reconstruction of the stimulus. As opposed to the forward filter, the reconstruction filter depends on the stimulus statistics (even without effects such as adaptation). We seek the filter $h$ that minimises the variance of the mean square error

37. ${\langle}n^2(t) {\rangle}= {\langle}(x(t) - \hat x(t))^2{\rangle}$, with $\hat x(t)=\int{\mathrm{d}}\tau\,h(\tau) y(t-\tau)$.

Taking the variational derivative of the error w.r.t.
the filter (coefficients) $h(\tau)$ and equating this to zero one obtains the orthogonality condition for the optimal Wiener filter (???)

38. ${\langle}n(t) y(t-\tau){\rangle}= 0$, $\forall\tau$.

Inserting the definition of the error, $n(t)=x(t)-\hat x(t)$, into Eq. (38) yields

${\langle}x(t)y(t-\tau){\rangle}-\int{\mathrm{d}}\tau_1\,h(l) {\langle}r(t-\tau_1)r(t-\tau){\rangle}= R_{xy}(\tau) - (h\ast R_{yy})(\tau) = 0$

In order to obtain $h$ we need to deconvolve the equation, which amounts to a division in the Fourier domain

39. $P_{xy}(\omega)= H(\omega) P_{yy}(\omega)\Longrightarrow H^\text{opt}(\omega) = \frac{P_{xy}(\omega)}{P_{yy}(\omega)}$.

To compute the mutual information rate, we now calculate the full auto-correlation of the noise when the filter is given by Eq. (39). For an arbitrary filter $h(t)$, we have

$R_{nn}(\tau)={\langle}n(t)n(t+\tau){\rangle}={\langle}n(t)x(t+\tau){\rangle}-\int{\mathrm{d}}\tau_1\,h(\tau_1) {\langle}n(t)y(t+\tau-\tau_1){\rangle}$.

When the orthogonality condition of Eq. (38) holds, the right-most term vanishes. Proceeding by expanding the first term algebraically leads to an expression for the noise correlations

$R_{nn}(\tau)={\langle}n(t)x(t+\tau){\rangle}=R_{xx}(\tau)-\int{\mathrm{d}}\tau_1\,h(\tau_1)R_{xy}(\tau-\tau_1)$.

This expression can be Fourier transformed in order to obtain the required noise spectrum

$P_{nn}(\omega)=P_{xx}(\omega)-H(\omega)P_{xy}(\omega)=P_{xx}(\omega)-\frac{|P_{xy}(\omega)|^2}{P_{yy}(\omega)}$,

where the definition of the optimal filter, Eq. (39), was utilised. This result can then be inserted into Eq. (36) to obtain the following well known bound on the information rate (???,lindner2005j:mi,holden1976b,stein1972j:coherenceInfo)

40. $\mathcal M[x,y]\geqslant-\frac12\int_{-\omega_c}^{\omega_c}\frac{{\mathrm{d}}\omega}{2\pi}\ln\left(1-\frac{|P_{xy}(\omega)|^2}{P_{xx}(\omega)P_{yy}(\omega)}\right)$.

This information bound involves only spectra and cross-spectra of the communication channel’s input and output processes which are experimentally measurable in macroscopic recordings . The channel, in this case the neuron, can remain a black box. But since we can bridge the divide between microscopic, biophysical models and their filtering properties, we will, in the following section, derive the mutual information rates.

Def. (spectral coherence):

The expression in Eq. (40) is termed the squared signal response coherence

41. ${c}^2(\omega)=\frac{|P_{xy}(\omega)|^2}{P_{xx}(\omega)P_{yy}(\omega)}$.

It quantifies the linearity of the relation between $x$ and $y$ in a way that it equals 1 if there is no noise and a linear filter transforms input to output. Both nonlinearities and noise reduce the coherence. The coherence can be estimated from data using the FFT algorithm and spectrum estimation. It is implemented in the free software packages scipy and matplotlib.

What renders the coherence a useful quanity? While the cross-spectrum informs us when stimulus and output have correlated power in a spectral band, the normlisation with the output auto-spectrum can be crucial. Say we find a particular high power in $P_{xy}(\omega)$, this may not be ralted to the stimulus but could just be the intrinsic frequency of the neuron itself.

The coherence is a quantity without mention of the explicite decoding filter, in fact it is symmetric in input and output just as the mutual information. This is beneficial because one can now take the encoding view in the next chapter.

6 Linear response filter

The stimulus spectral density is given by the environment or controlled by the experimental setup, while cross- and output spectra need to be measured or calculated from the model in question. In this lecture cross-spectral and spike train spectral densities are derived from phase oscillator, see Eq. (17), that are in turn derived from biophysical models. This means we do not treat the channel as a blackbox but assume a particular model.

The first quantity we need to calculate Eq. (41) is the cross-spactrum. On the one hand it is the Fourier of the cross-corrlation, on the other it can be written as averages of the Fourier transforms (FT and average are linear operation).

42. $P_{yx}(\omega)={\langle}{\langle}{{\tilde y}}(\omega){{\tilde x}}^*(\omega){\rangle}_{y|x}{\rangle}_x= {\langle}{\langle}{{\tilde y}}(\omega){\rangle}_{y|x}{{\tilde x}}^*(\omega){\rangle}_x$.

What has happened here? The cross-spectrum can be obtained by averaging the Fourier transformed quantities over trials and the stimulus ensemble. The average can be split into the conditinal average over trials ${\langle}\cdot{\rangle}_{y|x}$, given a fixed, frozen stimulus and the average over the stimulus ensemble, (${\langle}\cdot{\rangle}_x$). The former is essential an average over the encoding noise (Chacron, Lindner, and Longtin 2004,Lindner, Chacron, and Longtin (2005)).
Observe that ${\langle}\tilde y(\omega){\rangle}_{y|x}$ is Fourier transform of the trial averaged firing rate conditional on a frozen stimulus