Markovian Spiking Neural Network
Markovian Spiking Neural Network
Surrogate for linear integrate-and-fire (LIF) model
See original code.
Surrogate for FitzHugh-Nagumo (FHN) model
We refer to the system \[ \begin{cases} \frac{\mathrm{d}V}{\mathrm{d}t} &= -\frac{1}{\tau}\left((V-V_R)(V-V_{th})(V-V_M)+w \right) + I_{ext}(t) \\ \frac{\mathrm{d}w}{\mathrm{d}t} &= b(V-c)-\gamma w \end{cases} \tag{1} \]
while we'll see the parameters in the evolution of \(w\) does not matter in the Markovian surrogate. We set \(V_R=-70\)mV the resting potential, \(V_{th}=-50\)mV the firing threshold, \(V_M=40\)mV the magnitude of the peak of the spike.
The discretization of the state space is as follows: \(V\) space is discretized into \(V_M-V_I-1\) states (I found it difficult to deal with the boundary, or the stable fixed point), \([V_I+1,V_M-1]\). Here, \(V_I\) is the magnitude of hyperpolarization after a spike. It's not arbitrary and must be negative enough (so it's not ideal, and I think it's inherited from the deficit of the original FitzHugh-Nagumo model). \(w\) space is discretized into \(2\) states, \(0\) and \(w_M\). \(w_M\) is determined by \(V_I\), given by \[ w_{M} = -(V_{I}-V_{R})(V_{I}-V_{th})(V_{I}-V_{M})>0 \]
So we have \(V_{I}<V_{R}<V_{th}<V_{M}\). Now we explain how the membrane potential evolves.
- Each neuron receives independent Poisson kicks from an external source with rate \(\lambda^{E}\) or \(\lambda^{I}\).
- Each neuron receives synaptic input from within the population. Similar to the original Markovian model, denote \(H^{E}\) and \(H^{I}\) as numbers of kicks 'in play'.
- The membrane potential evolves as follows:
- If \(V_{I}< V < V_{th}\), then \(w=0\) and \(V\) evolves under the cubic term, external Poisson kicks and synaptic input from other neurons.
- If \(V_{th}\leq V < V_{M}\), then \(V\) evolves solely under the cubic term.
- If \(V=V_{M}\), then \(w\) is set to \(w_{M}\) with a rate \(\frac{1}{\tau_{w}}\), we think the neuron fires.
- If \(V_{I}<V\leqslant V_{M}\) and \(w=w_{M}\), then the neuron is in the phase of re-polarization and hyperpolarization. \(V\) evolves solely under the cubic term. We think the neuron is in a refractory period.
- If \(V=V_{I}\), then \(w\) is set to \(0\) with the rate \(\frac{1}{\tau_{w}}\).
The evolution of \(V\) is questionable and needs to be discussed.
Issue: questions about the original code
- Is the code really match the paper? Modelling the independent exponential clock of a synaptic input in play as a term in the rate matrix?
Issue: linear non-superimposability
The construction of the transition matrix involves calculation of the rate of a leaky/cubic term. Separating the leaky term is straightforward in leaky integrate-and-fire model, but it's more complicated if the leaky term is not linear.
In leaky integrate-and-fire model, the leaky term is \[ \frac{\mathrm{d}V}{\mathrm{d}t} = -\frac{1}{\tau}V \]
Notice that we deliberately set the resting potential to \(0\). In the case \(V>0\),
\[ V(t)=V(0)\exp \left( -\frac{t}{\tau} \right) \]
Set \(V(t)=V-1, V(0)=V\), then \[ t = \tau \ln \left( 1+\frac{1}{V-1} \right) \]
So the rate is given by \[ \frac{1}{t} = \frac{1}{\tau\ln \left( 1+\frac{1}{V-1} \right) } \tag{2} \]
Similarly, in the case \(V<0\), the rate of \(V\) to \(V+1\) is given by \[ \frac{1}{t} = \frac{1}{\tau \ln \left( 1-\frac{1}{V+1} \right) } \tag{3} \]
However, such rate is not linearly additive. Since the 'leaky' term in \((1)\) is a cubic term, how to calculate the rate it represents is a question.
To illustrate the linear non-superimposability, consider a quadratic 'leaky' term \[ \frac{\mathrm{d}V}{\mathrm{d}t} = - V(V+1) \]
and \(V>0\). Consider \(V(V+1)\) as a whole, the rate of \(V\) to \(V-1\) is given by \[ \frac{1}{2\ln V-\ln (V-1)-\ln (V+1)} \tag{4} \]
Separate \(V(V+1)\) into \(V^{2}\) and \(V\), the rate of \(V\) to \(V-1\) is given by \[ \frac{1}{\ln V -\ln (V-1)} + V(V-1) \tag{5} \]
Obviously, \((4)\) and \((5)\) are not equal though they are both quite close to \(V(V+1)\)
Several options to calculate the rate of the cubic term in \((1)\) are: - Consider \((V-V_R)(V-V_{th})(V-V_M)\) as a whole, separate it with \(w\). - Separate \(O(V^{3})\), \(O(V^{2})\), \(O(V)\), \(O(1)\) and \(w\) terms. - Directly use the absolute value of RHS of \((1)\) as the rate.
Issue: boundary & stable fixed points
We must be really careful of the boundary and fixed points, i.e., \(V=V_{I}, V=V_{R}, V=V_{th}, V=V_{M}\) though \(V=V_{I}\) and \(V=V_{M}\) are impossible due to our choices of state space.
In \((2)\), \(V\) cannot be \(0,1\). In \((3)\), \(V\) cannot be \(-1,0\). This is relatively simple. However, if we use the first option above to calculate the rate of the cubic term, then two functions are \[ \frac{1}{\tau \left( p \log (1+\frac{1}{V-V_{R}-1})+q\log (1+\frac{1}{V-V_{th}-1}) + r\log (1+\frac{1}{V-V_{M}-1}) \right) } \tag{6} \]
where \(p,q,r\) is the solution to \[ \begin{pmatrix} 1 & 1 & 1 \\ V_{th}+V_{M} & V_{R} + V_{M} & V_{R} + V_{th} \\ V_{th}V_{M} & V_{R}V_{M} & V_{R}V_{th} \end{pmatrix} \begin{pmatrix} p \\ q \\ r \end{pmatrix} = \begin{pmatrix} 0 \\ 0 \\ 1 \end{pmatrix} \]
Similarly, the rate of \(V\) to \(V+1\) is given by \[ \frac{1}{\tau \left( p \log (1-\frac{1}{V-V_{R}+1})+q\log (1-\frac{1}{V-V_{th}+1}) + r\log (1-\frac{1}{V-V_{M}+1}) \right) } \tag{7} \]
In \((6)\), \(V\) cannot be \(V_{R}, V_{R}+1, V_{th}, V_{th}+1, V_{M}\). In \((7)\), \(V\) cannot be \(V_{R}-1, V_{R}, V_{th}-1, V_{th}, V_{M}-1, V_{M}\). It's laborious to tackle all the cases.
Current approach to deal with the boundary: The state space of \(V\) is narrowed to \([V_{I}+1, V_{M}-1]\). The neuron evolves on the cubic term, Poisson kicks and synaptic inputs around \(V_{R}\).
- Around \(V=V_{R}\), the cubic term does not contribute if \(V=V_{R}, V_{R}+1\) when \(V\) is going to \(V-1\), and \(V=V_{R}, V_{R}-1\) when \(V\) is going to \(V+1\).
- Around \(V=V_{th}\), the cubic term does not contribute if \(V=V_{th}, V_{th}+1\) when \(V\) is going to \(V-1\), and \(V=V_{th}, V_{th}-1\) when \(V\) is going to \(V+1\). Moreover, since \(V=V_{th}\) is a fixed point when \(w=0\), the neuron still receives Poisson kicks and synaptic inputs at \(V=V_{th}\), and evolves solely under the cubic term when \(V\geqslant V_{th}+1\).
- Around \(V=V_{M}\), we think the neuron fires if \(V=V_{M}-1\), so there is no \(V\) to \(V+1\) issue in \((7)\).
- Around \(V=V_{I}\), there is actually no issue since \(V=V_{I}\) is not a singular point of \((6)\). Recall that we separate the cubic term with \(w=w_{M}\), and the rate of \(w_{M}\) is given by \(w_M/\tau\).
After discussion, we think the above two issues are not critical. We are seeking for an approximation of the original dynaimcs, the rate for the leaky/cubic term is not necessary to be accurate and various forms can be used.
Issue: scaling of the system
scaling of \(V\) and \(w\)
The original and variants of the FitzHugh-Nagumo model are not designed for modeling exact process of action potentials. They are mathematically ideal models to study the feature of periodic firing qualitatively. In specific, \(V_R, V_{th}, V_M\) are usually take small values, e.g., \(0,0.5,1.0\). The values can be seen as under unit 0.1V.
We typically look at voltage under unit mV. How can we switch between two different scalings? Denote \(V, w\) as the scale under unit 0.1V, and \(\bar{V}, \bar{w}\) as the scale under unit mV. The same treatment is applied to other parameters.
\[ \begin{cases} \frac{\mathrm{d}V}{\mathrm{d}t} &= -\frac{1}{\tau}\left((V-V_R)(V-V_{th})(V-V_M)+w \right) + I_{ext}(t) \\ \frac{\mathrm{d}w}{\mathrm{d}t} &= b(V-c)-\gamma w \end{cases} \]
\(\bar{V}=100V\), we may also assume \(\bar{w}=100w\). Then the system becomes \[ \begin{cases} \frac{\mathrm{d}\bar{V}}{\mathrm{d}t} = -\frac{1}{10000\tau}(\bar{V}-\bar{V}_{R})(\bar{V}-\bar{V}_{th})(\bar{V}-\bar{V}_{M}) -\frac{1}{\tau}\bar{w} + 100 I_{ext}(t)\\ \frac{\mathrm{d}\bar{w}}{\mathrm{d}t} = b(\bar{V}-100c) - \gamma \bar{w} \tag{8} \end{cases} \]
As can be seen, the system is not linear, so the scaling is not invariant.
scaling of time
Baltanas and Casado studied the bursting behaviour of the FitzHugh-Nagumo model subject to quasi-monochromatic noise. Specifically, the input \(I_{ext}(t)\) to the system is \(A_c + x(t)\), where \(A_c\) is a critical constant and \(x(t)\) is a noisy damped oscillation: \[ \frac{\mathrm{d}^{2}x}{\mathrm{d}t^{2}}+2\Gamma \frac{\mathrm{d}x}{\mathrm{d}t} + \omega^{2} x=\xi(t) \tag{9} \]
Here, \(\xi(t)\) is a Gaussian white noise with zero mean and correlation \(\langle \xi(t)\xi(s)\rangle = \Delta \delta(t-s)\).
The authors took \(\Gamma=0.03, \Delta=0.005, \omega=1\) and observed bursting behaviour in \((1)\). We empirically observed that with \(\tau=1\) and \(\tau_w=100\). However, the length of a spike is about 40 ms, inconsistent with a real spike, which is about 2~3 ms. How to scale \(\tau, \tau_{w}, \Gamma, \omega, \Delta\) to make the system more realistic? We simply scale \(\tau=0.1\) and \(\tau_{w}=10\).
It is a well-known fact that the solution to \((8)\) with \(\xi(t)=0\) is \[ x(t) = A_0 \exp \left( -\Gamma t \right) \cos \left( \sqrt{\omega^{2}-\Gamma^{2}} t+ \phi \right) \]
If we defined the Hamiltonian \(H(x,\dot{x})=\frac{1}{2}(\dot{x}^2+\omega^{2}x^{2})\), then the Hamiltonian equation is \[ \mathrm{d}H=-2\Gamma\left( H- \frac{\Delta}{4\Gamma} \right)\mathrm{d}t + \sqrt{\Delta H} \mathrm{d}W_t \tag{10} \]
In the long term the expected value of \(H\) tends towards \(\frac{\Delta}{4\Gamma}\). Combined with the definition of \(H\), we see that if we set \(\omega=10\) and want to keep the magnitude of oscillation, we should set \(\Delta=0.005\times 100=0.5\).
Issue: time alignment
We first discretize the state space as an ordered pair \((V,w)\). However, it's not a good idea if more variables are to be incorporated or even we have more states in \(w\) (recall that we only have \(2\) states of \(w\)). A better idea is to separate the state space of \(V\) and \(w\). Since the evolution of \(V\) and \(w\) depends on each other, the transition vector must be calculated each time the state changes.
The best we can do is it to precompute the transition probability due to the leaky/cubic term. Each time the state shall change, we grab the transition vector and add parts related to the other variable and external input.
One of the problems is that transition rate of \(V\) and \(w\) are not the same. How to align \(V\) and \(w\)? One way is to sample two random number \(\Delta t_V\), \(\Delta t_w\) from two exponential distributions with rates corresponding to \(V\) and \(w\). If \(\Delta t_V< \Delta t_w\), then we update \(V\), otherwise we update \(w\).
Issue: range of \(w\).
In the rescaled system \((8)\), the stationary value of \(\bar{w}\) is given by \(\frac{b}{\gamma}(\bar{V}-100c)\). However, we obviously don't need the range of \(\bar{w}\) to be that big from simple observation of the trajectory and the vector field.
To constrain the range of \(\bar{w}\), we set a threshold of the leaky/cubic term.