1 Preface

1.1 Content Summary

Nerve tissue is comprised of excitable cells, supporting cells and a surrounding bath. This seminar lecture reviews the biophysical basics of excitability, including the stochastic processes that govern voltage-gated ion-channels in nerve membranes. The composition of ion channels, among other factors, determines the functional properties of action potential pulses, which are analysed using dynamical systems theory. The statistics of pulse trains is derived from the microscopic level of ion-channel fluctuations. This will also serve the understanding of a self-consistent description of macroscopic neural-network dynamics.

1.2 (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 pressure to improve it. Please write to the lecturer if you find fault of any kind with it.

Although written in a somewhat formal fashion the book is by no means a text in applies mathematics. It completely lacks the mathematical rigor to meet such standards. The choice of wrapping statements in definitions, observations and propositions arose from the need to be able to easily refer to them and to render the logical buildings blocks of some arguments easier to identify and follow through the text. Let’s see how that worked out …

2 Overview

There are many different lectures on “the biophysics of computation”, or “computational neuroscience”. The topic by now is so vast, that only a small non-exhaustive subset of the subject can be discussed within the temporal confinement of a lecture. This one is no different.

The idea here is to approach questions like “What is a nerve impulse?” or “Are there different types of excitability?” from a mathematical and a biophysical angle.

From the biophysical point of view, one has to addresse the many buidling blocks of nervous systems: There is needs to understand the biophysics of lipid sheets and the chemical gradiente across them. Voltage-gated ion channels swimming in the membrane should be understood, including their stochastic dynamics. What is going on in the extracellular medium might be a topic. And of course what are the consequences of wiring neurons in networks via synapses or gap-junctions. Also the enourmous energetic demands of nerve tissue are important factors in understanding its structure and function.

Many of these components of nervous systems derive from preadaptations, meaning precurses have been arround long before the first nervous system emerged some 600 Myr ago. This is certainly true for excitability, which is found in singe celled Paramecium and even bacteria. Synaptic components like the SNAR complex is present in fungi and some of the transmitters are used in bacterial communication.

Naturally, there is a risc of getting overwelmed by these beautiful details. To explain the complexity of nervous systems ab initio seems a mindbogglingly difficult task. Yet, there are already all-atom simulation of ion channels in the membrane. Still, before and during discussing the structural aspects it is insightful to ask what are the exigences of a nervous system. Not in any teleological sense, theoretical biology is over that, but from an evolutionary standpoint. One may even formulate certain desiderata, without claiming that brains were design by anyone.

2.1 Levels of explanation

If by attending the lecture or goint through its notes the amount of confusion becomes unbearable Tinbergen’s levels of analysis might help

While, categories and schemes like the above are of utility in structuring ones thoughts about nature, they should not be taken as dogmatic fences that do not allow you to think, speculate outside of them. Nature is to beautiful to get involved into dogmatic fights over reified wordings, leave them to the humanities.

Roughly one might associated certain disciplins of applied math that can be utilised in the four categories from above.

3 The Floquet theory of the action potential

3.1 Periodic event

Of many natural systems we observe only some drastically obvious events, while the underlying complex evolution of all state variables is hidden. The malatonin falling level below threshold (in the suprachiasmatic nucleus) causes a “waking-up event”. Results the complex regulatory network causing gap phase 1 $\to$ synthsis $\to$ gap phase 2 $\to$ mitosis transition result in clearly visible cell division events. Sadly, extinction events can be the results of changes in complex foodwebs, the dynamics of which is not completely determined. Action potentials are no different. The dynamics involves charging the cell membrane and the nonlinear kinetics of several ion channels, but the electric impulse is an clearly detectable event (with a glass electrode at least).

3.2 Pulse-based communication

Why do nervous system use pulses, called action potentials or spikes to communicate? In fact, not all do. The information picked up by your photo receptors is processed by cells in the retina (amacrine cells, horizontal cells) still as graded potentials, which are then recoded into a impulse based representation. Conversely, the visusal system of flys starts out with pulses representation and converts it to graded potentials upstream. Last but not least, the small nervous system of nematode C. elegans (302 neurons) seemed to operate entirely on graded potentials. No spikes had been reportet for a long time. But in 2018 someone found Ca$^{2+}$ based spikes in one of their neurons.

One of the differences between passive graded potential and an active action potential is that the latter does not suffer from the same kind of siganl degradation over long distances. There could be active more or less graded potentials, but these are less temporally precise as a short pulse and it turns out that many computations rely on on exactly this temporal precession.

A neurobiological “dogma” states that action potentials follow the all-or-nothing principle (E.D. Adrian and Zotterman 1926). Fortunately, the jelly fish Aglanta digitalae does not care much for dogmas and encodes its swimming patterns in different action potential shapes. So its not really a dogma or even a principle. But most action potentials really are very stereotypical events, hence its a common property. The all-or-nothing property can be construed to mean that the exact shape of the action potential does not matter¹ and that information is instead stored in the exact spike times or the frequency of their occurence.

This can be formalised to mean that the voltage evolution ought to be represented by a shifted pulse wave form

1. $v(t) = \sum_k A(t-t_k) = A(t) * \sum_k\delta(t-t_k)$,

where $t_k$ are the spike times, which need to be defined later.

It also means that the action potential $A(t)$ should be stable under perturbations, but that inputs ought to be able to shift the occurrence of spikes in time. In mathematical terms, these desiderata read

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.

3.3 Tonic spikes are limit cycles

The mathematical object that suffices the desiderata (i) and (ii) are limit cycles. They are closed orbits in phase space. Naturally, the state space need to be at least two dimensional for it to contain a limit cycle. What are the state variables? The membrane voltage, $v(t)$, typically in milivolts sounds like a reasonable one, at least that is what electrophysiologist measure. The conformational state of the ion channels also changes with time, $c_k(t)$, and so do ion concentration, say $[\mathrm{K}^+]_o$. The exact governing equations are going to be derived in the following chapters. For now, denote a vector of state variables as

$\vec x=\begin{pmatrix}v\\c_1\\c_2\\ [\mathrm{K}^+]_o\\ \vdots\end{pmatrix}.$

A general differential equation could have the following structure

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

$\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\alpha$ are system’s parameter. Parameters could be cell size and other fixed properties, but also a constant applied bias current. The change of solutions in response to variations in system’s parameter will be studied later.

In the absence of perturbations, $\vec\eta=0$, Eq. (2) becomes a homogeneous ordinary differential equation

3. $\dot{\vec x} = \vec F(\vec x,\vec\alpha).$

We are interested in the cases in which, at least for cerctain parameter values $\vec\alpha$, $\exists$ a ${P}$-periodic limit cycle solution, with $\vec x_{\mathrm{lc}}(t)=\vec x_{\mathrm{lc}}(t+{P})$, also known as periodic orbit. This solution is identified as the action potential if it sufices the requirements of amplitude-stability and time-shiftability (i) and (ii).

3.4 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

4. $\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. (4) into Eq. (3) 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}})}_{{\boldsymbol{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 ${\boldsymbol{J}}(t)=\nabla\vec F(\vec x_{\mathrm{lc}}(t))$. Note that since the limit cycle solution is ${P}$-periodic, so is $\boldsymbol{J}(t)=\boldsymbol{J}(t+{P})$.

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

5. $\dot{\vec y}(t)={\boldsymbol{J}}(t)\vec y(t).$

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

6. $\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}}}={\boldsymbol{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. (5) can be written in the form of a ${P}$-periodic similarity matrix² and an matrix exponential

7. $y(t)=\boldsymbol{R}(t)e^{t\boldsymbol{\Lambda}}.$

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

Recall the matrix exponential

Def: Matrix Exp

Let $\boldsymbol{A}\in{I\!\!\!\!C}^{n\times n}$ then $\exp\boldsymbol{A}=\sum_{k=0}^\infty\frac1{k!}\boldsymbol{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 $\boldsymbol{A}$, i.e., $\boldsymbol{A}\vec w_i=\lambda_i\vec w_i$ then by using this identity $k$-times

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

Inserting the Floquet Ansatz into Eq. (5) yields

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

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

9. $\dot{\boldsymbol{R}}=\boldsymbol{J}(t)\boldsymbol{R}(t)-\boldsymbol{R}(t)\boldsymbol{\Lambda}.$

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

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

3.4.1 Eigensystem of the Floquet matrix

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

11. $\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}$.

The $\mu_k$ are called Floquet exponents, for reasons becomming obvious below. One can now define a “rotated” eigensystem

Def. (“Rotated” eigensystem and Floquet modes) #def:floqmode

The rotated eigensystem of a spike is $\vec W_i(t)=\boldsymbol{R}(t)\vec w_i$ and $\vec Z_i(t)=\vec z_i\boldsymbol{R}^{-1}(t)$. For which ipso facto obeys

12. $\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. (9) and (10) and use this definitions we get

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

and

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

As a final test of what has been chiefed by introducing the new time-dependent coordinate system of the right and left Floquet modes, project the general solution $y(t)$ from Eq. (7) on the the adjoint Floquet modes, $Z_k(t)$:

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

Along these rotating basis one observes a simple exponential contraction or expansion depending on the sign of the Floquet exponents. The maginitude of the Floquet exponents is the rate of convergence in units of 1/Time. What if a Floquet exponent is zero? In fact, by comparing Eq. (6) and (13) one may observe that the Goldtone mode is alwasys a solution to the

Note that if $\mu_0=0$,

15. $\vec W_0(t)=\frac{\mathrm{d}}{{\mathrm{d}}t}\vec x_{\mathrm{lc}}(t)$

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

16. $\vec Z_0(\phi)\cdot\frac{\mathrm{d}}{{\mathrm{d}}t}\vec x_{\mathrm{lc}}(\phi)= \vec Z_0(\phi)\cdot\vec F(\vec x_{\mathrm{lc}}(\phi))=1$

W.l.o.g, let the Floquet exponents of a stable limit cycle be ordert such that

$0=\mu_0\geqslant\mu_1\geqslant...\geqslant\mu_{n-1}$

What are left Floquet modes useful for? Here is an example. First, rewrite Eq. (14) for $k=0$ in the form

$(\frac{{\mathrm{d}}}{{\mathrm{d}}t}+J^\dagger(t))\vec Z_0=0$

Say one is interested in the change of the period of a limit cycle oscillator w.r.t. a parameter change. By rescaling time

$\frac{d\vec x_{\mathrm{lc}}}{dt}-{P}(\alpha)\vec F(\vec x_{\mathrm{lc}},\alpha)=0$

Taking the total derivative (and chain rule) one finds

$\frac{\partial}{\partial\alpha}\frac{d\vec x_{\mathrm{lc}}}{dt}-\frac{\partial {P}}{\partial\alpha}\vec F-{P}\frac{\partial \vec F}{\partial\alpha}-{P}\nabla\vec F\frac{\partial\vec x_{\mathrm{lc}}}{\partial\alpha}=0$

The derivatives of the first term can be excanged. Now project the equation from the left onto the zeroth Floquet mode

$-\frac{\partial {P}}{\partial\alpha}(\vec Z_0\cdot\vec F)-{P}(\vec Z_0\cdot\frac{\partial \vec F}{\partial\alpha}) +\vec Z_0\cdot(\frac{d}{dt}-{P}J(t))\frac{\partial\vec x_{\mathrm{lc}}}{\partial\alpha}=0$

The last term can be rewritten by transposing the operator

$\vec Z_0\cdot(\frac{d}{dt}-{P}J(t))\frac{\partial\vec x_{\mathrm{lc}}}{\partial\alpha}= -\underbrace{((\frac{d}{dt}+{P}J^\dagger(t))\vec Z_0)}_{=0}\cdot\frac{\partial\vec x_{\mathrm{lc}}}{\partial\alpha}$

Hence

17. $\frac{\partial {P}}{\partial\alpha}=-{P}(\vec Z_0\cdot\frac{\partial \vec F}{\partial\alpha})$

4 Spikes event time series and isochrons

4.1 Spike times and Isochrons

There is still a lack of a formal definnition of the spike time $t^{\text{sp}}$ in Eq. (1). Electrophysiologists typically envision a voltage threshold, but these would only exist if the voltage dynamics was a one dimensional equation, which it is not. In fact, it can not be if true limit cycles are to exist. In the general state-space of action potential dynamics there may exist a generalised threshold manifold. A reasonable alternative definition would take the (local) peak/maximum of the voltage

$v(t^{\text{sp}})=v_\mathrm{max}$.

Again, a definition based on properties of the state variables and it also leaves all other dimensions undefined. The general problem seems to be to uniquely associate time with states. If the dynamics was restricted to only on the limit cycle this is easy. As the maximum voltage value is only reached once per spike it also implies unique values for the other state variables, wich can be called the spike boundary condition: $\vec x_\mathrm{sp}$. Most importantly it also implies a unique time indices at which such events occur, $\vec x_{\mathrm{lc}}(t^{\text{sp}})=\vec x_\mathrm{sp}$. The normalised time index of points on the limit cycle is also called the phase, $\phi\in[0,1)$ s.t. $\vec x_{\mathrm{lc}}({P}\phi)=\vec x_{\mathrm{lc}}({P}(\phi+1))$. Given synaptic bombardment, intrinsic noise and other perturbations the dynamics, however, may be temporarily not exactly on the limit cycle. The simple idea is that the time phase of any point in the basin of attraction of a limit cycle can be defined as the phase it converges back to when asymptotically approaching the limit cycle. To make this precise, first define the flow:

Def. (Flow of an ODE) #def:flow

For a given initial condition $\vec x_0\in B\subset{I\!\!R}^n$, where $B$ is the basin of attraction of the limit cycle, the flow or trajectory induced by the ODE in Eq. (3) is denoted as $\vec x(t;x_0)$.

Numerically, the flow can be evaluated by forward-integration of the ODE starting at a given initial value.

Def. (Assymptotic phase) #def:assphase

For point $\vec y\in{B}$, in the basin of attraction of the limit cycle solution of Eq. (3), the asymptotic phase, $\psi$, is defined as

18. $\psi(\vec y)=\argmin_\psi\lim_{t\to\infty}\|\vec x(t-\psi {P};\vec y)-\vec x_{\mathrm{lc}}(t;\vec x_\mathrm{sp})\|$.

Note that the limit cycle solution was initiated at time $t=0$ at the spike boundary condition. Based on the asymptotic phase, the concept of isochrons can be introduced (Winfree 1974; Osinga and Moehlis 2010) as the level sets of constant phase:

Def. (Isochrons) #def:isochron

The isochron manifolds, $\mathcal I(\phi)$, are the level sets of the asymptotic phase variable

$\mathcal I(\phi)=\{\vec y\in{B}:\psi(\vec y)=\phi\}$.

An isochron defines a $n-1$ dimensional manifold of points that converge to the same asymptotic phase on the limit cycle.

Based on the isochron one can now define a spike as a crossing of the flow/trajectory with $\mathcal I(0)$.

Note: There are several other ways to define a phase (Hilbert transform, linear interpolation between spikes, …). In principle, there should be monotonous mappings between each other.

4.2 Neutral dimension and phase shifts

An alternative to studying the high dimensional set of coupled nonlinear biophysical equations is to find an (approximately) input/output(I/O)-equivalent equation for the phase evolution like in Fig. #. In the figure the phase is ploted as modulo 1. However, it is also convenient to consider the freely evolving phase since this is related to the spike count process,

$N(t)=\lfloor\phi(t)\rfloor$.

The $k^\mathrm{th}$ spike time is then defined as

$\phi(t^{\text{sp}}_k)=k$.

A spike-train can be written as

19. $r(t)=\sum_k\delta(t-t^{\text{sp}}_k)=\sum_k\delta(\phi(t)-k)$.

Def. (I/O-equivalent phase oscillator) #def:io_phase

An I/O-equivalent phase oscillator produces with its threshold crossing process the same spike times $\{t^{\text{sp}}_k\}_k$ as a particular conductance-based model (Lazar 2010).

Fig. #fig:phasered

The phase equation turns the spiking neuron into a one-dimensional level-crossing process with equivalent spike times. The phase is plotted modulo 1.

The evolution of the phase of Eq. (2) is given by the chain rule as

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

20. $\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)$.

The geometric interpretation of the phase gradient, $\nabla\phi$, is that it points perpendincular to the isochrons towards the steapest change of the asymptotic phase variable $\phi$.

4.3 The empirical infinitesimal phase response curve

Def. (empirical PRC) #def:empiricalPRC

Given a short $\delta$-pulse of amplitude $\epsilon$, delivered at a certain phase, $\phi$, (time point) in the inter-spike interval of a tonically firing neuron the linear phase response is defined as

$\prc(\psi)=\lim\limits_{\epsilon\to0}\frac1\epsilon\frac{{P}_0-{P}_\epsilon(\psi)}{{P}_0}$.

The underlying assumption of interpreting PRC is that it is a function, $\prc(\psi):[0,1)\mapsto{I\!\!R}$, that mapps phase or time of occurence of a perturbation to a persistent measurable phase shift. (Please appreciate how idealised and “noise free” these definitions are.) The convention is that a prolongation of intervall corresponds to a phase delay, while shortinging means phase advance.

The perturbed period can be written as the original period minus the phase shift converted back to time (by multiplying with the unperturbed period)

${P}_\epsilon={P}_0-\Delta\phi{P}_0$.

Given that the time evolution of $\phi(t)$ has been perturbed at some time point $\psi {P}_0<{P}_0$, The phase shift can be measured at the spike time of the unperturbed system as

$\Delta\phi = \phi({P}_0)-1$.

The empirical PRC in Def. # can then be expressed as

$\prc(\psi)=\frac1\epsilon\frac{{P}_0}{{P}_0}(1-(\phi({P}_0)-1)-1)=\frac1\epsilon(\phi({P}_0)-1)$.

Apply the identify operator of functions (convolution with a $\delta$-function) on the left and the fundamental theorem of calculus on the right hand side then to first order

$\epsilon{P}_0\int_0^{{P}_0}\prc(\phi(t))\delta(\psi-\phi(t)){\mathrm{d}}t=\int_0^{{P}_0}(\dot\phi-\frac1{{P}_0}){\mathrm{d}}t$

This step is not really allowed, but for well behaved integrals we can identify the integrands to see

$\dot\phi=\frac1{{P}_0}+\epsilon\prc(\phi)\delta(\psi-\phi)$.

One can bring this equation in agreement with the phase evolution in Eq. (20), if one identifies

$\nabla\phi(\psi)=\prc(\psi)$ and $\nabla\phi(\psi)\cdot\vec F(\psi)=\frac1{{P}_0}={f}_0$.

4.4 Constant phase shift and adjoint equation

The assumption underlying the experimental proceduce described in Def. # is that the phase shift stays constant in time, say $\Delta\phi$. Without perturbation there is not long-term shift of the phase. Hence it must be the perturbation off the limit cycle which causes the phase shift. The phase change must be due to the $\Delta\vec x(t)$ and is given by

$\Delta\phi=\nabla\phi\cdot\Delta\vec x(t)=$ const. From $\frac{{\mathrm{d}}\Delta\phi}{{\mathrm{d}}t}=0$ one gets

$0=\frac{\mathrm{d}}{{\mathrm{d}}t}\nabla\phi\cdot\Delta\vec x+ \nabla\phi\cdot\frac{\mathrm{d}}{{\mathrm{d}}t}\Delta\vec x$

Again the evolution of $\Delta x(t)$ can be linearised if it is small

$0=\frac{\mathrm{d}}{{\mathrm{d}}t}\nabla\phi\cdot\Delta\vec x+ \nabla\phi\cdot {\boldsymbol{J}}(t)\Delta\vec x= (\frac{\mathrm{d}}{{\mathrm{d}}t}\nabla\phi+{\boldsymbol{J}}^\dagger(t)\nabla\phi)\cdot\Delta\vec x$

Such that by Fredholm’s alternative

$\frac{\mathrm{d}}{{\mathrm{d}}t}\nabla\phi=-{\boldsymbol{J}}^\dagger(t)\nabla\phi$

This is the same ODE as for the zeroth left Floquet mode, rendering them identical up to scaling.

5 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, this was a chived by doing a diffusion approximation in the thermodynamic limit of many ion channels.

Finally, in 1998 the structure of a K⁺ ion channel from Streptomyces lividans, which is nearly identical to Drosophila’s shaker channel was discovered using X-ray analysis (Doyle 1998).

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)).

5.1 Empirical assertions about ion channels

Observation (Conserning the ion channels of excitable membranes) #obs:chan

1. Ion pores show different degrees of selectivity to ionic species. Understanding of the mechanisms of the ionic filter was achieved by structural analaysis (Doyle 1998). Potassiume channels conduct K⁺ 1000 fold more than Na⁺ and allow for conduction rate of 10⁸ ions per second.

2. In single channel patch experiments ion channels how discrete levels of conductivity, which look like a Telegraph process.

3. Independent ion channels (no cooperativity or spatial heterogeneity)

4. Membrane bound proteins like many others change conformation on various triggering signals:

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

To relate measured voltage changes to ion flux accorss the membrane note that according to Kirchofs law and electro neutrality

$C\frac{{\mathrm{d}}v}{{\mathrm{d}}t}=\frac{F\Omega}{A} \Big(\frac{{\mathrm{d}}[{\text{Na}^+}]_i}{{\mathrm{d}}t}+\frac{{\mathrm{d}}[{\text{K}^+}]_i}{{\mathrm{d}}t}+\frac{{\mathrm{d}}[{\text{Cl}^-}]_i}{{\mathrm{d}}t}\Big)$

The currents are assumed to follow Ohms law:

$\frac{F\Omega}{A}\frac{{\mathrm{d}}[{\text{K}^+}]_i}{{\mathrm{d}}t}=\gamma N_K(E_K-v)$

where $\gamma$ is the unitary conductance and $N$ is the total number of open pores.

5.2 First-order rate kinetics

Starting with a hypothetical ion pore of which there is a fixed total concentration, [Pore]=const. The pore an open conformation, $O$, in which, say only K⁺, ions can pass and a cloased states, $C$, which blocks ion flow completely:

Def. (Two-state ion-channel) #def:twochan

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

Based on the law of mass-action, the corresponding first-order rate equation is

$\frac{{\mathrm{d}}[O]}{{\mathrm{d}}t}=\alpha[C]-\beta[O]=\alpha([\text{Pore}]-[O])-\beta[O]$

The second equality is by conservation of mass $[O]+[C]=[\text{Pore}]=$ const.

has the typical strucutre of its stochastic counter parts the Master equation of a birth-death process.

Hodgkin and Huxley normalised the equation by defining $n(t) = \frac{[O]}{[\text{Pore}]}$.

$\frac{{\mathrm{d}}n}{{\mathrm{d}}t}=\alpha(1-n) -\beta n=\left(\frac\alpha{\alpha+\beta}-n\right)(\alpha+\beta) =(n_\infty-n)/\tau$

Under a voltage-clamp experiment, for each voltage level this equation predicts an exponential charge up to a level $n_\infty$ with the time constant $\tau$ that could be fitted.

Hodgkin and Huxley discovered that their measured traces could not be fitted well by the above equation.

Def. (A tetrameric K⁺ pore) #def:fourchan

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

5.3 The ion channel as a Markov model

Instead of writing down heuristic macroscopic laws for channel concentrations, one can attempt to derive them from microscopic dynamics of single pores. For this the state diagram in Def. # is reinterpreted as the state of a single molecule.

Based on Obs. # such single channel dynamics can be described by a finite state Markov model which is mathematically described as a Master equation.

A Markov model represented by a graph like Def. # is a special kind of chemical reaction network studied in systems biology. The connectivity/topology can be stored in a

Def. (Stoichiometry matrix) #def:stoichmat

The stoichiometry matrix, $\boldsymbol{N}\in{I\!\!N}^{m\times n}$ describes $n$ reactions (transitions) involving $m$ chemical species (states). It can be decomposed into the educt matrix $\boldsymbol{E}$ and the product matrix $\boldsymbol{P}$ by $\boldsymbol{N} = \boldsymbol{P} - \boldsymbol{E}$

The stoichiometry matrix of Def. # would look like

$\boldsymbol{N}=\begin{pmatrix} -1 & 0 & 0 & 0 & 1 & 0 & 0 & 0 \\ 1 & -1 & 0 & 0 & -1 & 1 & 0 & 0 \\ 0 & 1 & -1 & 0 & 0 & -1 & 1 & 0 \\ 0 & 0 & 1 & -1 & 0 & 0 & -1 & 1 \\ 0 & 0 & 0 & 1 & 0 & 0 & 0 & -1 \end{pmatrix}$

In addition one needs the rate vector

$\vec v = [4\alpha,3\alpha,2\alpha,\alpha,\beta,2\beta,3\beta,4\beta]^{\mathsf{T}}$

to define the complete reaction network for the channel.

To understand the dynamical evolution of occupation propability of different conformational states, take the two-state channel as an example. Define a vector of the probability of being closed and open, $\vec p=[p_C,p_O]^\dagger$, respectively.

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. One needs an update rule of how to get from the probability of being open at time $t$ to the probability of begin open at time $t+dt$. Given are the rates, so the propabilities of change in an interval of time $dt$ are the rates times the interval length, e.g. $dt\,\beta$. For the two state channel this cann be summarised in matrix notation

$\begin{pmatrix}p_C(t+dt)\\p_O(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(\boldsymbol{I}+\boldsymbol{Q} dt\right)\vec p(t)$

The infinitesimal limit yields a kind of Kolmogorov forward-equation

21. $\frac{\mathrm{d}}{{\mathrm{d}}t}\vec p=\boldsymbol{Q}\vec p$.

Def. (Transition rate matrix) #def:transmat

The transition rate matrix (or infinitesimal generator) is $\boldsymbol{Q}\in{I\!\!R}^{m\times m}$.

In general the transition rate matrix can be obtained from the graph

$\boldsymbol{Q} = \boldsymbol{N}\diag(\vec v)\boldsymbol{E}^{\mathsf{T}}$.

Exercise (Transition rate matrix) #ex:mc1

Derive $\boldsymbol{Q}$ of a three state channel

$\boxed{0} \underset{\beta}{\overset{3\alpha}{\rightleftharpoons}} \boxed{1} \underset{2\beta}{\overset{2\alpha}{\rightleftharpoons}} \boxed{2} \underset{3\beta}{\overset{\alpha}{\rightleftharpoons}} \boxed{3}$

Note that Eq. (21), though a consistent description of the occupation probabilities, does not give us realisations or sample paths comparable to the traces measured in single-channel patch clamp experiments.

5.4 Sample path of a telegraph process (simulation a jump process)

Given a total number of ion channels $N$, the goal is to track the actual number of channels in a particular state. Based on reaction rate theory and given a certain reaction rate $a$, the number $k$ of reaction events in a time window is

$k\sim\text{Poisson}(a)=\frac{a^ke^{-a}}{k!}.$

If several reactions can occur and several molecules are in volved then (because Poisson distribution is closed under convolution) one has

$k\sim\text{Poisson}(\underbrace{\sum_jN_ja_j}_\lambda)$.

So the rate of the occurence of any reaction is given by this linear combination of the rates of individual reactions. The waiting time distribution between reactions is exponential with the same rate:

$p(t)=\lambda e^{-\lambda t}$

Starting in a particular state $\vec N=[N_1,...,N_n]$ at time $t_0$ the life time of staying in that state until $t_0+t$ is

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_ja_j}{\sum_{k=1}^{N_{\mathrm{reac}}}N_ka_k}=N_ja_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

5.5 Statistical properties

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$

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

5.5.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

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

For the K⁺ channel with $n=4$ conformations states, of which one is conducting, let us further define

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

If there are several conducting states (with possibly different conductance level, one could have e.g., $G=\delta_{2K(t)}+2\delta_{4K(t)}$.

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

22. $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

23. $\frac{\mathrm{d}}{{\mathrm{d}}t}\vec p=\boldsymbol{Q}\vec p$

With formal solution

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

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

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

If Eq. (23) 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:\boldsymbol{Q}\vec v_1=\nu_0\vec v_1=0$. Therefore, the solution can be written as

25. $\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. (22) 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

26. $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.

5.6 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. (22)? Let us try

27. $\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. ((???)) 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. (26). 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

6 Continuation of fixpoints and orbits

Often small changes, say in an ion channel density, has only small quantitative effect on a neurons behaviour.

Often numerical solutions to nonlinear ordinary differential equations are found by (forward) time-integration. An interesting alternative is to track a found solution through parameter space, for which the solution must be persistent in a certain parameter. If it is not, then a bifurcation occurs and one observes a qualitative change of the solution. The study of system solutions under parameter perturbations is also the subject of sensitivity analysis.

For book on bifurcation theory consult (Izhikevich 2007) and numerical bifurcation analysis (Kielhöfer 2011).

6.1 Continuation of fixed points

Asume for

$\dot{\vec x}=\vec F(\vec x,\alpha)$

there is a steady state solution

28. $\vec F(\vec x, \alpha)=0$

In biological systems like cells, parameters are rarely stable over longer periods, because there is a continuous turnover (synthesis and degradation) of all molecules. Hence, so one should be interested in families of solutions for varying prarameters, $\vec x(\alpha)$.

Def. (solution banch) #def:bifbranch

A branch, familiy or curve of solutions with regard to an index parameter $\alpha$ is denoted as $\vec x(\alpha)$.

The existence of the solution as a function of the parameter is governed by the implicite function theorem:

Theorem (Implicite function theorem) #thm:implicitefoo

Consider a system of equations

$\vec F(\vec x,\alpha)=0$, with $\vec F\in I\!\!R^n$, $\vec x\in I\!\!R^n$, $\alpha\in I\!\!R$ and $\nabla_{x,\alpha}\vec F\in I\!\!R^{n\times n+1}$.

Let $f$ and $\nabla_{x,\alpha}\vec F$ be smooth near $x$. Then if the Jacobian $\nabla_{x,\alpha}\vec F$ is nonsingular, $\exists$ a unique, smooth solution family $\vec x(\alpha)$ such that

$\vec F(\vec x(\alpha),\alpha)=0$.

This establishes the existence of lines in a bifurcation diagram.
The method of continuation is a predictor-corrector method. In practice, assume the fixpoint is known for one particular parameter value $\alpha_0$, then for a small change in parameter the solution is is predicted

Predictor step:

Taylor’s linearisation

29. $\vec x(p+\delta p)\approx\vec x(p)+\delta p\frac{\partial\vec x}{\partial p}$.

To predict the new solution, the old solution is required and the derivative of the solution w.r.t. the parameter that is changed. How to compute the latter? Take the total derivative of Eq. (28) w.r.t. the parameter

$(\nabla_xf)\frac{\partial\vec x}{\partial p}+\frac{\partial\vec f}{\partial p}=0$.

with formal solution

$\frac{\partial\vec x}{\partial p}=-(\nabla_x\vec f)^{-1} \frac{\partial\vec f}{\partial p}$.

If $\nabla_x\vec f$ is full rank one can use some efficient linear algebra library to find the vector $\frac{\partial\vec x}{\partial p}$ and back insert it into Eq. (29). For too large $\delta p$ the predicted solution will be wrong. Yet, it is a good initial guess form which to find the correct version.

Corrector step:

Newton iterations to find the root of $\vec f(\vec x,p)$

30. $\vec x_{n+1}=\vec x_n-(\nabla_x\vec f)^{-1}\vec f(\vec x_n,p)$

Actually that Newton’s iterations are also obtained by linearisation

$0=\vec f(\vec x_{n+1},p)=\vec f(x_n)+(\nabla_x\vec f)(\vec x_{n+1}-x_{n})$,

which if solved for $\vec x_{n+1}$ yields Eq. (30).

Convergence analysis of Newton iterations yields that with each newton iterations the number of correct decimal places doubles. Hence often a low number of iterations (3-7) suffice to achieve sufficient numerical accuracy.

Example (Continuation of the QIF model) #ex:contQIF

$\dot v = I + v^2 = F(v,I)$ with the solution branches $v(I)=\pm\sqrt{-I}$.

From $\nabla_vF=2v$ and $\frac{\partial F}{\partial I}=1$, one gets $\frac{\partial v}{\partial I}=-\frac{1}{2v}=-\frac{1}{2\sqrt{-I}}$. What happens at $I=v=0$?

Example:

Say $f(x,p)=\sqrt{x^2+p^2}-1$, then the solution branches are $x(p)=\pm\sqrt{1-p^2}$. The linear stability analysis $\partial_xf=\frac x{\sqrt{x^2+p^2}}$ that

$x>0\to\partial_xf>0\to$ unstable
$x<0\to\partial_xf<0\to$ stable

Also $\partial_pf=\frac p{\sqrt{x^2+p^2}}$ and thus $\partial_px=\frac px=\frac p{\pm\sqrt{1-p^2}}$. What happens at $p=\pm1$ or $x=0$?

Or more general: What happens if we find that the condition number of $\nabla\vec f$ explodes?

In the example above the branch switches its stability and it bends back in a fold or turning point. In general folds can occure with and without stability change.

6.2 Local bifurcations: What happens if $\nabla\vec f$ is singular?

Bifurcatin analysis is one big “case discrimination” task. Luckily, and due to topology the cases are, given certain constraints, finite.

6.2.1 Folds and one-dimensional nullspaces

Recall some difinitions

Def. (nullspace, kernel):

The kernel or nullspace of a matrix $\boldsymbol{J}$ is

$N(\boldsymbol{J})=\{\vec x\in I\!\!R^n|\boldsymbol{J}\vec x=0\}$

Def. (range, image):

The range of a matrix $\boldsymbol{J}$ is

$R(\boldsymbol{J})=\{\vec y|\exists x\in I\!\!R^n:\boldsymbol{J}\vec x=\vec y\}$

Def:

Let $Q$ be the projection onto the range

Def:

Eigensystem at the fold. Let Jacobian matrix $\boldsymbol{J}=\nabla\vec f(\vec x(0),p(0))$

$\boldsymbol{J}\vec r_k=\lambda_k\vec r_k$ and $\vec l_k\boldsymbol{J}=\lambda_k\vec l_k$.

So that the nullspace is spanned by $\vec l_0$.

This section considers $\dim N(\boldsymbol{J})=1$. Hence, the implicite function theorem is not applicable. From the example above it is apparent that $\vec x'(p)=\infty$ at the bifurcation. The problem can be circumvented by defining a new “arclength parameter”, $p(s)$. The bifurcation branch is then a parametric curve, $(\vec x(s),p(s))$. Without loss of generality the bifurcation is to be at $s=0$.

If the Jacobian matrix $\boldsymbol{J}=\nabla\vec f$ is rank-deficient the Lyapunov-Schmidt reduction can be applied. Intutiefely the problem is reduced from a high-dimensional, possibly infinite-dimensional, one to one that has as many dimension as the deffect of $\nabla\vec f$.

The nullspace is spanned by the eigenvector, $r_0$ of $\boldsymbol{J}$, corresponding to the eigenvalue 0.

Assume that $f$ is twice differentiable w.r.t. $\vec x$, then differentiate $\vec f(\vec x(s),p(s))=0$ w.r.t. $s$ and evaluate at $s=0$

31. $\frac{\mathrm{d}}{{\mathrm{d}}s}\vec f(\vec x(s),p(s))=\nabla_x\vec fx'(s)+\partial_p\vec fp'(s)=\boldsymbol{J}\vec x'(s) + \partial_p\vec fp'(s)$

Let $\boldsymbol{H}=\nabla\nabla_x\vec f$ be the Hessian tensor.

$\frac{{\mathrm{d}}^2}{{\mathrm{d}}s}\vec f=\boldsymbol{H}\vec x'(s)\vec x'(s)+\boldsymbol{J}\vec x''(s)+\partial^2_p\vec fp'(s)+\partial_p\vec fp''(s)=0$

At $s=0$ one has $p'(0)=0$ and hence from Eq. (31) $\vec x'(0)=\vec r_0$. Projecting onto the left-eigenvector $\vec l_0$ to the eigenvalue 0. at $s=0$ one finds

$\vec l_0\boldsymbol{H}\vec r_0\vec r_0 + \vec l_0\partial_p\vec fp''(0)=0$

or with $\partial_p\vec f\not\in R(\boldsymbol{J})$

$p''(0)=-\frac{\vec l_0\boldsymbol{H}\vec r_0\vec r_0}{\vec l_0\partial_p\vec f}$.

This is a test function of wether the bifurcation is

1. subcritical ($p''(0)<0$)
2. superctirical ($p''(0)>0$)
3. transcritical

Def:

Let $P$ be the projection onto the null space

There are several cases to be distinguised.

The fold:

If the rank deficiency is one-dimensional, $\dim N(\boldsymbol{J})=1$.

The Andronov-Hopf:

If the rank deficiency is two-dimensional.

6.3 Stability exchange

At a simple fold one eigenvalue is zero. Study the eigen system of $\boldsymbol{J}=\nabla\vec f(\vec x(s),p(s))$ near $s=0$

$\boldsymbol{J}(\vec r_0+\vec w(s)) = \lambda(s)(\vec r_0+\vec w(s))$.

With a bit of analysis

$\lambda'(0)=-\vec l_0\cdot\partial_p\vec f(\vec x(0),p(0))p''(0)$

6.3.1 Extended system

Continue the extened system (linear and nonlinear)

$\vec f(\vec x(p),p)=0$

$\nabla\vec f(\vec x,p)\vec w=0$

$(\nabla\vec f(\vec x,p))^{\mathsf{T}}\vec z=0$

6.4 Continuation of boundary value problems and periodic orbits

The same procedure than above can be applied to the continuation of periodic orbits and boundary value problems. Define the implicit function to encompass the time derivative

$\vec g(\vec x,p)=\frac{\mathrm{d}}{{\mathrm{d}}t}\vec x-T\vec f(\vec x,p)=0$, with $t\in[0,1]$.

Then proceed as above. Note that the time derivative ${\mathrm{d}}/{\mathrm{d}}t$ is a linear operator which has a matrix representation just like the Jacobian. in that sense

$\nabla_x\frac{\mathrm{d}}{{\mathrm{d}}t}\vec x=\frac{\mathrm{d}}{{\mathrm{d}}t}\boldsymbol{I}$

Def. (saddle-node/fold bifurcation) #def:bif-sn

The fixpoint defined by $\vec F({\boldsymbol{x}^\mathrm{sn}},\alpha_\mathrm{sn}) = 0$ is a fold of a saddle and a node if

1. the Jacobian $J({\boldsymbol{x}^\mathrm{sn}},\alpha_\mathrm{sn})$ has one eigenvalue $\lambda_0({\boldsymbol{x}^\mathrm{sn}},\alpha_\mathrm{sn})=0$ with algebraic multiplicity one,
2. $\frac{\partial^2\lambda_0}{\partial\alpha^2}|_{{\boldsymbol{x}^\mathrm{sn}},\alpha_\mathrm{sn}}\neq0$

Def. (Hopf bifurcation) #def:bif-hopf

Assuming all but two eigenvalues have negative real parts, except for a conjugate pair of purely imaginary eigenvalues.

7 Normal forms, centre manifold and simple spiking models

7.1 Example: Dynamics in the centre manifold of a saddle-node

At arbitrary parameters the periodic soluiton adjoint to the first variational equation on the limit cycle yields the PRC. I can be calculated numerically with the continuation method. Near a bifurcation, however, if major parts of the dynamics happen in the centre manifold the PRC can be calculated analytically. As an example take the saddle-node on limit cycle bifurcation (SNLC). The spike in this case is a homoclinic orbit to a saddle-node, that enters and leaves via the semi-stable (centre) manifold that is associated with the eigenvalue $\lambda_0=0$ of the Jacobian at the saddle.

32. $\dot{\vec x}=\vec F(\vec x)+p\,\vec G(\vec x)$

${c\dot v\choose\tau\dot n}={I-I(n,v)\choose n_\infty(v)-n}$

Let there be a saddle-node (fold) bifurcation at some value $p_0$ and the saddle-node position is $\vec x_0$. In the previous lecture we noted that if the bifurcation is not transcritical (degenerated)

$\frac{{\mathrm{d}}^2 p(s)}{{\mathrm{d}}s^2}=-\frac{\vec l_0\boldsymbol{H}\vec r_0\vec r_0}{\vec l_0\partial_p\vec F}\neq0$

The change of the solution w.r.t. the pseudo-arclength parameter was

$\frac{{\mathrm{d}}\vec x}{{\mathrm{d}}s}=\vec r_0$.

Given that the linear part of the dynamics is degenerated at the bifurcation, one might look at the leading nonlinear term.

Expanding the right-hand side around saddle-node fixpoint, $\vec x_0$, yields (B. Ermentrout 1996)

$\vec F(\vec x_0+\delta\vec x)=\vec F(\vec x_0) + \boldsymbol{J}(\delta\vec x)+\boldsymbol{H}(\delta\vec x)(\delta\vec x)+...$

Establish the eigen system at the saddle-node $\boldsymbol{J}=\nabla\vec f(\vec x_0)$

$\vec l_k\boldsymbol{J}=\lambda_k\vec l_k$ and $\boldsymbol{J}\vec r_k=\lambda_k\vec r_k$ with $l_j\cdot r_k=\delta_{jk}$.

By assumption of a saddle-node the Jacobian has a simple zero with an associated eigenvector. All other eigenvalues are exponentially contracting $\lambda_k<0,\forall k>0$, while $\vec r_0$ is slow.

Def (centre subspace):

The subspace spanned by $\vec r_0$ is called the centre subspace or slow subspace.

Write the dynamics arround the saddle-node as $\vec x(t) = \vec x_0 + y\vec r_0$, and also use small changes in the bifurcation parameter $p_0+p$. Then Eq. (32) is

$\frac{{\mathrm{d}}\vec x _0}{{\mathrm{d}}t}+\dot y\vec r_0=\vec F(\vec x_0,p_0)+y\boldsymbol{J}\vec r_0+y^2\boldsymbol{H}\vec r_0\vec r_0 +p\frac{\partial\vec F}{\partial p}$.