MAP estimation algorithm for phase response curves based on analysis of the observation process

Abstract

Many research groups have sought to measure phase response curves (PRCs) from real neurons. However, methods of estimating PRCs from noisy spike-response data have yet to be established. In this paper, we propose a Bayesian approach for estimating PRCs. First, we analytically obtain a likelihood function of the PRC from a detailed model of the observation process formulated as Langevin equations. Then we construct a maximum a posteriori (MAP) estimation algorithm based on the analytically obtained likelihood function. The MAP estimation algorithm derived here is equivalent to the spherical spin model. Moreover, we analytically calculate a marginal likelihood corresponding to the free energy of the spherical spin model, which enables us to estimate the hyper-parameters, i.e., the intensity of the Langevin force and the smoothness of the prior.

Appendices

Appendix 1

Here, we derive the probability density describing the stochastic process in observing phase responses. The probability density of the stochastic process described by the Langevin phase equation (3) exactly obeys the following Fokker-Planck equation (Risken 1989):

$$ \begin{array}{lll} \dot P (\phi,t) &=& - \frac{\partial}{\partial \phi} \left(D^{(1)}(\phi,t) P(\phi,t)\right)\\&& + \frac{\partial^2}{\partial \phi^2}\left(D^{(2)}(\phi,t) P(\phi,t)\right), \label{eq:AA1} \end{array} $$

(23)

where, according to Stratonovich’s definition, the drift coefficient D ⁽¹⁾(φ,t) and diffusion coefficient D ⁽²⁾(φ,t) are given as

$$ \begin{array}{lll} D^{(1)}(\phi,t) &=& 1 + Z(\phi) \cdot G(t - t_0) + Z(\phi)\cdot\frac{d Z(\phi)}{d \phi}, \\ D^{(2)}(\phi,t) &=& \sigma^2 Z(\phi)^2. \end{array}$$

Because G and σ are infinitesimal, the advance or delay of φ is vanishingly small for a short time. For that reason, we can approximate φ = t within a period T. The zero-order approximation for the drift and diffusion coefficients are

$$ \begin{array}{lll} D^{(1)}(t) &=& 1 + Z(t) \cdot G(t - t_0) + Z(t)\cdot\frac{d Z(t)}{d t}, \\ D^{(2)}(t) &=& \sigma^2 Z(t)^2. \end{array}$$

This approximation corresponds to the averaging approximation for deriving the convolution version of the phase equation (Winfree 1967; Kuramoto 1984).

The solution of Eq. (23) under the zero-order approximation is found by executing a Fourier transform with respect to φ, i.e.,

$$ \begin{array}{lll} \int^{+\infty}_{-\infty} P(\phi,t) e^{-i k \phi} d\phi = P(k,t). \end{array}$$

The equation for the Fourier transform representation is given as (replace $\partial/\partial \phi$ by ik)

$$ \begin{array}{lll} \dot P(k,t) &=& - D^{(1)}(t) i k P(k,t) - D^{(2)}(t) k^2 P(k,t) \\ &=& - \left( i k D^{(1)}(t) + k^2 D^{(2)}(t) \right) P(k,t), \label{eq:AA2} \end{array} $$

(24)

which is simpler than Eq. (23) because the partial differential equation is transformed to a set of independent ordinary differential equations. The solution of Eq. (24) is

$$ P(k,t) \!=\! C_k \exp\left(-\! i k \int_0^t D^{(1)}(t') dt' - k^2 \int_0^{t} D^{(2)}(t') dt' \right) , $$

where C _k is an initial condition. To measure the spike- timing responses, the initial distribution at t = 0 becomes P(φ,0) = δ(φ − 0). Therefore, the initial condition for the Fourier transform representation is

$$ C_k = \int^{+\infty}_{-\infty} \delta(\phi) e^{-i k \phi} d\phi = 1 . $$

By performing an inverse Fourier transform, we obtain the Gaussian (t > 0)

$$ \begin{array}{lll} P(\phi,t) &=& \frac{1}{\sqrt{4\pi\sigma^2 \int_0^{t} Z(t')^2 dt'} }\\ &\times& \exp\left(-\frac{\left(\phi\! -\! t - \int_0^{t} \left(Z(t')\cdot \frac{\partial Z(t')}{\partial t'} + Z(t') \cdot G(t' - t_0)\right) dt'\right)^2}{4 \sigma^2 \int_0^{t} Z(t')^2 dt'} \right). \end{array}$$

When t = T,

$$ \int_0^T dt' Z(t')\cdot \frac{d Z(t')}{d t'} = \frac{1}{2} \left(Z(T)^2 - Z(0)^2\right) =0 $$

is satisfied because Z(t) is a periodic function with period T. Therefore, we can ultimately derive Eq. (4).

When neurons are perturbed by periodic inputs, the periodic boundary condition should be introduced for solving the Fokker-Planck equation. However, in the formulation of the PRC observation process, as shown in Fig. 1(A) we focused on an interspike interval perturbed by a one-shot artificial input instead of periodic inputs. In this case, the free boundary conditions are suitable for the derivation of Eq. (4).

Appendix 2

Here, by inserting Fourier series representations of Z(t), G(t) and φ(t ₀) described in Eqs. (10), (11) and (12) in the posterior density (9), we derive the Fourier transformed forms of the posterior density (Tsuzurugi and Okada 2002).

By inserting Eq. (10) into the power constraint term, we get

$$ \begin{array}{lll} \int_0^{T}\!\!\! Z(t')^2 dt' \!&\!\!=\!&\! \sum\limits_{l=-M}^{M} \sum\limits_{l'=-M}^{M} \!\left(Z(l)_c \!-\! i Z(l)_s \right) \left(Z(l')_c \!+ \!i Z(l')_s \right) \\ &&\times \int_0^{T} \exp\left(i l \omega t' \right) \exp\left(-i l' \omega t' \right) dt'. \end{array}$$

The last integral is the δ function

$$ \int_0^{T} \exp\left(i l \omega t' \right) \exp\left(-i l' \omega t' \right) dt' = T \delta_{l l'}. \label{eq:AB1} $$

(25)

Therefore, we obtain

$$ \begin{array}{lll} &&{\kern-7.5pt}\int_0^{T} Z(t')^2 dt' \\ &&{\kern5pt} = T \sum\limits_{l=-M}^{M} \sum\limits_{l'=-M}^{M} \left(Z(l)_c - i Z(l)_s \right) \left(Z(l')_c + i Z(l')_s \right) \delta_{l l'} \\ &&{\kern5pt} = T \sum\limits_{l=-M}^{M} \left(Z(l)_c^2 + Z(l)_s^2 \right) . \end{array}$$

By inserting Eq. (10) into the smoothing term and replacing $\partial/\partial t$ with i l ω, we obtain

$$ \int_0^{T} \left(\frac{\partial^n Z(t')}{\partial {t'}^n}\right)^2 dt' = T \sum\limits_{l=-M}^{M} (l \omega)^{2n} \left(Z(l)_c^2 + Z(l)_s^2 \right) $$

By inserting Eqs. (10) and (11) into the convolution term, we get

$$ \begin{array}{lll} {\kern-2pt}\int_0^{T}\!\! Z(t') G(t'\!-\!t_0) dt' &\!=\!& \left(Z(l)_c \!-\! i Z(l)_s \right) \left(G(l')_c \!+\! i G(l')_s \right)\\ && \times\!\! \int_0^T \!\!\exp\!\left(i l \omega t' \right)\! \exp\left(\!-i l' \omega (t'\!\!-\!t_0) \!\right) \! dt' \!. \end{array}$$

According to Eq. (25), we can rewrite this equation as follows:

$$ \begin{array}{lll} &&{\kern-8pt}\int_0^{T}\!\! Z(t') G(t' - t_0) dt' \\ &&{\kern.2em}= T\sum\limits_{l=-M}^{M} \sum\limits_{l'=-M}^{M} \left(Z(l)_c - i Z(l)_s \right)\\ && {\kern22pt} \times\left(G(l')_c + i G(l')_s \right) \exp\left(i l' \omega t_0 \right) \delta_{l l'} \\ &&{\kern.2em}= \!T\!\!\sum\limits_{l=-M}^{M}\! \!\left(Z(l)_c \!-\! i Z(l)_s \right) \left(G(l)_c \!+\! i G(l)_s \right) \! \exp\left(i l \omega t_0 \right)\!.\\ \label{eq:AB2} \end{array} $$

(26)

Then, by inserting Eqs. (12) and (26) into the likelihood term, we get

$$ \begin{array}{lll} &&{\kern-7pt}\sum\limits_{j=0}^{2L-1} \left(\phi(\Delta t_p j) - \int_0^{T} Z(t') G(t' - \Delta t_p j) dt'\right)^2 \\[-2pt] &&{\kern12pt} =\sum\limits_{j=0}^{2L-1} \Biggl(\sum\limits_{k=-L}^{L} \left(\Phi_i(k)_c - i \Phi_i(k)_s \right) \exp\left(i k \omega \Delta t_p j \right) \\[-2pt] &&{\kern48pt} - T\sum\limits_{l=-M}^{M} \left(Z(l)_c - i Z(l)_s \right) \\[-2pt] &&{\kern48pt} \times\left(G(l)_c + i G(l)_s \right) \exp\left(i l \omega \Delta t_p j \right) \Biggr)^2. \end{array}$$

Here, we truncate terms of Φ _i (l) _c and Φ _i (l) _s higher than the L-th term as Φ _i (L + 1) _c  = Φ _i (L + 2) _c  = ⋯ = Φ _i (M) _c  = 0 and Φ _i (L + 1) _s  = Φ _i (L + 2) _s  = ⋯ = Φ _i (M) _s  = 0. In so doing, we can simplify it as follows:

$$ \begin{array}{lll} &&{\kern-7pt} \sum\limits_{j=0}^{2L-1} \left(\phi_(\Delta t_p j) - \int_0^{T} Z(t') G(t' - \Delta t_p j) dt'\right)^2 \\[-2pt] &&{\kern18pt}=\! \sum\limits_{j=0}^{2L-1} \!\!\left(\, \sum\limits_{l=-M}^{M} \biggl( \left( \Phi_i(l)_c \!-\! i \Phi_i(l)_s \right) \!-\! T \left( Z(l)_c \!-\! i Z(l)_s \right)\right. \\[-2pt] &&{\kern35pt} \left.\phantom{\sum\limits_{l=-M}^{M}} \times\left( G(l)_c + i G(l)_s \right) \biggr) \exp\left(i l \omega \Delta t_p j \right) \right)^2 \\[-2pt] &&{\kern18pt} = 2l \sum_{l=-M}^{M} \biggl| \left(\Phi_i(l)_c\! -\! i \Phi_i(l)_s \right)\! -\! T \left(Z(l)_c - i Z(l)_s \right) \\[-2pt] &&{\kern65pt}\times\left(G(l)_c\! +\! i G(l)_s \right) \biggr|^2 \end{array}$$

Appendix 3

The partition function ${\cal Z}_{ps}$ can be obtained by performing a multiple integration of the following equation.

$$ \begin{array}{lll} &&{\kern-6pt}{\cal Z}_{ps} \propto \int_{-\infty}^{\infty} dZ(0)_c \prod\limits_{l=1}^M dZ(l)_c \prod\limits_{l=1}^M Z(l)_s \\[-2pt] &&{\kern14pt}\times \!\exp\left(- \mathcal{H}\right) \!\delta\left(\!\frac{1}{\beta} \!-\! 2 \sigma^2 T \left(Z(0)_c^2 \vphantom{\sum_{l=1}^M}\right.\right.\\[-2pt] &&{\kern124pt}\left.\left.+ 2 \sum\limits_{l=1}^M \!\big(Z(l)_c^2 \!+\! Z(l)_s^2 \big) \!\right)\!\right) \end{array}$$

We next introduce an auxiliary variable x for creating a Fourier transformed representation of the delta function:

$$ \delta(t) = \frac{1}{\sqrt{2\pi}}\int_{-\infty}^{\infty} dx \exp\left(i x t \right). $$

Accordingly, the posterior density is factorized into independent Gaussians.

$$ \begin{array}{lll} {\cal Z}_{ps} &\!\propto\!& \exp\left(\!-\beta L \!\sum\limits_{i=1}^N \Phi_i(0)_c^2 \!-\! 2 \beta L\! \sum\limits_{l=1}^{M} \!\sum\limits_{i=1}^N \!\left(\Phi_i(l)_c^2 \!+\! \Phi_i(l)_s^2 \right)\!\right) \\ && \times \int_{-\infty}^{\infty} dx e^{\frac{ix}{\beta}} \int_{-\infty}^{\infty} dZ(0)_c e^{- \overline{\mathcal{H}}(0)_c} \prod\limits_{l=1}^M\\&&\times \int_{-\infty}^{\infty} dZ(l)_c e^{- \overline{\mathcal{H}}(l)_c} \prod\limits_{l=1}^M \int_{-\infty}^{\infty} dZ(l)_s e^{- \overline{\mathcal{H}}(l)_s} \\ & &\overline{\mathcal{H}}(0)_c = \frac{1}{2} \left( a(0) + i 4 \sigma^2 T x \right) Z(0)_c^2 - b_c(0) Z(0)_c, \\ & &\overline{\mathcal{H}}(l)_c = \frac{1}{2} \left(2 a(l) + i 8 \sigma^2 T x \right)Z(l)_c^2 - 2 b_c(l)Z(l)_c, \\ & &\overline{\mathcal{H}}(l)_s = \frac{1}{2} \left(2 a(l) + i 8 \sigma^2 T x \right)Z(l)_s^2 - 2 b_s(l) Z(l)_s, \\ &&{\kern19.5pt} l=1,2,\cdots,M \end{array}$$

Here, we need to integrate a Gaussian with a complex variance. The Gaussian integral formula can be extended to the following complex case:

$$ \begin{array}{lll} &&{\kern-7pt}\int_{-\infty}^{\infty} dx \exp\left(- \frac{1}{2} (A+iB) x^2 + C x\right) \\[-2pt] &&{\kern14pt} = \frac{\sqrt{2\pi}}{(A + iB)^{1/2}} \exp\left(\frac{C^2}{2\left(A + iB\right)} \right), \label{eq:AD1} \end{array} $$

(27)

in which A, B, and C are real numbers.

By transforming variables, z = (A + iB)^1/2 x, we can rewrite it into the form

$$ \begin{array}{lll} &&{\kern-7pt}\frac{(A + iB)^{1/2}}{\sqrt{2\pi} } \int_{-\infty}^{\infty} dx \exp\left(- \frac{1}{2} (A+iB) x^2 + C x\right) \\ &&{\kern12pt}= \frac{1}{\sqrt{2\pi}}\int_{\rm path} dz f(z), \\ &&{\kern-7pt}f(z)=\exp\left(- \frac{1}{2}\left(z^2 + 2 C (A + iB)^{-1/2} z \right) \right), \end{array}$$

where the path of integration becomes z = (A + iB)^1/2 x (x = − ∞ ~ + ∞ ). The angle of (A + iB)^1/2 is lower than π/4 even in the limit B→ + ∞. Let us consider the two closed loops shown in Fig. 8. Both integrations along these loops become zero because f(z) does not have singularities within these loops. We obtain the following relation:

$$ \begin{array}{lll} &&{\kern-7pt}\frac{1}{\sqrt{2\pi}}\int_{\rm path} dz f(z) \\ &&{\kern6pt}= \frac{1}{\sqrt{2\pi}}\int_{-\infty}^{\infty} dx \exp\left(- \frac{1}{2}\left( x^2 + 2 C (A + iB)^{-1/2} x \right) \right) \\ &&{\kern6pt}= \exp \left( -\frac{C^2}{2\left(A + iB\right)}\right) \end{array}$$

We can prove the complex version of the Gaussian integral formula in Eq. (27).

Fig. 8

Path of integration on the complex plane to prove Eq. (27)

According to Eq. (27), we can write the following:

$$ \begin{array}{lll} &&{\kern-35pt}{\cal Z}_{\!ps} \!\propto \exp\left(\!-\beta L \!\sum\limits_{i=1}^N \!\Phi_i(0)_c^2 \!-\! 2 \beta L \!\sum\limits_{l=1}^{M} \!\sum\limits_{i=1}^N \!\left(\Phi_i(l)_c^2 \!+\! \Phi_i(l)_s^2 \right)\!\right) \\ &&{\kern-14pt} \times \!\!\int_{-\infty}^{\infty} \!dx e^{\frac{ix}{\beta}} \!\frac{\sqrt{2\pi}}{ (a(0) + i 4 \sigma^2 T x)^{1/2}}\exp\!\left(\!\frac{b_c(0)^2}{a(0) \!+ \!i 4 \sigma^2 T x} \!\right) \\ &&{\kern-14pt} \times \prod\limits_{l=1}^M \frac{2\pi}{2 a(l) + i 8 \sigma^2 T x}\exp\left(\frac{b(l)^2}{a(l) + i 4 \sigma^2 T x} \right) \\ b(l)^2 &=& b_c(l)^2 + b_s(l)^2 \\ &=& 4 \beta^2 L^2 T^2 \left(G(l)_c^2 + G(l)_s^2 \right)\\&&\times \left(\left(\sum\limits_{i=1}^{N}\Phi_i(l)_c \right)^2 + \left(\sum\limits_{i=1}^{N}\Phi_i(l)_s \right)^2\right) \end{array} $$

(28)

According to Euler’s formula, we finally obtain Eq. (17).