Spatio-temporal spike train analysis for large scale networks using the maximum entropy principle and Monte Carlo method

We consider a network of N neurons. We assume there is a minimal timescale δ, such that a neuron can fire at most one spike within a time window of size δ. To each neuron k and discrete time n, we associate a spike variable: ω_k(n) = 1 if neuron k fires at time n, and ω_k(n) = 0 otherwise. The state of the entire network in time bin n is thus described by a vector $\omega (n)\stackrel{\mathrm{def}}{=}\left [\hspace{0.167em} {\omega }_{k}(n)\hspace{0.167em} \right ]_{k=1}^{N}$, called a spiking pattern.

A spike block, which describes the activity of the whole network between moments of time n₁ and n₂, is a finite ordered list of such vectors, written as:

$$\begin{eqnarray*} &&{\omega }_{{n}_{1}}^{{n}_{2}}=\left \{ \hspace{0.167em} \omega (n)\hspace{0.167em} \right \} _{\{ {n}_{1}\leq n\leq {n}_{2}\} }. \end{eqnarray*}$$

The range of a block is n₂ − n₁ + 1, the number of time steps from n₁ to n₂. Here is an example of a spike block of range 5 with N = 4 neurons.

$$\begin{eqnarray*} \left (\begin{array}{@{}ccccc@{}} \displaystyle 1&\displaystyle 1&\displaystyle 0&\displaystyle 1&\displaystyle 0\\ \displaystyle 0&\displaystyle 1&\displaystyle 0&\displaystyle 1&\displaystyle 0\\ \displaystyle 0&\displaystyle 0&\displaystyle 0&\displaystyle 1&\displaystyle 0\\ \displaystyle 0&\displaystyle 1&\displaystyle 0&\displaystyle 1&\displaystyle 1 \end{array}\right ). \end{eqnarray*}$$

A spike train or raster is a spike block ${\omega }_{0}^{T-1}$ from some initial time 0 to some final time T − 1. To simplify notation we simply write ω for a spike train. We note Ω = { 0,1 }^NT for the set of all possible spike trains.

We call an observable a function:

$$\begin{eqnarray} &&\mathcal{O}(\omega )=\prod _{u=1}^{r}{\omega }_{{k}_{u}}({n}_{u}), \end{eqnarray} \tag{ 1 }$$

i.e. a product of binary spike events where k_u is a neuron index and n_u a time index, with u = 1,...,r, for some integer r > 0. Typical choices of observables are ω_k1(n₁), which is 1 if neuron k₁ fires at time n₁ and which is 0 otherwise; ω_k1(n₁) ω_k2(n₂), which is 1 if neuron k₁ fires at time n₁ and neuron k₂ fires at time n₂ and which is 0 otherwise. Another example is ω_k1(n₁) (1 − ω_k2(n₂)), which is 1 if neuron k₁ fires at time n₁ and neuron k₂ is silent at time n₂. This example emphasizes that observables are able to consider events where some neurons are silent.

We say that an observable has range R if it depends on R consecutive spike patterns, e.g. $\mathcal{O}(\omega )=\mathcal{O}({\omega }_{0}^{R-1})$. We consider here that observables do not depend explicitly on time (time-translation invariance of observables). As a consequence, for any time n, $\mathcal{O}({\omega }_{0}^{R-1})=\mathcal{O}({\omega }_{n}^{n+R-1})$ whenever ${\omega }_{0}^{R-1}={\omega }_{n}^{n+R-1}$.

It is common in the study of spike trains to attempt to detect some statistical regularity. Spike train statistics is assumed to be summarized by a hidden probability μ characterizing the probability of spatio-temporal spike patterns: μ is defined as soon as the probability $\mu \left [\hspace{0.167em} {\omega }_{{n}_{1}}^{{n}_{2}}\hspace{0.167em} \right ]$ of any block ${\omega }_{{n}_{1}}^{{n}_{2}}$ is known. We assume that μ is time-translation invariant: for any time n, $\mu \left [\hspace{0.167em} {\omega }_{0}^{R-1}\hspace{0.167em} \right ]=\mu \left [\hspace{0.167em} {\omega }_{n}^{n+R-1}\hspace{0.167em} \right ]$, whenever ${\omega }_{0}^{R-1}={\omega }_{n}^{n+R-1}$.

Equivalently, μ allows the computation of the average of the observables. We denote μ[  ] as the average of the observable under μ. If (ω) = ω_k1(n₁), then μ[  ] is the firing rate of neuron k₁ (it does not depend on n₁ from the time-translation invariance hypothesis); if  = ω_k1(n₁) ω_k2(n₂), then μ[  ] is the probability that neurons k₁ and k₂ fire during the time span n₂ − n₁. Additionally, μ[ ω_k1(0) ω_k2(0) ] − μ[ ω_k1(0) ]μ[ ω_k2(0) ] represents the instantaneous pairwise correlation between the neurons k₁ and k₂.

There are several methods which allow the computation or estimation of μ. In the following we shall assume that neural activity is described by a Markov process with memory depth D and positive time-translation invariant transition probabilities $P\left [\hspace{0.167em} \omega (D)\hspace{0.167em} \left \vert \hspace{0.167em} {\omega }_{0}^{D-1}\right .\hspace{0.167em} \right ]\gt 0$. From the assumption $P\left [\hspace{0.167em} \omega (D)\hspace{0.167em} \left \vert \hspace{0.167em} {\omega }_{0}^{D-1}\right .\hspace{0.167em} \right ]\gt 0$, this chain has a unique invariant probability μ such that, for any n > D, and any block ${\omega }_{0}^{n}$:

$$\begin{eqnarray} &&\mu \left [\hspace{0.167em} {\omega }_{0}^{n}\hspace{0.167em} \right ]=\prod _{l=0}^{n-D}P\left [\hspace{0.167em} \omega (D+l)\hspace{0.167em} \left \vert \hspace{0.167em} {\omega }_{l}^{D+l-1}\right .\hspace{0.167em} \right ]\mu \left [\hspace{0.167em} {\omega }_{0}^{D-1}\hspace{0.167em} \right ]. \end{eqnarray} \tag{ 2 }$$

Therefore, knowing the transition probabilities (corresponding to blocks ${\omega }_{0}^{D}$ of range D + 1) and μ (which can be determined as well from the transition probabilities as exposed in section 2.3.1), the probability of larger blocks can be computed. Equation (2) makes explicit the role of memory in statistics of spike blocks, via the product of transition probabilities and the probability of the initial block $\mu \left [\hspace{0.167em} {\omega }_{0}^{D-1}\hspace{0.167em} \right ]$.

In contrast, if D = 0, the probability to have the spike pattern ω(D) does not depend on the past activity of the network (memory-less case). In this case $P\left [\hspace{0.167em} \omega (D+l)\hspace{0.167em} \left \vert \hspace{0.167em} {\omega }_{l}^{D+l-1}\right .\hspace{0.167em} \right ]$ becomes μ[ ω(l) ] and the probability (2) of a block becomes:

$$\begin{eqnarray} &&\mu \left [\hspace{0.167em} {\omega }_{0}^{n}\hspace{0.167em} \right ]=\prod _{l=0}^{n}\mu \left [\hspace{0.167em} \omega (l)\hspace{0.167em} \right ]. \end{eqnarray} \tag{ 3 }$$

Therefore, in the memory-less case, spikes occurring at different times are independent.

This emphasizes the deep difference between the case D = 0 and the case D > 0.

Let us assume that we are given an experimental raster of length T, such that ${\omega }_{0}^{T-1}$. The estimation of spike statistics has to be done on this sample. In the context of the maximum entropy principle, where statistics is assumed to be time-translation invariant, statistics of events is obtained via the time-average. The time-average or empirical average of an observable in a raster ω of length T is denoted by ${\pi }_{\omega }^{(T)}\left [\hspace{0.167em} \mathcal{O}\hspace{0.167em} \right ]$. For example, if  = ω_k(0) the time-average ${\pi }_{\omega }^{(T)}\left [\hspace{0.167em} \mathcal{O}\hspace{0.167em} \right ]=1/(T-1)\hspace{0.167em} {\mathop{\sum }\nolimits }_{n=0}^{T-1}{\omega }_{k}(n)$ is the firing rate of neuron k, estimated on the experimental raster ω.

The empirical average is a random variable, depending on the raster ω as well as on the time length of the sample, and it has Gaussian fluctuations whose amplitude tends to 0 as T → +∞ like $1/\sqrt{T}$. This is the case, e.g. for the empirical averages obtained from several spike train acquired with several repetitions.

If one has N neurons and wants to consider spike block events within R time steps, one has 2^NR possible states. For a reasonable multielectrodes array (MEA) sample, N = 100, R = 3 (for a time lag of 30 ms with a 10 ms binning), this gives 2³⁰⁰ ∼ 4 × 10¹⁸⁰, which is considerably more than the expected number of particles in the (visible) universe. Taking into account the huge number of states in the set of blocks, it is clear that any method requiring the extensive description of the state space will fail as NR grows. Additionally, while the accessible state space is huge, the observed state space (e.g. in an experimental raster) is rather small. For example, in a MEA raster for a retina experiment, the sample size T is about 10⁶–10⁷, which is quite a bit less than 2^NR. As a matter of fact, any reasonable estimation method must take this small-sample constraint into account. As we show in section 2.2, the maximum entropy principle and the related notion of Gibbs distributions allows us to take these aspects into account.

Following (2.1) the goal is to find a probability distribution μ such that:

• μ is inferred from an empirical raster ω, by computing the empirical average of a set of ad hoc observables _k, k = 1,...,K. One asks that the average of _k with respect to μ satisfies:
  $$\begin{eqnarray} &&\mu \left [\hspace{0.167em} {\mathcal{O}}_{k}\hspace{0.167em} \right ]={\pi }_{\omega }^{(T)}\left [\hspace{0.167em} {\mathcal{O}}_{k}\hspace{0.167em} \right ],\tqs k=1,\ldots ,K. \end{eqnarray} \tag{ 4 }$$
  The mean of _k predicted by μ is equal to the mean computed on the experimental raster. μ is called a 'model' in the following. The set of observables _k defines the model.
• μ has to be 'as simple as possible', with the least structure and a minimum number of tunable parameters. In the maximum entropy paradigm [27] these issue are (partly) solved by seeking a probability distribution μ which maximizes the entropy under the constraints (4). The entropy is defined explicitly below (see equations (5) and (21)).
• From the knowledge of μ one can compute the probability of arbitrary blocks (e.g. via equation (2)) and the average of observables other than the _k.

Remark Assume that we want to select observables _k in the set of all possible observables with range R. For N neurons there are 2^NR possible choices. When NR increases, the number of possible observables will quickly exceed the number of samples available in the recording. Including all of them in the model would overfit the data. Therefore, one has to guide the choice of observables by additional criteria. We now review some of the criteria which have been used by other authors.

In a seminal paper, Schneidman et al [55] aimed at unravelling the role of instantaneous pairwise correlations in retina spike trains. Although these correlations are weak, researchers investigated whether they play a more significant role in spike train statistics than firing rates.

Firing rates correspond to the average of observables of the form ω_i(0),i = 1,...,N (the time index 0 comes from the assumed time-translation invariance) while instantaneous pairwise correlations correspond to averages of observables of the form ω_i(0)ω_j(0),1 ≤ i < j ≤ N. Analysing the role of pairwise correlations in spike train statistics, compared to the firing rate, amounts therefore to comparing two models, defined by two different types of observables.

Note that all of these observables correspond to spatial events occurring at the same time. They give no information on how the spike patterns at a given time depend on the past activity. This situation corresponds to a memory-less model (D = 0 in section 2.1.3), where transition probabilities do not depend on the past. As a consequence the sought probability μ weights blocks of range 1, and the probability of blocks with larger range is given by (3): spike patterns at successive time steps are independent in spatial models.

In this case, the entropy of μ is given by:

$$\begin{eqnarray} &&S(\mu )=-{\sum }_{\omega (0)}\mu \left [\hspace{0.167em} \omega (0)\hspace{0.167em} \right ]\log \mu \left [\hspace{0.167em} \omega (0)\hspace{0.167em} \right ]. \end{eqnarray} \tag{ 5 }$$

The natural log could be replaced by the logarithm in base 2.

Now, the maximum entropy principle of Jaynes [27] corresponds to seeking a probability μ which maximizes S(μ) under the constraints (4). It can be shown (see section 2.3 for the general statement) that this maximization problem is equivalent to the following Lagrange problem: maximizing the quantity S(μ) + μ[ _β ], where _β, called a potential, is given by:

$$\begin{eqnarray} &&{\mathcal{H}}_{\boldsymbol{\beta }}=\sum _{k=1}^{K}{\beta }_{k}{\mathcal{O}}_{k}. \end{eqnarray} \tag{ 6 }$$

The β_k are real numbers and free parameters. μ[ _β ] is the average of _β with respect to μ. Since _β is a linear combination of observables we have $\mu \left [\hspace{0.167em} {\mathcal{H}}_{\boldsymbol{\beta }}\hspace{0.167em} \right ]={\mathop{\sum }\nolimits }_{k=1}^{K}{\beta }_{k}\mu \left [\hspace{0.167em} {\mathcal{O}}_{k}\hspace{0.167em} \right ]$. If the _k have finite range and are > − ∞ and if the β_k are finite then it can be shown (see section 2.3) that there is only one probability μ, depending on the β_k, which solves the maximization problem. It is called a Gibbs distribution.

In this context (D = 0) it reads:

$$\begin{eqnarray} &&\mu \left [\hspace{0.167em} \omega (0)\hspace{0.167em} \right ]=\frac{{\mathrm{e}}^{{\mathcal{H}}_{\boldsymbol{\beta }}(\omega (0))}}{{Z}_{\boldsymbol{\beta }}}, \end{eqnarray} \tag{ 7 }$$

where the normalization factor

$$\begin{eqnarray} &&{Z}_{\boldsymbol{\beta }}={\sum }_{\omega (0)}\hspace{0.167em} {\mathrm{e}}^{{\mathcal{H}}_{\boldsymbol{\beta }}(\omega (0))} \end{eqnarray} \tag{ 8 }$$

is the so-called partition function.

To match (4) the parameters β_k have to be tuned. This can be done thanks to the following property of Z_β:

$$\begin{eqnarray} &&\mu \left [\hspace{0.167em} {\mathcal{O}}_{k}\hspace{0.167em} \right ]=\frac{\partial \log {Z}_{\boldsymbol{\beta }}}{\partial {\beta }_{k}}. \end{eqnarray} \tag{ 9 }$$

Thus the β_k have to be tuned so that μ matches (4) as well as (9). It turns out that logZ_β is convex with respect to the β_k, so the problem has a unique solution.

Note that logZ_β does not only allow us to obtain the averages of the observables, it also allows us to characterize fluctuations. If a raster is distributed according to the Gibbs distribution (7), then, as pointed out in section 2.1.4, the empirical average of an observable has fluctuations. One can show that these fluctuations are Gaussian (central limit theorem). The joint probability of ${\pi }_{\omega }^{(T)}\left [\hspace{0.167em} {\mathcal{O}}_{k}\hspace{0.167em} \right ]$, k = 1,...,K is Gaussian, with mean μ[ _k ] given by (9) and covariance matrix Σ/T, where the matrix Σ has entries:

$$\begin{eqnarray} &&{\Sigma }_{k l}=\frac{{\partial }^{2}\log {Z}_{\boldsymbol{\beta }}}{\partial {\beta }_{k}\partial {\beta }_{l}}. \end{eqnarray} \tag{ 10 }$$

Let us now discuss what this principle gives in the two cases considered by Schneidman et al

• (i)  
  Only firing rates are constrained. Then:

  $$\begin{eqnarray*} &&{\mathcal{H}}_{\boldsymbol{\beta }}(\omega (0))=\sum _{k=1}^{N}{\beta }_{k}{\omega }_{k}(0). \end{eqnarray*}$$

  It can be shown that the corresponding probability μ is:

  $$\begin{eqnarray*} &&\mu \left [\hspace{0.167em} \omega (0)\hspace{0.167em} \right ]=\prod _{k=1}^{N}\frac{{\mathrm{e}}^{{\beta }_{k}\hspace{0.167em} {\omega }_{k}(0)}}{1+{\mathrm{e}}^{{\beta }_{k}}}. \end{eqnarray*}$$

  Thus, the corresponding statistical model is such that spikes emitted by distinct neurons at the same time are independent. The parameter β_k is directly related to the firing rate r_k since r_k = μ[ ω_k(0) = 1 ] = e^β_k/(1 + e^β_k), so that we have:

  $$\begin{eqnarray*} &&\mu \left [\hspace{0.167em} {\omega }_{0}^{n}\hspace{0.167em} \right ]=\prod _{l=0}^{n}\hspace{0.167em} \prod _{k=1}^{N}\hspace{0.167em} {r}_{k}^{{\omega }_{k}(l)}\hspace{0.167em} (1-{r}_{k})^{1-{\omega }_{k}(l)}, \end{eqnarray*}$$

  the classical probability of coin tossing with independent probabilities (Bernoulli model). Thus, fixing only the rates as constraints, the maximum entropy principle leads us to analyse spike statistics as if each spike were thrown randomly and independently, as with coin tossing. This is the 'most random model', which has the advantage of making as few hypotheses as possible. However, when only constrained with mean firing rates, the prediction of even small spike blocks in the retina was not successful. This was expected since this model assumes independence between neurons, an assumption that has been proven wrong in earlier studies (e.g. [50]).
• (ii)  
  Firing rates and pairwise correlations are constrained. In the second model, Schneidman et al constrained the maximum entropy model with both mean firing rates and instantaneous pairwise correlations between neurons. In this case,

  $$\begin{eqnarray*} &&{\mathcal{H}}_{\boldsymbol{\beta }}(\omega (0))=\sum _{k=1}^{N}{\beta }_{k}\hspace{0.167em} {\omega }_{k}(0)+\sum _{k,l=1}^{N}{\beta }_{k l}\hspace{0.167em} {\omega }_{k}(0)\hspace{0.167em} {\omega }_{l}(0). \end{eqnarray*}$$

  Here the potential can be identified with the Hamiltonian of a magnetic system with binary spins. It is thus often called the 'Ising model' in the spike train analysis literature, although the original Ising model has constant and positive couplings [21]. The corresponding statistical model is the least structured model respecting these first-order and second-order pairwise instantaneous constraints. The number of parameters is of the order³ of N², to be compared with the 2^N possible patterns.

Schneidman et al showed that the Ising model successfully predicts spatial patterns, a result which was confirmed by [61] (see [39] for a review). Other works have used the same method and found also a good prediction in cortical structure in vitro [69], and in the visual cortex in vivo [79]. Later on, several authors considered higher-order terms still corresponding to D = 0 [40, 55, 73, 18]. Note that these results have been obtained on relatively small subsets of neurons (usually groups of 10). An interesting challenge is to test how these results hold for larger subsets of neurons, and if other constraints have to be added [19](Tkacik et al, in preparation).

These models are only designed to predict the occurrence of 'spatial' patterns, lying within one time bin. The use of spatial observables naturally leads to a time independence assumption where the probability of occurrence of a spatio-temporal pattern is given by the product (3). Tang et al [69] tried to predict the temporal statistics of the neural activity with such a model and showed that it does not give a faithful description of temporal statistics. The idea to consider spatio-temporal observables then naturally emerges with the problem of generalizing the probability equation (7) to that case.

From the statistical mechanics point of view, a natural extension consists of considering the space of rasters Ω as a lattice where one dimension is 'space' (neurons index) and the other is time. The idea is then to consider a potential, still of the form (6), but where the observables correspond to spatio-temporal events. We assume that _β has range R = D + 1, 0 ≤ D <+ ∞. The potential of a spike block ${\omega }_{0}^{n}$, n ≥ D is:

$$\begin{eqnarray} &&{\mathcal{H}}_{\boldsymbol{\beta }}\left (\hspace{0.167em} {\omega }_{0}^{n}\hspace{0.167em} \right )=\sum _{l=0}^{n-D}{\mathcal{H}}_{\boldsymbol{\beta }}\left (\hspace{0.167em} {\omega }_{l}^{D+l}\hspace{0.167em} \right ). \end{eqnarray} \tag{ 11 }$$

On this basis, restricting to the case where D = 1 (one time step memory depth) Marre et al have proposed in [35] to construct a Markov chain where transition probabilities P[ ω(l + 1) | ω(l) ] are proportional to ${\mathrm{e}}^{{\mathcal{H}}_{\boldsymbol{\beta }}({\omega }_{l}^{l+1})}$. If μ is the invariant probability of that chain, the application of (2) leads to probability of blocks $\mu \left [\hspace{0.167em} {\omega }_{0}^{n}\hspace{0.167em} \right ]$, proportional to ${\mathrm{e}}^{{\mathcal{H}}_{\boldsymbol{\beta }}\left (\hspace{0.167em} {\omega }_{0}^{n}\hspace{0.167em} \right )}$: the probability of a block is proportional to the exponential of its potential ('energy'). This approach is therefore quite natural from the statistical mechanics point of view.

The main problem, however, is 'what is the proportionality coefficient?' As shown in [35], the normalization of conditional probabilities does not reduce to the mere division by a constant partition function. This normalization factor is itself dependent on the past activity.

To overcome this dependency, Marre et al assumed that the activity respected a detailed balance. In this particular case, it can be shown that the normalization factor becomes, again, a constant. But this is an important reduction that could have implications for the interpretation of the data: for example, with this simplification, it is not possible to give an account of asymmetric cross correlograms.

In this section we show how one can generate a Markov chain where transition probabilities are proportional to ${\mathrm{e}}^{{\mathcal{H}}_{\boldsymbol{\beta }}({\omega }_{l}^{l+D})}$, for a potential _β corresponding to spatio-temporal events. We also solve the normalization problem. This construction is well known and is based on the so-called transfer matrix (see e.g. [21] for a presentation in the context of statistical physics; [44] for a presentation in the context of ergodic theory and [77] for a presentation in the context of spike train analysis).

This matrix is constructed as follows. Consider two spike blocks w₁,w₂ of range D ≥ 1. The transition w₁ → w₂ is legal if w₁ has the form $\omega (0){\omega }_{1}^{D-1}$ and w₂ has the form ${\omega }_{1}^{D-1}\omega (D)$. The vectors ω(0),ω(D) are arbitrary but the block ${\omega }_{1}^{D-1}$ is common. Here is an example of a legal transition:

$$\begin{eqnarray*} {w}_{1}=\left [\begin{array}{@{}ccc@{}} \displaystyle 0&\displaystyle 0&\displaystyle 1\\ \displaystyle 0&\displaystyle 1&\displaystyle 1 \end{array}\right ];\tqs {w}_{2}=\left [\begin{array}{@{}ccc@{}} \displaystyle 0&\displaystyle 1&\displaystyle 1\\ \displaystyle 1&\displaystyle 1&\displaystyle 0 \end{array}\right ]. \end{eqnarray*}$$

Here is an example of a forbidden transition

$$\begin{eqnarray*} {w}_{1}=\left [\begin{array}{@{}ccc@{}} \displaystyle 0&\displaystyle 0&\displaystyle 1\\ \displaystyle 0&\displaystyle 1&\displaystyle 1 \end{array}\right ];\tqs {w}_{2}=\left [\begin{array}{@{}ccc@{}} \displaystyle 0&\displaystyle 1&\displaystyle 1\\ \displaystyle 0&\displaystyle 1&\displaystyle 0 \end{array}\right ]. \end{eqnarray*}$$

Any block ${\omega }_{0}^{D}$ of range R = D + 1 can be viewed as a legal transition from the block ${w}_{1}={\omega }_{0}^{D-1}$ to the block ${w}_{2}={\omega }_{1}^{D}$ and in this case we write ${\omega }_{0}^{D}\sim {w}_{1}{w}_{2}$.

The transfer matrix is defined as:

$$\begin{eqnarray} {\mathcal{L}}_{{w}_{1},{w}_{2}}=\left \{ \begin{array}{@{}cc@{}} \displaystyle {\mathrm{e}}^{{\mathcal{H}}_{\boldsymbol{\beta }}({\omega }_{0}^{D})}&\displaystyle \text{if}\hspace{0.167em} {w}_{1},{w}_{2}\hspace{0.167em} \text{is legal with }{\omega }_{0}^{D}\sim {w}_{1}{w}_{2}\\ \displaystyle 0,&\displaystyle \text{otherwise}. \end{array}\right . \end{eqnarray} \tag{ 12 }$$

From the matrix the transition matrix of a Markov chain can be constructed, as we now show. Since observables are assumed to be bounded from below, ${\mathcal{H}}_{\boldsymbol{\beta }}({\omega }_{0}^{D})\gt -\infty$, thus ${\mathrm{e}}^{{\mathcal{H}}_{\boldsymbol{\beta }}({\omega }_{0}^{D})}\gt 0$ for each legal transition. As a consequence of the Perron–Frobenius theorem [20, 56], has a unique real positive eigenvalue s_β, strictly larger than the modulus of the other eigenvalues (with a positive gap), and with right, R, and left, L, eigenvectors: R = s_βR, L = s_βL, or, equivalently⁴:

$$\begin{eqnarray*} {\sum }_{\omega (D)\in \left \{ \hspace{0.167em} 0,1 \hspace{0.167em} \right \} ^{N}}{\mathrm{e}}^{{\mathcal{H}}_{\boldsymbol{\beta }}\left (\hspace{0.167em} \omega _{0}^{D}\hspace{0.167em} \right )}R\left (\hspace{0.167em} \omega _{1}^{D}\hspace{0.167em} \right )&=&{s}_{\boldsymbol{\beta }}R\left (\hspace{0.167em} \omega _{0}^{D-1}\hspace{0.167em} \right );\\ {\sum }_{\omega (0)\in \left \{ \hspace{0.167em} 0,1 \hspace{0.167em} \right \} ^{N}}L\left (\hspace{0.167em} \omega _{0}^{D-1}\hspace{0.167em} \right ){\mathrm{e}}^{{\mathcal{H}}_{\boldsymbol{\beta }}\left (\hspace{0.167em} \omega _{0}^{D}\hspace{0.167em} \right )}&=&{s}_{\boldsymbol{\beta }}L\left (\hspace{0.167em} \omega _{1}^{D}\hspace{0.167em} \right ). \end{eqnarray*}$$

These eigenvectors have strictly positive entries R( ⋅ ) > 0, L( ⋅ ) > 0, functions of blocks of range D. They can be chosen so that the scalar product 〈 L,R 〉 = 1. We define:

$$\begin{eqnarray} &&\mathcal{P}({\mathcal{H}}_{\boldsymbol{\beta }})=\log {s}_{\boldsymbol{\beta }}. \end{eqnarray} \tag{ 13 }$$

called 'topological pressure'. We discuss the origin of this term and its properties in section 2.3.2.

To define a Markov chain from the transfer matrix (equation (12) ) we introduce the normalized potential:

$$\begin{eqnarray} &&\phi ({\omega }_{0}^{D})={\mathcal{H}}_{\boldsymbol{\beta }}({\omega }_{0}^{D})-{\mathcal{G}}_{\boldsymbol{\beta }}({\omega }_{0}^{D}) \end{eqnarray} \tag{ 14 }$$

with:

$$\begin{eqnarray} &&{\mathcal{G}}_{\boldsymbol{\beta }}({\omega }_{0}^{D})=\log R\left (\hspace{0.167em} {\omega }_{0}^{D-1}\hspace{0.167em} \right )-\log R\left (\hspace{0.167em} {\omega }_{1}^{D}\hspace{0.167em} \right )+\log {s}_{\boldsymbol{\beta }}, \end{eqnarray} \tag{ 15 }$$

and a family of transition probabilities:

$$\begin{eqnarray} &&P\left [\hspace{0.167em} \omega (D)\hspace{0.167em} \left \vert \hspace{0.167em} {\omega }_{0}^{D-1}\right .\hspace{0.167em} \right ]\stackrel{\mathrm{def}}{=}{\mathrm{e}}^{\phi ({\omega }_{0}^{D})}\gt 0. \end{eqnarray} \tag{ 16 }$$

These transition probabilities define a Markov chain which admits a unique invariant probability:

$$\begin{eqnarray} &&\mu ({\omega }_{0}^{D-1})=R\left (\hspace{0.167em} {\omega }_{0}^{D-1}\hspace{0.167em} \right )L\left (\hspace{0.167em} {\omega }_{0}^{D-1}\hspace{0.167em} \right ). \end{eqnarray} \tag{ 17 }$$

From the general form of block probabilities (2) the probability of blocks of depth n ≥ D is, in this case:

$$\begin{eqnarray} &&\mu \left [\hspace{0.167em} {\omega }_{0}^{n}\hspace{0.167em} \right ]={\mathrm{e}}^{\sum _{l=0}^{n-D}\phi \left (\hspace{0.167em} {\omega }_{l}^{D+l}\hspace{0.167em} \right )}\mu \left [\hspace{0.167em} {\omega }_{0}^{D-1}\hspace{0.167em} \right ]. \end{eqnarray} \tag{ 18 }$$

thus, from (17),(14),(15):

$$\begin{eqnarray} &&\mu \left [\hspace{0.167em} {\omega }_{0}^{n}\hspace{0.167em} \right ]=\frac{{\mathrm{e}}^{{\mathcal{H}}_{\boldsymbol{\beta }}\left (\hspace{0.167em} {\omega }_{0}^{n}\hspace{0.167em} \right )}}{{s}_{\boldsymbol{\beta }}^{n-D+1}}R\left (\hspace{0.167em} {\omega }_{n-D+1}^{n}\hspace{0.167em} \right )L\left (\hspace{0.167em} {\omega }_{0}^{D-1}\hspace{0.167em} \right ), \end{eqnarray} \tag{ 19 }$$

where ${\mathcal{H}}_{\boldsymbol{\beta }}\left (\hspace{0.167em} {\omega }_{0}^{n}\hspace{0.167em} \right )$ is given by (11).

• (1)  
  We have been able to compute the probability of any blocks ${\omega }_{0}^{n}$. It is proportional to ${\mathrm{e}}^{{\mathcal{H}}_{\boldsymbol{\beta }}\left (\hspace{0.167em} {\omega }_{0}^{n}\hspace{0.167em} \right )}$ and the proportionality factor has been computed. In the general case of spatio-temporal events, it depends on ${\omega }_{0}^{D-1}$ and ${\omega }_{n-D+1}^{n}$.The same arises in statistical mechanics when dealing with boundary conditions. The forms (18), (19), are reminiscent of Gibbs distributions on spin lattices, with lattice translation invariant probability distributions given specific boundary conditions. Given a spin potential of spatial range n, the probability of a spin block depends upon the state of the spin block, as well as the spin states in a neighbourhood of that block. The conditional probability of this block given a fixed neighbourhood is the exponential of the energy characterizing physical interactions, within the block, as well as interactions with the boundaries. In (18), spins are replaced by spiking patterns; space is replaced by time. Spatial boundary conditions are here replaced by the dependence upon ${\omega }_{0}^{D-1}$ and ${\omega }_{n-D+1}^{n}$.As a consequence, as soon as one is dealing with spatio-temporal events the normalization of conditional probabilities does not reduce to the mere division by:

  $$\begin{eqnarray} &&{Z}_{n}={\sum }_{\omega _{0}^{n}}{\mathrm{e}}^{{\mathcal{H}}_{\boldsymbol{\beta }}\left (\hspace{0.167em} \omega _{0}^{n}\hspace{0.167em} \right )}, \end{eqnarray} \tag{ 20 }$$

  as easily checked in (19).
• (2)  
  The topological pressure obeys nevertheless:

  $$\begin{eqnarray*} &&\mathcal{P}({\mathcal{H}}_{\boldsymbol{\beta }})={\lim }_{n\rightarrow +\infty }\frac{1}{n}\log {Z}_{n}, \end{eqnarray*}$$

  and is analogous to a thermodynamic potential density (free energy, free enthalpy, pressure). This analogy is also clear in the variational principle (23) below. To our best knowledge the term 'topological pressure' has its roots in the thermodynamic formalism of hyperbolic (chaotic) maps [54, 44, 4]. In this context, this function can be computed as the grand potential of the grand canonical ensemble, as a cycle expansion over unstable periodic orbits. It is therefore equivalent to a pressure⁵ depending on topological properties (periodic orbits).
• (3)  
  In the case D = 0 the Gibbs distribution reduces to (7). One can indeed easily show that:

  $$\begin{eqnarray*} &&\exp {\mathcal{G}}_{\boldsymbol{\beta }}={s}_{\boldsymbol{\beta }}={\sum }_{\omega (0)}\hspace{0.167em} {\mathrm{e}}^{{\mathcal{H}}_{\boldsymbol{\beta }}(\omega (0))}={Z}_{\boldsymbol{\beta }}, \end{eqnarray*}$$

  where Z_β is the partition function (8). Additionally, since spike patterns occurring at distinct time are independent in the D = 0 case, Z_n in (20) can be written as ${Z}_{n}={Z}_{\boldsymbol{\beta }}^{n}$ so that (_β) = logZ_β.
• (4)  
  In the general case of spatio-temporal constraints, the normalization requires the consideration of a normalizing function _β depending as well on the blocks ${\omega }_{0}^{D}$. Thus, in addition to function _β normalization introduces a second function of spike blocks. This consequently increases the complexity of Gibbs potentials and Gibbs distributions compared to the spatial (D = 0) case where _β reduces to a constant.

We now show that the probability distribution defined this way solves the variational problem 'maximizing entropy under constraints'.

We define the entropy rate (or Kolmogorov–Sinai entropy):

$$\begin{eqnarray} &&h\left [\hspace{0.167em} \mu \hspace{0.167em} \right ]\hspace{0.167em} =\hspace{0.167em} -\hspace{0.167em} {\limsup }_{n\rightarrow \infty }\frac{1}{n+1}\hspace{0.167em} {\sum }_{\omega _{0}^{n}}\hspace{0.167em} \mu \left [\hspace{0.167em} \omega _{0}^{n}\hspace{0.167em} \right ]\hspace{0.167em} \log \mu \left [\hspace{0.167em} \omega _{0}^{n}\hspace{0.167em} \right ], \end{eqnarray} \tag{ 21 }$$

where the sum holds over all possible blocks ${\omega }_{0}^{n}$. Note that in the case of a Markov chain h[ μ ] also reads [16]:

$$\begin{eqnarray} &&h\left [\hspace{0.167em} \mu \hspace{0.167em} \right ]\hspace{0.167em} =\hspace{0.167em} -{\sum }_{\omega _{0}^{D}}\hspace{0.167em} \mu \left [\hspace{0.167em} \omega _{0}^{D}\hspace{0.167em} \right ]\hspace{0.167em} P\left [\hspace{0.167em} \omega (D)\hspace{0.167em} \left \vert \hspace{0.167em} \omega _{0}^{D-1}\right .\hspace{0.167em} \right ]\hspace{0.167em} \log P\left [\hspace{0.167em} \omega (D)\hspace{0.167em} \left \vert \hspace{0.167em} \omega _{0}^{D-1}\right .\hspace{0.167em} \right ], \end{eqnarray} \tag{ 22 }$$

whereas, when D = 0, h[ μ ] reduces to the (5).

As a general result from ergodic theory [54, 29, 12] and mathematical statistical physics [21], there is a unique⁶ probability distribution μ such that [54, 29, 12]:

$$\begin{eqnarray} &&\mathcal{P}\left [\hspace{0.167em} {\mathcal{H}}_{\boldsymbol{\beta }}\hspace{0.167em} \right ]={\sup }_{\nu \in {\mathcal{M}}_{\mathrm{inv}}}\left (\hspace{0.167em} h\left [\hspace{0.167em} \nu \hspace{0.167em} \right ]\hspace{0.167em} +\hspace{0.167em} \nu \left [\hspace{0.167em} {\mathcal{H}}_{\boldsymbol{\beta }}\hspace{0.167em} \right ]\hspace{0.167em} \right )=h\left [\hspace{0.167em} \mu \hspace{0.167em} \right ]\hspace{0.167em} +\hspace{0.167em} \mu \left [\hspace{0.167em} {\mathcal{H}}_{\boldsymbol{\beta }}\hspace{0.167em} \right ], \end{eqnarray} \tag{ 23 }$$

where [ _β ] is given by (13). _inv is the set of all possible time-translation invariant probabilities on the set of rasters with N neurons and $\nu \left [\hspace{0.167em} {\mathcal{H}}_{\boldsymbol{\beta }}\hspace{0.167em} \right ]={\sum }_{{\omega }_{0}^{D}}{\mathcal{H}}_{\boldsymbol{\beta }}({\omega }_{0}^{D})\hspace{0.167em} \nu ({\omega }_{0}^{D})$ is the average value of _β with respect to the probability ν.

Looking at the second equality, the variational principle (23) selects, among all possible probabilities ν, a unique one realizing the supremum. This is exactly the invariant distribution of the Markov chain and is the sought Gibbs distribution. It is clear from (23) that the topological pressure is the formal analogue to a thermodynamic potential density, where _β somewhat fixes the 'ensemble': $\nu \left [\hspace{0.167em} {\mathcal{H}}_{\boldsymbol{\beta }}\hspace{0.167em} \right ]={\mathop{\sum }\nolimits }_{k=1}^{K}{\beta }_{k}\nu \left [\hspace{0.167em} {\mathcal{O}}_{k}\hspace{0.167em} \right ]$ plays the role of βE (canonical ensemble), βE − μN (grand canonical ensemble), ... in thermodynamics  [4].

The inverse problem of finding the β_k values from the observables average measured on the data is a hard problem with no exact analytical solution. However, in the context of spatial models with pairwise interactions, the wisdom coming from statistical physics, and especially the Ising model and spin glasses, as well as from the Boltzmann machine learning community, can be used. As a consequence, in this context, several strategies have been proposed. Ackley et al [2] proposed a technique to estimate the parameters of a Boltzmann machine. This technique is effective for small networks but it is time consuming. In practice, the time necessary to learn the parameters increases exponentially with the number of units. To speed up the parameter estimation, analytical approximations of the inverse problem have been proposed, which express the parameters β_k as a nonlinear function of the correlations of the activity (see for example [68, 51, 57, 45, 53, 2, 25, 28]).

These methods do not give an exact result, but are computationally fast. We do not pretend to review all of them here, but we quote a few prominent examples. In [57], Sessak and Monasson proposed a systematic small-correlation expansion to solve the inverse Ising problem. They were able to compute couplings up to the third order in the correlations for generic magnetizations, and to the seventh order in the case of zero magnetizations. Their resulting expansion outperforms existing algorithms on the Sherrington–Kirkpatrick spin-glass model [60].

Based on a high-field expansion of the Ising thermodynamic potential, Cocco et al [14] designed an algorithm to calculate the parameters in a time polynomial with N, where the couplings are expressed as a weighted sum over the power of the correlations. They did not obtained a closed analytical expression, but their algorithm could run in a time that was polynomial in the number of neurons.

Other methods, based on Thouless–Anderson–Palmer equations [71] and linear response [28], or information geometry [68], initially proposed in the field of spin glasses, have been adapted and applied to spike train analysis (see e.g. the work done by Roudi et al [52]).

The success of these approximations depends on the dataset, and there is no a priori guarantee about their efficiency at finding the right values of the parameters. However, by getting closer to the correct solution, they can potentially speed up the convergence of the learning by starting with a seed much closer to the real solution than if taking a random starting point.

Note also that all the techniques mentioned above have been designed for the case where there is no temporal interaction (except [14, 52] which are discussed in section 2.3.5). Now, we explain how the parameter estimation can be done in the spatio-temporal models.

In the general case the parameters β_k can be determined thanks to the following properties.

• [ _β ] is a log generating function of cumulants. First:
  $$\begin{eqnarray} &&\frac{\partial \mathcal{P}\left [\hspace{0.167em} {\mathcal{H}}_{\boldsymbol{\beta }}\hspace{0.167em} \right ]}{\partial {\beta }_{k}}=\mu \left [\hspace{0.167em} {\mathcal{O}}_{k}\hspace{0.167em} \right ]. \end{eqnarray} \tag{ 24 }$$
  This is an extension of (9) to the time-dependent case.
• Second:
  $$\begin{eqnarray} &&\frac{{\partial }^{2}\mathcal{P}\left [\hspace{0.167em} {\mathcal{H}}_{\boldsymbol{\beta }}\hspace{0.167em} \right ]}{\partial {\beta }_{k}\partial {\beta }_{l}}=\frac{\partial \mu \left [\hspace{0.167em} {\mathcal{O}}_{k}\hspace{0.167em} \right ]}{\partial {\beta }_{l}}=\sum _{n=0}^{+\infty }{C}_{{\mathcal{ O}}_{k}{\mathcal{O}}_{l}}(n), \end{eqnarray} \tag{ 25 }$$
  where C_k l(n) is the correlation function between the two observables _k and _l at time n. Note that correlation functions decay exponentially fast whenever _β has finite range. So that ${\mathop{\sum }\nolimits }_{n=0}^{+\infty }{C}_{{\mathcal{O}}_{k}\hspace{0.167em} {\mathcal{O}}_{l}}(n)\lt +\infty$.Equation (25) characterizes the variation in the average value of _k when varying β_l (linear response). The corresponding matrix is a susceptibility matrix. It controls the Gaussian fluctuations of observables around their mean (central limit theorem) [54, 44, 12]. This is the generalization of (10) to the time-dependent case. As a particular case, the fluctuations of the empirical average ${\pi }_{\omega }^{(T)}\left [\hspace{0.167em} {\mathcal{O}}_{k}\hspace{0.167em} \right ]$ of _k around its mean μ[ _k ] are Gaussian with a mean-square deviation $\sqrt{\mu \left [\hspace{0.167em} {\mathcal{O}}_{k}\hspace{0.167em} \right ](1-\mu \left [\hspace{0.167em} {\mathcal{O}}_{k}\hspace{0.167em} \right ])}/\sqrt{T}$.It is clear that the structure of the linear response in the case of spatio-temporal constraints is quite a bit more complex than the case D = 0 (see equation (10)). Actually, for D = 0, all correlations C_k l(n) vanish for n > 0 (distinct times are independent).
• (_β) is a convex function of β. As a consequence, if there is a set of β values, β*, such that
  $$\begin{eqnarray} &&\frac{\partial \mathcal{P}\left [\hspace{0.167em} {\mathcal{H}}_{\boldsymbol{\beta }}\hspace{0.167em} \right ]}{\partial {\beta }_{k}^{\ast }}=\mu \left [\hspace{0.167em} {\mathcal{O}}_{k}\hspace{0.167em} \right ]={C}_{k}, \end{eqnarray} \tag{ 26 }$$
  then this set is unique. Thus, the solution of the variational problem (23) is unique.

Basically, equations (24)–(26), tell us that techniques based on free energy expansion in spatial models can be extended as well to spatio-temporal cases, where the free energy is replaced by the topological pressure. Obviously, estimating (not to speak of computing) the topological pressure can be a formidable task. Although the transfer matrix technique allows the computation of the topological pressure, the use of this method for large N is hopeless (see section 3.1). However, techniques based on periodic orbit expansion and zeta functions could be useful [44]. Additionally, cumulant expansions of the pressure, equations (24) and (25) corresponding to the two first orders, suggest that the extension of methods based on free energy expansion could be used. In addition to the works quoted above, we can also think of constraint satisfaction problems by Mézard and Mora [37] and approaches based on Bethe free energy [78]. Finally, as we checked, the properties of spatio-temporal Gibbs distributions allows one to extend the parameter estimation methods developed for the spatial case in [17, 7] to spatio-temporal distributions (to be published).

Here we shortly review alternative spatio-temporal models. We essentially refer to approaches attempting to construct a Markov chain and related invariant probability by proposing a specific form for the transition probabilities.

A prominent example is provided by the so-called linear–nonlinear (LN) models and generalized linear models (GLM) [5, 36, 43, 74, 48, 46, 3, 47]. Shortly, the idea is to model spike statistics by a point process where the instantaneous firing rate of a neuron is a nonlinear function of the past network activity, including feedbacks and interaction between neurons [63]. This model has been applied in a wide variety of experimental settings [6, 13, 70, 8, 42, 74, 46]. Typically, referring e.g. to [3], the rate r_i has the form (adapting to our notations):

$$\begin{eqnarray} &&{r}_{i}=f\left [\hspace{0.167em} {b}_{i}+{K}_{i}.x+{\sum }_{j}{H}_{i j}{r}_{j}\hspace{0.167em} \right ] \end{eqnarray} \tag{ 27 }$$

where the kernel K_i represents the ith cell's linear receptive field and x is an input. H_ij characterizes the effect of spikes emitted in the past by pre-synaptic neuron j on post-synaptic neuron i. In this approach, neurons are assumed to be conditionally independent given the past. The probability to have a given spike response to a stimulus, given the past activity of the network, reads as the product of firing rates (see e.g. equation (eq2.4) in  [3]).

In [3] the authors use several Monte Carlo approaches to learn the parameters of the model for a Bayesian decoding of the rasters. Comparing to the method presented in the previous sections, the main advantages of the GLM are: (i) the transition probability is known (postulated) from the beginning and does not require the heavy normalization imposed by potentials of the form (6); (ii) the model parameters have a neurophysiological interpretation, and their number grows at most as a power law in the number of neurons, as opposed to (6), where the parameters are delicate to interpret and whose number can become quite large, depending on the set of constraints.

Note, however, that a model of the form (27) can be written as well in the form (6): this is a straightforward consequence of the Hammersley–Clifford theorem [22]. The parameters β_k in (6) are then nonlinear functions of the parameters in (27) (see [10] for an example).

The main drawback of this approach is the assumption of conditional independence between neurons: neurons are assumed independent at time t when the past, which appears in the function H_ij in (27), is given, and the probability of a spiking pattern at time t is the product of neuron firing rates. In contrast, the maximal entropy principle does not require this assumption.

It is interesting to note that the conditional independence assumption can be rigorously justified in conductance-based integrate-and-fire models [10, 11] and the form of the function f can be explicitly found (this a sigmoid function instead of an exponential as usually postulated in GLM). This result holds true if only chemical synapses are involved (this is also implicit in the kernel form H_ij in (27) [3]), but conditional independence breaks down, for example, as soon as electric synapses (gap junctions) are involved: this can be mathematically shown in conductance-based integrate-and-fire models [15]. Note that, in this case, a large fraction of the correlations are due to dynamical interactions between neurons: as a consequence they persist even if there is no shared input.

Recently, Macke et al [34] extended the GLM model to fix the lack of instantaneous correlations between neurons in the GLM. They added a common input function that has a linear temporal dynamics. However, one of the disadvantages of this technique is that its likelihood is not unimodal, and thus computationally expensive expectation–maximization algorithms have to be used to fit parameters.

The GLM model is usually used to model both the stimulus–response dependence as well as the interaction between neurons, while the MaxEnt models usually focus on the latter (but see [72]).

To finish this subsection, we would like to quote two important works dealing with spatio-temporal events too. First, in [14] Cocco and co-workers consider the spiking activity of retinal ganglion cells with a dual approach: on one hand they consider an Ising model (and higher-order spatial terms) where they propose an inverse method based on a cluster expansion to find efficiently the coupling in the Ising model from data; on the other hand, they consider the problem of finding the parameters (synaptic couplings) in a integrate-and-fire model with noise from its spike trains. In the weak noise limit the conditional probability of a spiking pattern, given the past, is given by a least action principle. This probability is a Gibbs distribution whose normalized potential is characterized by the action computed over an optimal path. This second approach allows the characterization of spatio-temporal events. Especially it gives a very good fit of the cross correlograms.

Second, in [52], the authors consider a one step memory Markov chain where the conditional probability has a time-dependent potential of the Ising type. Adapting a Thouless–Anderson–Palmer [71] approach used formerly in the Sherrington–Kirkpatrick mean-field model of spin glasses [60] they propose an inversion algorithm to find the model parameters. As in the GLM, their model assumes conditional independence given the past (see equation (1) in [52]).

The Kullback–Leibler divergence between μ,μ^(∗) is given by:

$$\begin{eqnarray} &&d({\mu }^{(\ast )},\mu )={\limsup }_{n\rightarrow \infty }\frac{1}{n+1}{\sum }_{\omega _{0}^{n}}{\mu }^{(\ast )}\left [\hspace{0.167em} \omega _{0}^{n}\hspace{0.167em} \right ]\log \left [\hspace{0.167em} \frac{{\mu }^{(\ast )}\left [\hspace{0.167em} \omega _{0}^{n}\hspace{0.167em} \right ]}{\mu \left [\hspace{0.167em} \omega _{0}^{n}\hspace{0.167em} \right ]}\hspace{0.167em} \right ], \end{eqnarray} \tag{ 28 }$$

which provides some notion of asymmetric 'distance' between μ and μ^(∗). The KL divergence accounts for the discrepancy between the predicted probability $\mu \left [\hspace{0.167em} {\omega }_{0}^{n}\hspace{0.167em} \right ]$ and the exact probability ${\mu }^{(\ast )}\left [\hspace{0.167em} {\omega }_{0}^{n}\hspace{0.167em} \right ]$ for all blocks of range n.

This quantity is not numerically computable from (28). However, for μ a Gibbs distribution and μ^(∗) a time-translation invariant probability, the following holds:

$$\begin{eqnarray*} &&{d}_{\mathrm{KL}}({\mu }^{(\ast )},\mu )=\mathcal{P}\left [\hspace{0.167em} {\mathcal{H}}_{\boldsymbol{\beta }}\hspace{0.167em} \right ]\hspace{0.167em} -\hspace{0.167em} {\mu }^{(\ast )}\left [\hspace{0.167em} {\mathcal{H}}_{\boldsymbol{\beta }}\hspace{0.167em} \right ]\hspace{0.167em} -\hspace{0.167em} h({\mu }^{(\ast )}). \end{eqnarray*}$$

The topological pressure [ _β ] is given by (13) while μ^(∗)[ _β ] is estimated by ${\pi }_{\omega }^{(T)}\left [\hspace{0.167em} {\mathcal{H}}_{\boldsymbol{\beta }}\hspace{0.167em} \right ]={\mathop{\sum }\nolimits }_{k=1}^{K}{\beta }_{k}{\pi }_{\omega }^{(T)}\left [\hspace{0.167em} {\mathcal{O}}_{k}\hspace{0.167em} \right ]={\mathop{\sum }\nolimits }_{k=1}^{K}{\beta }_{k}{C}_{k}$.

Since μ^(∗) is unknown, h(μ^(∗)) is unknown, and can only be estimated from data, i.e. one estimates the entropy of the empirical probability, $h({\pi }_{\omega }^{(T)})$. There exist efficient methods for that. Note that the entropy of a Markov chain is readily given by equation (22), so the entropy $h({\pi }_{\omega }^{(T)})$ is obtained by replacing the exact probability P in equation (22), by the empirical probability $h({\pi }_{\omega }^{(T)})$. As T → +∞, $h({\pi }_{\omega }^{(T)})\rightarrow h({\mu }^{(\ast )})$, at exponential rate⁷. For finite T, finite size corrections exist, see e.g. Strong et al [67]. In figure 1 is plotted an example. For a potential _β with N = 5 neurons and range R = 2, containing all possible observables, we have plotted the difference between the exact probability (known from (22) and the explicit form (16), (17) of transition probabilities and invariant probability) and the approached entropy $h({\pi }_{\omega }^{(T)})$ obtained by replacing the exact probability P by the empirical probability $h({\pi }_{\omega }^{(T)})$, as a function of raster size T. We have also plotted the finite corrections method proposed by Strong et al in  [67].

Zoom In Zoom Out Reset image size

Now, if one wants to compare how two Gibbs distributions μ₁,μ₂ approximate data, one compares the divergence d_KL( μ^(∗),μ₁ ), d_KL( μ^(∗),μ₂ ), where h(μ^(∗)) is independent of the model choice. Therefore, the comparison of two models can be done without computing h(μ^(∗)).

Another criterion, easier to compute, is to compare the expected value of the observables average, μ^(∗)[_k], known from (24) to the empirical average ${\pi }_{\omega }^{(T)}\left [\hspace{0.167em} {\mathcal{O}}_{k}\hspace{0.167em} \right ]$. Error bars are expected to follow the central limit theorem where fluctuations are given by equation (25). Examples are given in figure 2. Note that the comparison of observables average is less discriminant than minimizing the Kullback–Leibler divergence, since there are infinitely many possible models matching the observables average.

Zoom In Zoom Out Reset image size

In figure 2 we show the comparison between the exact values of the observable averages (ordinate) and the estimated Monte Carlo values (abscissa). The error bars have been computed with N_seed samples. The tests were performed with N_seed = 20, N_times = 10 000 and N_flip = 100 000. In this example, N goes from 3 to 8 and R = 3. For larger values of NR the numerical computation with the transfer matrix method is no longer possible (NR = 24 corresponds to matrices of size 16 777 216 × 16 777 216).

In this section, we show how the Kullback–Leibler divergence varies as a function of N_times. In figure 3 we show the evolution of the Kullback–Leibler divergence between the real distribution and its estimation with the Monte Carlo method.

Zoom In Zoom Out Reset image size

We also consider the error

$$\begin{eqnarray} &&\mathrm{error}={\max }_{k=1\cdots K}\left \vert \hspace{0.167em} 1-\frac{\pi _{}^{(T)}\left [\hspace{0.167em} {\mathcal{O}}_{k}\hspace{0.167em} \right ]}{\mu \left [\hspace{0.167em} {\mathcal{O}}_{k}\hspace{0.167em} \right ]}\hspace{0.167em} \right \vert . \end{eqnarray} \tag{ 32 }$$

As developed in section 3.3, this quantity is expected to converge to 0 if N_times → +∞ when N_flip grows proportional to N × N_times. For finite N_times, N_flip → +∞ the error is controlled by the central limit theorem. The probability that the error on the average of observable _k is bigger than (equation (30)), behaves like exp(  − N_times × ²/σ(_k) ).

In contrast, if N_flip stays constant while N_times grows, the error is expected to first decrease to a minimum, after which it increases. This is because the number of flips is insufficient to reach the equilibrium distribution of the Monte Carlo–Hastings Markov chain. This effect is presented in figure 4. It shows the error (32) as a function of N_times for N = 3 to 7 neurons with 3 different N_flip values (1000, 10 000, 100 000).

Clearly, the number of flips N_flip should be at least more than N × N_times in order to give a chance to all spikes in the raster to be flipped at least once. A value N_flip = 10 × N × N_times seems to be enough and computationally reasonable to perform the estimations. With an N_seed = 20, we have results with a reasonable error around the mean values.

Zoom In Zoom Out Reset image size

Here we compare the CPU time for a Monte Carlo simulation and the time for a computation with the transfer matrix. We illustrate this in figure 5. We have plotted the CPU time necessary to obtain the observables average presented in figure 2, for the Monte Carlo average and for the exact average, as a function of N. The CPU time for the Monte Carlo method increases slightly more than linearly while the CPU time for the transfer matrix method increases exponentially fast: note that R = 3 here, so that Rlog2 = 2.08, close to the exponential rate found by fit: 2.24.

Zoom In Zoom Out Reset image size

We also plot in figure 6 the CPU times as a function of N_times with three N_flip values (1000, 10 000, 100 000), for 3 and 7 neurons. The CPU time increases in a linear fashion with the N_times value. The CPU time also increases linearly with the N_flip value (for the same N value). The simulations have been done on a computer with the following characteristics: 7 Intel(R) Xeon(R) 3.20 GHz processors with a 31.5 Gb RAM.

Zoom In Zoom Out Reset image size

In figure 7 we compare the CPU time increase with N_times for several N values. It shows that the CPU time increases linearly with N_times (as in figure 6). However, the slope increases with the number of neurons N.

Zoom In Zoom Out Reset image size

We fix the number of neurons N and the range R and we choose L distinct pairs (i_l,t_l), l = 1⋯L, i_l∈{ 1,...,N }, t_l∈{ 0,...,D − 1 }. To this set is associated a set of K = 2^L events _k = ( ω_i1(t₁),...,ω_il(t_l) ), k = 0,...,K − 1. For example, if L = 2, there are 4 possible events ₀ = (0,0) : neuron i₁ is not firing at time t₁ and neuron i₂ is not firing at time t₂; ₁ = (0,1) : neuron i₁ is not firing at time t₁ and neuron i₂ is firing at time t₂; and so on. It is convenient to have a label k corresponding to the binary code of the event.

We define K observables _k of range D taking binary values 0,1. For a block ${\omega }_{0}^{D}$, ${\mathcal{O}}_{k}\left [\hspace{0.167em} {\omega }_{0}^{D}\hspace{0.167em} \right ]=0$ if the event _k is not realized in the bloc ${\omega }_{0}^{D}$ and is 1 otherwise. In the example above, ${\mathcal{O}}_{0}\left [\hspace{0.167em} {\omega }_{0}^{D}\hspace{0.167em} \right ]=1$ if neuron i₁ is not firing at time t₁ and neuron i₂ is not firing at time t₂ in the block ${\omega }_{0}^{D}$. Thus, for N = 3, R = 4, (i₁,t₁) = (1,0);(i₂;t₂) = (1,1),

$$\begin{eqnarray*} {\mathcal{O}}_{0}\left [\hspace{0.167em} \left (\hspace{0.167em} \begin{array}{@{}cccc@{}} \displaystyle 0&\displaystyle 0&\displaystyle 0&\displaystyle 1\\ \displaystyle 0&\displaystyle 1&\displaystyle 0&\displaystyle 1\\ \displaystyle 0&\displaystyle 1&\displaystyle 0&\displaystyle 1 \end{array}\hspace{0.167em} \right )\hspace{0.167em} \right ]={\mathcal{O}}_{0}\left [\hspace{0.167em} \left (\hspace{0.167em} \begin{array}{@{}cccc@{}} \displaystyle 0&\displaystyle 0&\displaystyle 0&\displaystyle 1\\ \displaystyle 1&\displaystyle 1&\displaystyle 0&\displaystyle 1\\ \displaystyle 1&\displaystyle 1&\displaystyle 0&\displaystyle 1 \end{array}\hspace{0.167em} \right )\hspace{0.167em} \right ]=1, \end{eqnarray*}$$

while

$$\begin{eqnarray*} {\mathcal{O}}_{0}\left [\hspace{0.167em} \left (\hspace{0.167em} \begin{array}{@{}cccc@{}} \displaystyle 1&\displaystyle 0&\displaystyle 0&\displaystyle 1\\ \displaystyle 0&\displaystyle 1&\displaystyle 0&\displaystyle 1\\ \displaystyle 0&\displaystyle 1&\displaystyle 0&\displaystyle 1 \end{array}\hspace{0.167em} \right )\hspace{0.167em} \right ]={\mathcal{O}}_{0}\left [\hspace{0.167em} \left (\hspace{0.167em} \begin{array}{@{}cccc@{}} \displaystyle 0&\displaystyle 1&\displaystyle 0&\displaystyle 1\\ \displaystyle 1&\displaystyle 1&\displaystyle 0&\displaystyle 1\\ \displaystyle 1&\displaystyle 1&\displaystyle 0&\displaystyle 1 \end{array}\hspace{0.167em} \right )\hspace{0.167em} \right ]=0. \end{eqnarray*}$$

We finally define a potential as in (6), ${\mathcal{H}}_{\boldsymbol{\beta }}\hspace{0.167em} =\hspace{0.167em} {\mathop{\sum }\nolimits }_{k=1}^{K}{\beta }_{k}{\mathcal{O}}_{k}$.

For this type of potential, whatever ${\omega }_{0}^{D-1}$, ${\sum }_{{\omega }_{0}^{D-1}}{\mathrm{e}}^{{\mathcal{H}}_{\boldsymbol{\beta }}\left (\hspace{0.167em} {\omega }_{0}^{D}\hspace{0.167em} \right )}$ is independent of ω(D). As a consequence of the Perron–Frobenius theorem ${s}_{\boldsymbol{\beta }}={\sum }_{{\omega }_{0}^{D-1}}{\mathrm{e}}^{{\mathcal{H}}_{\boldsymbol{\beta }}\left (\hspace{0.167em} {\omega }_{0}^{D}\hspace{0.167em} \right )}$ and therefore:

$$\begin{eqnarray*} &&\mathcal{P}({\mathcal{H}}_{\boldsymbol{\beta }})=(N-L)\log (2)+\log \left [\hspace{0.167em} \sum _{k=1}^{K}{\mathrm{e}}^{{\beta }_{k}}\hspace{0.167em} \right ]. \end{eqnarray*}$$

As a consequence, from (24), giving the average of observable _k as the derivative of with respect to β_k:

$$\begin{eqnarray*} &&\mu \left [\hspace{0.167em} {\mathcal{O}}_{k}\hspace{0.167em} \right ]=\frac{{\mathrm{e}}^{{\beta }_{k}}}{\sum _{n=1}^{K}{\mathrm{e}}^{{\beta }_{n}}}. \end{eqnarray*}$$

The fluctuations of observables can also be estimated as well. They are Gaussian with a covariance matrix given by the Hessian of (see equation (25)) and the central limit theorem.

Remark An important assumption here is that observables do not depend on ω(D). This important simplification as well as the specific form of observables makes the computation of tractable. Note however that, although does not depend on ω(D) as well, the normalized potential and therefore the conditional probability $P\left [\hspace{0.167em} \omega (D)\hspace{0.167em} \left \vert \hspace{0.167em} {\omega }_{0}^{D-1}\right .\hspace{0.167em} \right ]$ depend on ω(D) thanks to the normalization factor and its dependence in the right eigenvector R of the Perron–Frobenius matrix.

Let us now use this example as a benchmark for our Monte Carlo method. We have considered a case with range R = 4 and L = 6 pairs, corresponding to 2⁶ = 64 terms in the potential. We have analysed the convergence rate as a function of N, the number of neurons, and N_times, the time length of the Monte Carlo raster. The number of flips, N_flip is fixed to 10 × N × N_times, so that, on average, each spike of the Monte Carlo raster is flipped ten times in one trial.

In figure 8 (left) we have shown the relative error (equation (32)) as a function of N_times for several values of N. The empirical average ${\pi }_{}^{(T)}\left [\hspace{0.167em} {\mathcal{O}}_{k}\hspace{0.167em} \right ]$ is computed on 10 Monte Carlo rasters. That's why we don't write the index ω in the empirical probability. We stopped the simulation when the error is lower than 5%. As expected from the central limit theorem (CLT), the error decreases as a power law $C{N}_{\mathrm{times}}^{-1/2}$, where the constant C has been obtained by fit.

In figure 8 (right) we have drawn the CPU time as a function of N_times for several values of N. It increases linearly with N_times, with a coefficient depending on N. The simulation is relatively fast: it takes 150 mins for N = 60, N_times = 8000 (with 10 Monte Carlo trials) on a 7-processor machine (each of these processors has the following specifications: Intel(R) Xeon(R) CPU 2.27 GHz, 1.55 MB of cache memory size) with a 17.72 GB RAM.

Zoom In Zoom Out Reset image size