Thermodynamic Formalism in Neuronal Dynamics and Spike Train Statistics

Abstract

The Thermodynamic Formalism provides a rigorous mathematical framework for studying quantitative and qualitative aspects of dynamical systems. At its core, there is a variational principle that corresponds, in its simplest form, to the Maximum Entropy principle. It is used as a statistical inference procedure to represent, by specific probability measures (Gibbs measures), the collective behaviour of complex systems. This framework has found applications in different domains of science. In particular, it has been fruitful and influential in neurosciences. In this article, we review how the Thermodynamic Formalism can be exploited in the field of theoretical neuroscience, as a conceptual and operational tool, in order to link the dynamics of interacting neurons and the statistics of action potentials from either experimental data or mathematical models. We comment on perspectives and open problems in theoretical neuroscience that could be addressed within this formalism.

1. Introduction

Initiated by Boltzmann [1,2], the goal of statistical physics was to establish a link between the microscopic mechanical description of interacting particles in a gas or a fluid, and the macroscopic description that is provided by thermodynamics [3,4].

Although this program is, even nowadays, far from being completed [1,5], the work of Boltzmann and his successors opened new avenues of research, not only in physics but also in mathematics. Especially the term “ergodic”, which was coined by Boltzmann [1], inaugurated an important branch of mathematics that provides a rigorous link between the description of dynamical systems in terms of their trajectories and the description in terms of statistics of orbits and, more generally, between dynamical systems theory and probability theory. At the core of the ergodic theory, there is a set of "natural" dynamically invariant probability measures in the phase space, somewhat generalising the Liouville distribution for conservative systems with strong analogies with Gibbs distributions in statistical physics [6,7]. This strong connection, in particular, gave birth to the so-called Thermodynamic Formalism.

The introduction of Thermodynamic Formalism that occurred in the 1970s was primarily due to Yakov Sinai, David Ruelle, and Rufus Bowen [8,9,10]. The development of Thermodynamic Formalism initially served to derive rigorous criteria characterising the existence and uniqueness of the Gibbs states in the infinite volume limit. Although Gibbs states and equilibrium states (see Section 2.1.1) are naturally defined in finite volume systems, the extension to infinite volume (“thermodynamic limit”) is far from straightforward. Indeed, it does not follow from the Carathéodory or Kolmogorov extension theorems [11,12], that the equilibrium states of the infinite volume define a measure, as there is no way to express the marginals associated to an infinite-volume Gibbs measure without making explicit reference to the measure itself [12]. When considering conditional probabilities, rather than marginals, Dobrushin, Lanford, and Ruelle led to a different consistency condition that affords for the building of infinite volume Gibbs measures [7].

In the context of dynamical systems, Sinai, Ruelle, and Bowen were able to connect the theory of hyperbolic (Anosov) dynamical systems to results in statistical mechanics. Indeed, Sinai found an unexpected link between the equilibrium statistical mechanics of spin systems and the ergodic theory of Anosov systems by a codification using Markov partitions (see Section 2.1.1 for details). This idea was later extended for a much more general class of hyperbolic systems [10,13,14]. While the Thermodynamic Formalism started as a branch of rigorous statistical mechanics, nowadays it is viewed from different communities as a branch of dynamical systems or stochastic processes.

There have been a few attempts to use Thermodynamic Formalism in different ways other than as a natural mathematical foundation of statistical mechanics, for example, studying population dynamics [15,16], self-organised criticality [17], the relative abundance of amino acids across diverse proteomes [18], analyse the difference between introns and exons in genetic sequences [19,20], coding sequence density estimation in genomes [21], and statistics of spike trains in neural systems [22,23,24,25,26,27,28], which is the main topic of this review. Neuronal networks are biological systems whose components, such as neurons, synapses, ionic channels, ..., are ruled by the laws of physics and are written in terms of differential equations; hence, are dynamical systems. On the other hand, because of their large dimensionality, it is natural to attempt to characterise neuronal networks using methods inspired by statistical physics, for example, mean-field methods [29,30,31,32], density methods [33] or Fokker-Planck equations [34]. Most neurons produce short electrical pulses, called action potentials or spikes, and it is widely believed that the collective spike trains emitted by neuronal networks encode information about the underlying dynamics and response to stimuli [35,36,37]. Thus, researchers have devoted a lot of effort to understanding the correlations structure in the statistics of spike trains [38,39,40,41,42,43,44,45].

Because spikes can be represented as binary variables, it is natural to adapt methods and concepts from statistical physics, and more specifically the statistical physics of spin systems, to analyse spike trains statistics. There have been many successful attempts in this direction. All of the approaches we know about are based on variational principles. The most direct connection from statistical physics to spike train statistics is done via the maximum entropy principle which has attracted a lot of attention during the past years [38,39,43,46,47] (see [47] for a physicists-oriented review). Unfortunately, most of these articles are limited to the original form of an Ising spin-glass potential (pairwise interactions with random couplings) or variants of it with higher order interactions [40,41,42], where successive times are independent, thereby neglecting the time correlations and causality one may expect from a network of neurons with interactions (exceptions can be found in [48,49,50]). We focus on this approach and its extension to causal networks in the present review as it is natural to link with the Thermodynamic Formalism. Another approach, which actually appeared earlier in mathematical neuroscience, is the dynamic mean-field theory that takes into account dynamics and time correlations. In this approach, originated in quantum field theory and Martin–Siggia–Rose formalism (Statistical dynamics of classical systems [51]), the variational principle is expressed via the minimisation of an effective action containing the equations for the dynamics. It was introduced in the field of theoretical neuroscience by Sompolinsky who initially applied it to the study of spin-glasses dynamics [52,53,54], before analysing neuronal networks dynamics [29]. Here, the effective action can be summarised, in the limit when the number of neurons tends to infinity, by dynamic mean-field equations depending on neuron’s firing rates and pairwise time-correlations. Thus, here, the information about spike statistics is contained in the two first moments (Gaussian limit). This approach has inspired a lot of important works, where temporal correlations are taken into account, see, for example [55,56,57], or the recent work by M. Helias and collaborators (i.e., the review on dynamic mean-field theory in [58] and the link to from dynamic mean-field theory to large deviations from Ben-Arous and Guionnet [59] and [60]). Another trend attempts to relate neuronal dynamics to spike statistics, expressed via a maximum-likelihood approach [61] (we apologise if we have forgotten other approaches that we ignore). To our taste, the most achieved work in this direction is Amari’s information geometry [62] making a beautiful link between probability measures (e.g., exponential, like Gibbs measures) and differential geometry.

In this review, we show how Thermodynamic Formalism can be used as an alternative way to study the link between the neuronal dynamics and spike statistics, not only from experimental data, but also from mathematical models of neuronal networks, properly handling causal and memory dependent interactions between neurons and their spikes. It is also based on a variational principle (and, actually, the four methods that are discussed in this introduction certainly have a common “hat”, the large deviations theory [63]), although we depart from this principle at some point where we discuss extensions of this approach to non-stationary situations through a linear response formula. Additionally, Thermodynamic Formalism provides a rigorous mathematical framework to study phase transitions that may illuminate the current discussions regarding signatures of criticality observed in many examples of neuroscience [64,65,66,67,68,69,70,71,72] and, especially, in Gibbs distributions inferred while using the maximum entropy principle from experimental data [73,74,75]. The aim of this review is twofold. On the one hand, to bridge the gap between mathematicians working in the field of Thermodynamic Formalism and scientists interested in characterising spike statistics, especially those that apply the maximum entropy principle, which is placed here in a broader context. On the other hand, to show new perspectives in the field of mathematical neuroscience related to Thermodynamic Formalism, including phase transitions.

The rest of the article is organised, as follows. In Section 2, we introduce the main tools and ideas of Thermodynamic Formalism that we use later. In Section 3 and Section 4, we present the uses of this formalism in neuroscience, for spike train statistics. In Section 5, we end with a discussion and present perspectives for future works.

2. Mathematical Setting of Thermodynamic Formalism

Our goal in this section is to present the basic tools and ideas of Thermodynamic Formalism to the unfamiliar reader (detailed accounts of this subject can be found in [9,10,76,77,78]).

2.1. General Properties

In order to study a dynamical system from the perspective of the Thermodynamic Formalism, we need, first of all, to describe the set of elements that need to be either finite or countable. Thus, continuous particle characteristics such as position, speed, or momentum do not enter in the setting we analyse here, unless one “coarse grain“ the phase space with a finite or countable partition. Discrete particle characteristics, like spin, or symbols attached to specific features of a dynamical system, constitute the set of symbols denoted by A also referred to as the alphabet. Let $A^M$ be the set of blocks of M symbols of the alphabet, that is, sequences of the form $x_0{x}_1\cdots x_{M-1}$, where $x_i\in A$, $i=0\cdots M-1$ and $A^{\mathbb{N}}$, is the set of right-infinite sequences of the form $x=x_0{x}_1\cdots$, with $x_i\in A$ for all $i\in \mathbb{N}$. One may also consider the bi-infinite sequences $A^{\mathbb{Z}}$, but we will mostly stick to $A^{\mathbb{N}}$ in the sequel.

One can equip the space $A^{\mathbb{N}}$ with a distance. We associate to $\theta \in (0,1)$ the distance ${\mathrm{d}}_{\theta }$ such that ${\mathrm{d}}_{\theta }(x,y)={\theta }^m$, where m is the largest non-negative integer, such that $x_i=y_i$ for every $0\le i<m$. In particular, if $x=y$ then $m=\infty$ and ${\mathrm{d}}_{\theta }(x,x)=0$, and, by convention, if $x_0\ne y_0$ then $m=0$. When considering this distance, the metric space $(A^{\mathbb{N}},{\mathrm{d}}_{\theta })$ is compact [10].

In order to consider time order, we introduce the time evolution in the form of a left-shift, $\sigma :A^{\mathbb{N}}\to A^{\mathbb{N}}$, defined by ${\left(\sigma x\right)}_i=x_{i+1}$, for all $i\in \mathbb{N}$ and any $x\in A^{\mathbb{N}}$, where the index i refers to time. That is, the i-th symbol of $\sigma x$ (the image of x under $\sigma$), corresponds to the $(i+1)$-th symbol of the original x, for all i.

Let us define continuous functions $f:A^{\mathbb{N}}\to \mathbb{R}$. We introduce the modulus of continuity of f:

$${\mathrm{var}}_k\left(f\right):=sup\left\{\right|f\left(x\right)-f\left(y\right)|:x_i=y_i,\mathrm{for}i=0,\dots ,k-1.\},$$

characterising the maximal variation of f on the set of infinite sequences agreeing on their first k symbols. The function f is continuous if ${\mathrm{var}}_k\left(f\right)\to 0$ when $k\to \infty$, and the continuity is exponential if ${\mathrm{var}}_k\left(f\right)\to 0$ exponentially fast. A function $f:A^{\mathbb{N}}\to \mathbb{R}$ has range R if $f\left(x\right)\equiv f(x_0,\cdots ,x_{R-1})$, i.e., f only depends on the first R coordinates $x_0,\cdots ,x_{R-1}$. Functions with finite range are obviously continuous. The notion of continuity, and, especially, exponential continuity is essential when studying the Thermodynamic Formalism of functions with infinite range.

A useful concept for this formalism is that of the Markov partitions. Let $\Omega$ be a general (continuous) state space and a general dynamics (a flow or a discrete time mapping) $\mathcal{G}$ on $\Omega$. Consider a finite covering of $\Omega$, made by sets called rectangles, $\mathcal{R}=\{R_1,\dots ,R_l\}$ with the property that, for every pair $i\ne j$, $\mathrm{int}\left(R_i\right)\cap \mathrm{int}\left(R_j\right)=\varnothing$, and such that the closure of the image set $G\left(R_i\right)$ equals a union of closures of sets in $\mathcal{R}$. There is a specific construction for different dynamical systems (for details we refer the reader to [10]). This construction allows one to make a natural coding for the state space, where the trajectories of the dynamical systems correspond to the sequences of symbols of the labels identifying the sets $R_i$. Rectangles allow for us to partition a continuous phase space into a discrete finite partition. For hyperbolic dynamical systems, each point in the phase space admits local stable or unstable manifolds with a non zero diameter. The edges of rectangles in this case are made of these local stable and unstable manifolds and the corresponding partition is called a “Markov partition”, because the image of a rectangle under the dynamics can be represented by a Markov transition matrix [79].

2.1.1. Gibbs Measures

We now define the measurable sets on $A^{\mathbb{N}}$. Let us denote the sequence $x_0\dots ,x_{n-1}$ by $x_{0,n-1}$. The set $\left[x_{0,n-1}\right]:=\{y\in A^{\mathbb{N}}:y_0=x_0,\dots ,y_{n-1}=x_{n-1}\}$ is called a cylinder. Cylinders define a Borel sigma algebra $\mathcal{B}$ on $A^{\mathbb{N}}$ (see, for instance, Chapter 1, in [80]). We are interested in probability distributions on $(A^{\mathbb{N}},\mathcal{B})$ also referred as (macro) states in physics. In the sequel we will skip the prefix “macro” and deal with “states” as probability measures on $(A^{\mathbb{N}},\mathcal{B})$. There is a special class of states $\nu$, such that, for any measurable set $B\subset A^{\mathbb{N}}$, $\nu \left({\sigma }^{-1}\left(B\right)\right)=\nu \left(B\right)$, where ${\sigma }^{-1}\left(B\right)$ stands for the set of all the pre-images of $\sigma$ of elements of B. These distributions are the set of shift-invariant probability measures denoted by ${\mathcal{M}}_{\sigma }\left(A^{\mathbb{N}}\right)$. In general, it is possible to consider systems where there exists some forbidden transitions between symbols. In that case, we need to consider a subset of A invariant under the dynamics, i.e., $X\subset A^{\mathbb{N}}$, such that $\sigma \left(X\right)=X$.

The analogy with the spin systems in statistical mechanics at the root of the terminology “Thermodynamic Formalism”, goes, as follows. Let $\varphi :A^{\mathbb{N}}\to \mathbb{R}$ be a continuous function, also called “energy” or “potential”. Two important examples are finite range potentials, where $\varphi \left(x\right)\equiv \varphi \left(x_{0,R-1}\right)$, or exponentially continuous potentials with infinite range. Subsequently, the “energy” of a configuration of n sites based on the sequence $x\in X$, is given by the “Birkhoff sums”:

$$S_n\varphi \left(x\right):=\sum _{i=0}^{n-1}\varphi \circ {\left(\sigma x\right)}_i.$$

We define the measures that assign probability to the cylinder sets based on the potential function, with the so called “Gibbs property”. There exist constants $C>1$ and $F\left(\varphi \right)$, such that, for all $x\in X$ and for all $n\ge 1$:

$$C^{-1}\le \frac{{\mu }_{\varphi }\left(\left[x_{0,n-1}\right]\right)}{exp({\sum }_{i=0}^{n-1}\varphi \circ {\left(\sigma x\right)}_i-nF\left(\varphi \right))}\le C.$$

(1)

The measures satisfying the condition (1) are called “Gibbs measures”. The quantity:

$$F\left(\varphi \right)=\underset{n\to \infty }{lim}\frac{1}{n}log\sum _{x_0\cdots x_{n-1}}{e}^{S_n\varphi \left(y\right)},$$

(2)

for all $y\in \left[x_{0,n-1}\right]$ [10], is called the “pressure” (or free energy) of the potential $\varphi$. Observe that it does not depend on the measure ${\mu }_{\varphi }$ itself, only on the potential. Thus, two Gibbs measures that are associated with the same potential have the same pressure.

Given a continuous potential $\varphi$, such that ${\mathrm{var}}_k\left(\varphi \right)\le b{\theta }^k$ for some constants $0<\theta <1$ and, $b>0$, there is a unique shift-invariant probability measure satisfying the Gibbs property (1), [10]. Furthermore, under the same assumption for the potential, the associated Gibbs measure is mixing i.e., ${lim}_{n\to \infty }\mu \left(A\cap {\sigma }^{-n}B\right)=\mu \left(A\right)\mu \left(B\right)$, for any measurable set $A,B$ and thus ergodic (A measure $\mu$ is said to be ergodic for the dynamics $\mathcal{G}$ if for any measurable $\mathcal{G}$-invariant set A, (${\mathcal{G}}^{-1}\left(A\right)=A$), its measure is either $\mu \left(A\right)=0$ or $\mu \left(A\right)=1$. See, for instance, Chapter 3 of [79] for a detailed introduction.) (see [10,80] for definitions and details). Moreover, two continuous potentials $\varphi$ and $\psi$ are called cohomologous with respect to the shift $\sigma$ (and denoted by $\varphi \sim \psi$), if there is a continuous function u and a constant K, such that:

$$\varphi =\psi +u-u\circ \sigma -K.$$

Cohomologous potentials have the same Gibbs measure, ${\mu }_{\varphi }={\mu }_{\psi }$.

2.1.2. Entropy and the Variational Principle

The Shannon entropy of a probability distribution quantifies the average level of “uncertainty” in its possible outcomes.

$$h\left(p\right):=-\sum _{x\in A}p\left(x\right)logp\left(x\right).$$

(3)

More generally, let $\nu$ be a shift-invariant probability measure on X, we introduce the block entropy:

$$H_n\left(\nu \right):=-\sum _{x_{0,n-1}\in A^n}\nu \left(\left[x_{0,n-1}\right]\right)log\nu \left(\left[x_{0,n-1}\right]\right).$$

For finite alphabets and because $\nu$ is shift-invariant, the following limit exists,

$$h\left(\nu \right):=\underset{n\to \infty }{lim}\frac{1}{n}{H}_n\left(\nu \right),$$

and is called the Kolmogorov–Sinai entropy or simply entropy of $\nu$ [76,81]. The key is that the sequence $H_n$ is sub-additive and the possible values lie in the compact set $[0,log|A\left|\right]$ (see [76]).

Another important quantity is the Kullback–Leibler (KL) divergence, which quantifies the difference between probability distributions. Consider two probability distributions p and q on the same probability space X, the KL divergence is:

$$D\left(p\right|\left|q\right):=\sum _{x\in X}p\left(x\right)log\frac{p\left(x\right)}{q\left(x\right)}.$$

This divergence is finite whenever p is absolutely continuous with respect to q, and it is only zero if $p=q$.

Following the analogy with systems in statistical mechanics, an equilibrium state is defined as the measure that satisfies the so-called variational principle, namely:

$$sup\left\{h\left(\nu \right)+\int \varphi \mathrm{d}\nu :\nu \in {\mathcal{M}}_{\sigma }\left(X\right)\right\}=h\left({\mu }_{\varphi }\right)+\int \varphi \mathrm{d}{\mu }_{\varphi }=F\left(\varphi \right),$$

(4)

where $\int \varphi \mathrm{d}\nu$ represents the expected value of $\varphi$ with respect to the measure $\nu$ (it would be noted ${\mathbb{E}}_{\nu }\left(\varphi \right)$ in a probabilistic context). Let us comment on this result. The first equality establishes that, among all shift-invariant probability measures, there is a unique one, ${\mu }_{\varphi }$, which maximises $h\left(\nu \right)+\int \varphi \mathrm{d}\nu$. If $\int \varphi \mathrm{d}\nu$ is fixed to some value, this corresponds to maximising the entropy under constraints. This is the Maximum Entropy Principle, but (4) is more general. The second equality states that the maximum, $h\left({\mu }_{\varphi }\right)+\int \varphi \mathrm{d}{\mu }_{\varphi }$ is exactly the pressure defined in Equation (2). The equation $F\left(\varphi \right)=h\left({\mu }_{\varphi }\right)+\int \varphi \mathrm{d}{\mu }_{\varphi }$ establishes a link with thermodynamics, as $F\left(\varphi \right)$ is equivalent to the free energy (In contrast to thermodynamics where free energy or pressure refer to different thermodynamic ensembles, we will not make such distinction here as it is irrelevant. Also, note that we do not keep the minus sign and take the Boltzmann constant equal to one. Thus, (4) coincides with the principle of free energy minimisation in statistical mechanics, but it corresponds to a maximising measure in our formalism because of a sign convention. Gibbs measures associated to potentials satisfying that ${\mathrm{var}}_k\left(\varphi \right)\le b{\theta }^k$, as stated above, are equilibrium states. For more general potentials, the notion of Gibbs states and equilibrium states may not coincide.

2.2. Observables and Fluctuations of Their Time Averages

An observable f is a function $f:X\to \mathbb{R}$, such that $\left|f\right|<\infty$ ($|·|$ denote the absolute value), and which is time-translation invariant, i.e., for any time i, $f\left(x_{0,R-1}\right)=f\left(x_{i,i+R-1}\right)$ whenever $x_{0,R-1}=x_{i,i+R-1}$.

An important observable, considered in the following sections, is the potential:

$$U_{\lambda }\left(x\right)=\sum _k{\lambda }_k{f}_k\left(x\right)k\in \{1,...,K\},$$

(5)

that is a linear combination of K observables $f_k$, where the parameters ${\lambda }_k$’s $\in \mathbb{R}$. Here, we want to make a few important remarks. First, this decomposition is similar to the definition of potentials in thermodynamics, which are written in terms of a linear combination of intensive variables (e.g., temperature, chemical potential, and so on) and extensive functions of the configurations of the system (physical energy, number of particles). In order to emphasise this analogy, we use the symbol U for the potential, instead of $\varphi$. The function $U_{\lambda }$ parametrically depends on the ${\lambda }_k$’s; hence, we use the lower index $\lambda$, to denote $U_{\lambda }$ as well as the corresponding Gibbs measure, denoted here ${\mu }_{\lambda }$.

Now, we denote $A_n\left(f\right)$ the empirical average of the observable f in a sample of size n, that is,

$$A_{n}f\left(x\right)=\frac{1}{n}\sum _{i=0}^{n-1}f\left(\sigma x\right).$$

This quantity is important for empirical purposes. If $\mu$ is ergodic, the convergence of the empirical average to the actual expected value is $\mu$-almost-sure, that is $\stackrel{}{A_n\left(f\right)\stackrel{\mathrm{a}.\mathrm{s}.}{\to }\int fd\mu }$ as $n\to \infty$. For an example of application of this theorem in the context of spike train statistics, see (Section 2.6).

2.3. Correlations

Let us consider a pair of observables $f,g\in L^2\left(\mu \right)$, (square integrable functions with respect to $\mu$). We define their correlation at time n by

$$C_{f,g}\left(n\right):=\int f·g\circ {\sigma }^n\mathrm{d}\mu -\int f\mathrm{d}\mu \int g\mathrm{d}\mu .$$

One also might be interested in the auto-correlation (or auto-covariance) of f at time n,

$$C_f\left(n\right):=\int f·f\circ {\sigma }^n\mathrm{d}\mu -{\left(\int f\mathrm{d}\mu \right)}^2.$$

Observe that these quantities might decay fast. The mixing property implies that the correlations go to zero with $n\to +\infty$.

2.3.1. Properties of the Pressure

From the pressure (2), important statistical information regarding the system can be obtained, in particular, correlations. A distinguished case corresponds to the potentials of the form (5). When the corresponding pressure is differentiable to any order, taking the successive derivatives of the pressure with respect to their conjugate parameters gives the average values of the observables, their correlations, and their high-order cumulants with respect to the Gibbs measure. That is, in general:

$$\frac{{\partial }^{n}F\left(U_{\lambda }\right)}{\partial {\lambda }_k^n}={\kappa }_n\mathrm{for} \mathrm{all}k\in \{1,...,K\},$$

(6)

where ${\kappa }_n$ is the cumulant of order n with respect to ${\mu }_{\lambda }$. In particular, ${\kappa }_1$ is the mean of $f_k$, ${\kappa }_2$ is its variance, ${\kappa }_3$ the skewness, and ${\kappa }_4$ the kurtosis. Partial derivatives with respect to pairs of parameters can also be considered [82]:

$$\begin{array}{cc}\hfill \frac{{\partial }^{2}F\left(U_{\lambda }\right)}{\partial {\lambda }_k\partial {\lambda }_j}& =C_{f_k,f_j}\left(0\right)+\sum _{n=1}^{\infty }{C}_{f_k,f_j}\left(n\right)+\sum _{n=1}^{\infty }{C}_{f_j,f_k}\left(n\right)\hfill \\ \hfill & =\frac{\partial \mu \left[f_k\right]}{\partial {\lambda }_j}=\frac{\partial \mu \left[f_j\right]}{\partial {\lambda }_i}.\hfill \end{array}$$

(7)

Observe, we differentiate $C_{f_k,f_j}$ from $C_{f_j,f_k}$ as the dynamics may be irreversible in time. For an example of application of these formulas in the context of spike train statistics, see (Section 2.6).

Remark 1.

The last two equations are fundamental. They establish a link between the variations in the average of the observable $f_k$, when slightly varying the parameter ${\lambda }_j$ and the sum of time correlations between the pair of observables $f_j,f_k$. This result is known, in statistical physics as the fluctuation-dissipation theorem [83,84]. It relates, for example, in a ferromagnetic model, to the variation of the magnetisation of the spin k to the variations of the local magnetic field $h_k$, via the magnetic susceptibility, which is the second derivative of the free energy. This is also the context of the linear response theory, which quantifies how a small perturbation of a parameter affect the average values of observables in terms of the unperturbed measure. Equation (7) can also be used to extend results of Thermodynamic Formalism to non-stationary situations, as we discuss in Section 4.5.

Now, in the classical formulation in statistical physics and the Maximum Entropy models, only the correlation $C_{f_k,f_j}\left(0\right)$ appears in (7), because successive times are independent (correlations $C_{f_k,f_j}\left(n\right),C_{f_j,f_k}\left(n\right),n>0$ vanish). When handling memory (thus, potentials with range $R>1$), there is an infinite sum (series) of correlations appearing in the linear response. This infinite sum converges whenever correlations are decaying sufficiently fast with time n (exponentially). In contrast, when they do not converge fast enough (e.g., power-law with a small exponent), the series diverges, leading to a divergence of the second derivative of the pressure, corresponding to a second-order phase transition. Here, follow the Ehrenfest classification of phase transitions [85]. There is a phase transition of order k if the pressure (free energy) is $C^{k-1}$ but not $C^k$. The known examples are first-order, second-order, or infinite order (Kosterlitz-Thouless) phase transitions. Note that phase transitions in memory-less models can also happen if the instantaneous correlation function $C_{f_k,f_j}\left(0\right)$ diverges, e.g., when the number of degrees of freedom in the system (number of spins, neurons) tends to infinity. Therefore, in the present setting, second-order phase transitions can arise either when the number of degrees of freedom tends to infinity, or (not exclusive), when time correlations decay slowly.

Thus, Equation (7) connects the second derivative of the pressure, variations in static average of observables, and dynamical correlations. In the next sections, we discuss theorems that relate the dynamical evolution to a criteria ensuring that time-correlations are exponentially decaying, preventing the possibility of second order phase transitions for systems with a finite number of degrees of freedom.

2.3.2. Ruelle-Perron-Frobenius Operator

Let $C\left(X\right)$ be the set of continuous functions on X. Consider the potential $\varphi$ and a continuous function $f\in C\left(X\right)$, so one can define a bounded linear operator associated to $\varphi$ (transfer operator), called the Ruelle–Perron–Frobenius (RPF), as follows (There is a close analogy between this operator and the propagator in quantum field theory or the Koopman operator in classical dynamics [5]. All of these operators characterise how measures or observables evolve ruled by the dynamics.):

$${\mathcal{L}}_{\varphi }f\left(x\right)=\sum _{y\in {\sigma }^{-1}x}{e}^{\varphi \left(y\right)}f\left(y\right).$$

(8)

The spectral properties of (8) yields information to characterise the pressure and study the ergodic properties of the system, in particular, the rate of decay of their correlation functions [80]. For instance, if 1 is a simple eigenvalue and the modulus of each of the other eigenvalues is smaller than one, this is equivalent to be mixing [80]. When the potential considered is of finite range, then the transfer operator corresponds to a matrix and the whole formalism is equivalent to Markov chains defined on finite alphabets. A potential $\varphi$ is called normalised if ${\mathcal{L}}_{\varphi }\left(1\right)=1$. The log of a normalised potential of range $R+1$, corresponds to the transition probabilities of a Markov chain with memory depth R. Moreover, in this case $F\left(\varphi \right)=0$. For Lipschitz observables in the finite dimensional case, the Perron–Frobenius theorem assures a unique eigenvector associated to the maximal eigenvalue, from which the unique invariant measure (Markov) is obtained. This measure has mixing properties, exponential decay of correlations, central limit theorem, and a large deviations principle (see Section 2.3.4). When the operator acts on an infinite dimensional space (such as the space of continuous functions), then the spectrum of a bounded linear operator $\mathcal{L}$ is given by the set $\mathrm{spec}\left(\mathcal{L}\right)=\{\lambda \in \mathbb{C}:\mathrm{such} \mathrm{that} (\lambda I-\mathcal{L}\left)\mathrm{has} \mathrm{no} \mathrm{bounded} \mathrm{inverse}\right\}$, this set may contain points that are not necessarily eigenvalues (see, for instance, [86]). In this case, the strategy is to find a proper subspace where the spectrum of $\mathcal{L}$ has a finite number of such complex numbers whose norm is the spectral radius, say $\rho$, and the rest of the spectrum has norm strictly less than $\rho$ (spectral gap). In this scenario, it is known that there is exponential decay of correlations for a sufficiently regular class of observables (such as Lipschitz), and the central limit theorem holds. In the absence of the spectral gap, then one has sub-exponential decay of correlations, which breaks down the central limit theorem, and the phase transition phenomenon appears (for further details and precise definitions, see [80] and the references therein).

Note that, given a potential $\psi$, one can explicitly find a normalised potential $\varphi$ cohomologous to $\psi$, as follows,

$$\varphi :=\psi +logR-logR\circ \sigma -log\rho ,$$

(9)

where R is the right eigenvector (real and positive) associated to the unique maximum eigenvalue $\rho$ that is associated to ${\mathcal{L}}_{\psi }$.

Remark 2.

Note that the normalisation of the potential ψ does not require a partition function. In fact, as discussed below, the classical normalisation by a partition function is a particular case of (9), holding for memory-less potentials that does not generalise to range $R>1$ potentials.

2.3.3. Time Averages and Central Limit Theorem

We have seen that if the measure $\mu$ is ergodic the time averages $A_n\left(f\right)$ converge $\mu$-almost surely to the expected value $\int f\mathrm{d}\mu$. Now, we can ask about fluctuations around the expected value. The observable f satisfies the central limit theorem (CLT), with respect to $(\sigma ,\mu )$ if:

$$\frac{A_n\left(f\right)-n\int f\mathrm{d}\mu }{\sqrt{n}}\stackrel{\mathrm{law}}{⟶}{\mathcal{N}}_{0,{\sigma }_f^2}$$

(10)

where ${\mathcal{N}}_{0,{\sigma }_f^2}$ is the Gaussian distribution with zero mean and covariance ${\sigma }_f^2$, which is given by the following expression involving temporal correlations:

$${\sigma }_{f,g}^2=C_{f,g}\left(n\right)\left(0\right)+\sum _{n=1}^{\infty }{C}_{f,g}\left(n\right)+\sum _{n=1}^{\infty }{C}_{g,f}\left(n\right).$$

that is a particular case of (7). We illustrate an application of this theorem in the context of spike train statistics later in Section 2.6.

Strong properties of convergence and exponential decay of correlations are ensured for Hölder continuous potentials in finite dimension. These properties are associated with the spectral gap property and they do not (necessarily) hold for less regular potentials or in non-compact spaces [78,80].

2.3.4. Large Deviations

The central limit theorem describes small fluctuations in the limit when n goes to infinity. Rare events that are exponentially small are the object of study of the large deviations theory.

An empirical average $A_n\left(f\right)$ satisfies a large deviation principle (LDP) with rate function $I_f$, if the following limit exists:

$$I_f\left(s\right):=-\underset{n\to \infty }{lim}\frac{1}{n}log\mathbb{P}\left(\left\{A_n\left(f\right)>s\right\}\right),$$

(11)

for $s\in \mathbb{R}$. The above condition for large n implies that $\mathbb{P}\left(\{A_n\left(f\right)>s\}\right) \approx e^{-n{I}_f\left(s\right)}$. In particular, if $s>\int f\mathrm{d}\mu$, then $\mathbb{P}\left(\{A_n\left(f\right)>s\}\right)$ should tend to zero as n increases. The rate function tells us precisely how fast this probability goes to zero. Computing the rate function from Equation (11), may be a laborious task. The Gärtner-Ellis theorem provides a way to compute $I_f$ more easily [63]. To this end, let us introduce the scaled cumulant generating function (SCGF) associated with the observable f, by

$${\Gamma }_f\left(k\right)=:\underset{n\to \infty }{lim}\frac{1}{n}log\int e^{nk{A}_n\left(f\right)}\mathrm{d}\mu k\in \mathbb{R},$$

whenever the limit exists. The name comes from the fact that the n-th cumulant of f can be obtained by successive differentiation operations over ${\Gamma }_f\left(k\right)$ with respect to k, and then evaluating the result at $k=0$. If ${\Gamma }_f$ is differentiable, then the Gärtner-Ellis theorem ensures that the average $A_n\left(f\right)$ satisfies a LDP with rate function given by the Legendre transform of ${\Gamma }_f$, that is

$$I_f\left(s\right)=\underset{k\in \mathbb{R}}{max}\{ks-{\Gamma }_f\left(k\right)\}.$$

Therefore, one can study the large deviations of empirical averages $A_n\left(f\right)$ by first computing their SCGF, characterise its differentiability, and then find the Legendre transform. We compute this function in the context of spike train statistics later in Section 2.6.

If ${\Gamma }_f\left(k\right)$ is differentiable then $I_f\left(s\right)$ is convex [87], thus has a unique global minimum $s^{*}$ such that $I_f\left(s^{*}\right)=0$, then $I_f^{\prime }\left(s^{*}\right)=0$. Assume that $I_f\left(s\right)$ admits a Taylor expansion around $s^{*}$, then for s close to $s^{*}$,

$$I_f\left(s\right)=I_f\left(s^{*}\right)+I_f^{\prime }\left(s^{*}\right)\left(s-s^{*}\right)+\frac{I_f^{″}\left(s^{*}\right){\left(s-s^{*}\right)}^2}{2}+O{\left(s-s^{*}\right)}^3.$$

Because $I_f\left(s^{*}\right)=0$ and $I_f^{\prime }\left(s^{*}\right)=0$, for large values of n we obtain from (11)

$$\begin{array}{cc}\hfill \mathbb{P}\left(\left\{A_n\left(f\right)>s\right\}\right)& \approx e^{-n{I}_f\left(s\right)}\hfill \\ & \approx e^{-n\left(\frac{I_f^{″}\left(s^{*}\right){\left(s-s^{*}\right)}^2}{2}\right)}\hfill \end{array}$$

Therefore, the small deviations of $A_t\left(f\right)$ around $s^{*}$ are Gaussian with variance $1/n{I}_f^{″}\left(s^{*}\right)$. In this way, the LDP can be regarded as an extension of the CLT.

The large deviation principle plays an important role in statistical mechanics, in particular in spin glass dynamics [59]. A large deviation principle can be used in order to relate entropy and free energy (here pressure) through a Legendre transform and to explain why variational principles arise in statistical mechanics [63,88]. As mentioned in the introduction, large deviations is the common theoretical principle linking dynamic mean field theory, maximum entropy principle, maximum likelihood, and Thermodynamic Formalism, although this link has not been studied in detail, to our best knowledge.

2.4. Potentials of Range One

A specific case where the variational principle (4) holds, is when the potential has the form (5). Then, equilibrium states are probability distributions ${\mu }_{\lambda }$, that maximise the entropy (3), under the constraints of expected values of K observables ${\mathbb{E}}_{{\mu }_{\lambda }}\left(f_k\right):={\sum }_x{f}_k\left(x\right){\mu }_{\lambda }\left(x\right)=C_k$ for $k=1,\cdots ,K$. This problem can be solved by introducing a Lagrange multipliers ${\lambda }_k$ in the potential (5):

$$F\left(U_{\lambda }\right):=\underset{p}{max}\{H\left(p\right)+{\mathbb{E}}_p\left(U_{\lambda }\right)\}=H\left({\mu }_{\lambda }\right)+{\mathbb{E}}_{{\mu }_{\lambda }}\left(U_{\lambda }\right).$$

(12)

There exists a unique maximum entropy distribution ${\mu }_{\lambda }$ (equilibrium state) satisfying the constraints. It turns out that the maximising distribution can be explicitly found for range one potentials and the distribution satisfies the Gibbs property (1), which, in this particular case, reduces to,

$${\mu }_{\lambda }\left(x\right)=exp\left(-F\left(U_{\lambda }\right)+U_{\lambda }\left(x\right)\right)=\frac{exp\left(U_{\lambda }\left(x\right)\right)}{Z},$$

(13)

for all $x\in X$. Equation (13) is obtained by considering $F\left(U_{\lambda }\right)=logZ$, where Z is the “partition function”. From Equation (12), the expression for the entropy (3) and the Jensen inequality, one can obtain the formula for the pressure in this case:

$$F\left(U_{\lambda }\right)=log\sum _{x\in X}{e}^{U_{\lambda }\left(x\right)}.$$

The constrained problem can be uniquely solved, because the map $\lambda \mapsto {\mathbb{E}}_{{\mu }_{\lambda }}\left(U\right)$ maps the real line monotonically onto the interval $(minU,maxU)$ [76].

For range one potentials, the measure of a block becomes a product distribution, as given by:

$${\mu }_{\lambda }\left(\left[x_{0,n-1}\right]\right)=\prod _{i=0}^{n-1}\frac{exp\left(U_{\lambda }\left(x_i\right)\right)}{Z}.$$

(14)

As the index n corresponds to time, having an interaction depending on one single coordinate implies that configurations at distinct times are independent.

2.5. Finite Range Potentials

Equation (9) can be used to find the unique Markov measure that is associated with a finite range potential. As an example, consider a potential U of range two representing the pairwise interactions in a graph with incidence matrix I. The entries $I(y,x)=1$ represent the allowed transitions between symbols $y\to x$ and $I(y,x)=0$ the forbidden. We introduce the finite $\mid A\mid \times \mid A\mid$ transfer matrix ${\mathcal{L}}_U$, which corresponds to the RPF “operator” (8) restricted to a finite space.

$${\mathcal{L}}_U(y,x)=I(y,x)e^{U\left(y\right)},y,x\in A,y\in {\sigma }^{-1}x$$

(15)

As anticipated in Section 2.3.2, calling $\rho$ the unique maximal positive eigenvalue of ${\mathcal{L}}_U$ guaranteed by the Perron–Frobenius theorem, and $R\left(x\right)$ and $L\left(x\right)$ the x-th entry in the right and left eigenvectors associated with $\rho$, respectively, we define a normalized potential $\varphi (y,x)=U\left(y\right)+logR\left(y\right)-logR\left(x\right)-log\rho$, such that the matrix

$$P(y,x)=I(y,x)e^{\varphi (y,x)}=\frac{I(y,x)R\left(y\right)e^{U\left(y\right)}}{\rho R\left(x\right)}$$

(16)

is stochastic, i.e., ${\sum }_{x}P(y,x)=1$, and represents the transition probabilities of a Markov chain $P(y\to x)=P(x\mid y)$. The invariant measure p associated to the matrix P satisfying $pP=p$ is

$$p\left(x\right)=\frac{R\left(x\right)L\left(x\right)}{\langle R,L\rangle },$$

(17)

where $\langle R,L\rangle ={\sum }_{x}R\left(x\right)L\left(x\right)$. Note that normalisation is done without defining a partition function. The Markov measure $\mu (p,P)$ of a block is given by $\mu \left[x_{0,n}\right] = p\left(x_0\right)P\left(x_1,x_2\right)\cdots P\left(x_{n-1},x_n\right)$ for $x_k\in A,k=0,..,n.$ Here, we have a nice way to show that this measure satisfies the Gibbs property while using the Markov property $\mu \left[x_{1,n}|x_0\right]=e^{{\sum }_{k=1}^{n}logP\left(x_{k-1},x_k\right)}$ where we see that the conditioning upon the first time is similar to left boundary conditions in statistical physics.

It follows from (9) and (17) that $\mu \left[x_{0,n}\right]$ obeys the variational principle and satisfies Equation (1), where $F\left(U\right)=log\rho .$ The Gibbs measure $\mu \left[x_{0,n}\right]$ gives an exponential weight to each cylinder set, depending on the “energy” depending on smaller blocks.

2.6. Example

To illustrate the maximum entropy principle and the statistical analysis that can be performed while using tools and ideas from Thermodynamic Formalism, we include here a toy example. Consider the state space of all the binary blocks of size $2\times 2$ and one step transitions between them. We associate to each block en integer (23), and index a matrix using this representation of blocks we built the RPF matrix (15). There are allowed and forbidden transitions as explained in Section 3.3.1 (see Figure 3). Assume that we obtain from data (T samples) the empirical average value of the observable $A_T(x_0^1·x_1^2)=0.1$ and $A_T(x_1^1·x_0^2)=0.4$ and we want to find the maximum entropy Markov chain compatible with these constraints. Using Equations (6), (15) and (16), and, we obtain the maximum entropy Markov chain, defined by the following Markov transition matrix:

$$\begin{array}{cc}\begin{array}{c}\\ 0\\ 1\\ 2\\ 3\\ 4\\ 5\\ 6\\ 7\\ 8\\ 9\\ 10\\ 11\\ 12\\ 13\\ 14\\ 15\end{array}& \begin{array}{c}\begin{array}{ccccccccccccccccccccccccccccccc}0 & & 1 & & 2 & & 3 & & 4 & & 5 & & 6 & & 7 & & 8 & & 9 & & 10 & & 11 & & 12 & & 13 & & 14 & & 15\end{array}\\ \left(\begin{array}{cccccccccccccccc}0.16& 0.04& 0.64& 0.16& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0.64& 0.16& 0.16& 0.04& 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0& 0.04& 0.16& 0.16& 0.64& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0.16& 0.64& 0.04& 0.16\\ 0.16& 0.04& 0.64& 0.16& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0.64& 0.16& 0.16& 0.04& 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0& 0.04& 0.16& 0.16& 0.64& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0.16& 0.64& 0.04& 0.16\\ 0.16& 0.04& 0.64& 0.16& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0.64& 0.16& 0.16& 0.04& 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0& 0.04& 0.16& 0.16& 0.64& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0.16& 0.64& 0.04& 0.16\\ 0.16& 0.04& 0.64& 0.16& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0.64& 0.16& 0.16& 0.04& 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0& 0.04& 0.16& 0.16& 0.64& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0.16& 0.64& 0.04& 0.16\end{array}\right)\end{array}\end{array}$$

From this Markov transition matrix, we can compute the fluctuations that are associated to each observable either using numerical simulations or analytically. We illustrate in Figure 1, the limit theorems and fluctuations introduced in Section 2 applied to this example.

The entropy maximisation for this toy example can be explicitly solved, and the simulations can also be performed directly from the transition matrix. However, large scale networks require sophisticated Montecarlo sampling methods to fit maximum entropy models that include non-synchronous interactions [89]. In the first column of Figure 1, we directly sample from the Markov transition matrix for different sample sizes and average the empirical frequency of both observables considered in the toy example. In the second column we plot the fitted Gaussian distributions of the empirical averages for different sample sizes. The third row correspond to the large deviations rate function. As explained in Section 2.3.4, the second derivative at the minimum of $I_f$ characterise the Gaussian fluctuations around the expected value of f. The last column represents the auto-correlations (7).

2.7. Systems with Infinite Range Potentials, Chains with Infinite Memory and Gibbs Distributions

In this section, we somewhat depart from the strict setting of Thermodynamic Formalism, switching to the perspective of Markov chains and their extension to infinite memory. Although Thermodynamic Formalism allows for one to consider infinite memory (infinite range potentials) the advantage of the approach presented here is to allow considering non stationary dynamics, i.e., escape from the variational principle (4) constrained by the entropy definition, which requires stationarity.

A general class of stochastic processes to deal with infinite memory are called Chains with complete connections [90,91]. These chains define non-markovian processes. However, Markovian approximations are possible and useful [92]. This section follows closely from [90].

Definition 1.

A system of transition probabilities is a family ${\left\{P_n\right\}}_{n\in \mathbb{Z}}$ of functions with $P_n(·\mid ·):A\times A_{-\infty ,n-1}\to [0,1]$, such that the following conditions hold for every $n\in \mathbb{Z}$:

(a) For every $x_n\in A$, the function $P_n\left(x_n\mid ·\right)$ is measurable with respect to the filtration ${\mathcal{F}}_{\le n-1}$.

(b) For every $x_{-\infty ,n-1}\in A_{-\infty ,n-1}$,

$$\sum _{x_n\in A}{P}_n\left(x_n\mid x_{-\infty ,n-1}\right)=1.$$

Definition 2.

A probability measure μ in $\mathcal{P}(A_{-\infty ,n},\mathcal{F})$ is consistent with a system of transition probabilities ${\left\{P_n\right\}}_{n\in \mathbb{Z}}$ if:

$$\int h\left(x_{-\infty ,n}\right)\mu \left(dx\right)=\int \sum _{x_n\in A}h\left(x_{-\infty ,n-1}{x}_n\right)P_n\left(x_n\mid x_{-\infty ,n-1}\right)\mu \left(dx\right).$$

for all $n\in \mathbb{Z}$ and all ${\mathcal{F}}_{\le n}$-measurable functions h. The probability measure μ, when it exists, is called a chain with complete connections consistent with the system of transition probabilities ${\left\{P_n\right\}}_{n\in \mathbb{Z}}.$ It is possible that multiple measures are consistent with the same system of transition probabilities.

We now give conditions ensuring the existence and uniqueness of a probability measure consistent with the system of transition probabilities [90].

Theorem 1.

A system of continuous transition probabilities (${\mathrm{var}}_m\left[P_n\left(x_n\mid ·\right)\right]\to 0$ as $m\to +\infty$) on a compact space has at least one probability measure consistent with it.

Definition 3.

A system of transition probabilities is non-null, if, for all $n\in \mathbb{Z}$ and all $x_{-\infty }^n\in A_{-\infty ,n}$:

$$P\left[x_n\left|x_{-\infty ,n-1}\right.\right]>0$$

Definition 4.

A normalized potential has bounded squared variations if, for all $n\in \mathbb{Z}$ and all $x_{-\infty ,n}\in A_{-\infty ,n}$:

$$\sum _{k\ge 0}{var}_k^2\left(logP\left[x_n\left|x_{-\infty ,n-1}\right.\right]\right)<+\infty .$$

There exists a unique probability measure consistent with the system of transition probabilities if these are non-null and the associated normalised potential has bounded squared variations [90].

There is a mathematically well-founded correspondence between chains with complete connections and Gibbs distributions presented up to now [7,90,93]. Let us now discuss the formal analogy.

Define $\varphi \left(n,x\right) : \mathbb{Z}\times A\to \mathbb{R}$ by:

$$\varphi \left(n,x\right) \equiv logP\left[x_n\left|x_{-\infty ,n-1}\right.\right],$$

(18)

and:

$$\Phi (m,n,x) = \sum _{r=m}^n\varphi \left(r,x\right).$$

(19)

Then:

$$P\left[x_{m,n}\left|x_{-\infty ,m-1}\right.\right] = e^{\Phi (m,n,x)} = e^{{\sum }_{r=m}^n\varphi \left(r,x\right)},$$

and:

$$\mu \left[x_{m,n}\right] = {\int }_{A_{-\infty ,m-1}}{e}^{\Phi (m,n,x)}\mu \left(dx\right).$$

These last equations highlight the connection between chains with complete connections and Gibbs distributions in statistical physics. Indeed, the conditional probability $P\left[x_{m,n}\left|x_{-\infty ,m-1}\right.\right]$ has a “Gibbs” form where $\varphi$ acts as an “energy” [90]. The correspondence is obtained considering “time” as a one-dimensional “lattice” and the “boundary conditions” as the past of the stochastic process. In contrast to statistical physics, there is no need to define a partition function (the potential is defined via transition probabilities, and is thus normalised).

While, for chains with complete connections defined through transition probabilities, the present is conditioned upon the past, Gibbs distributions, in general, also allow for conditioning “upon the future”. More generally, Gibbs distributions in statistical physics extend to probability distributions on ${\mathbb{Z}}^d$ where the probability to observe a certain configuration of spins in a restricted region of space is constrained by the configuration at the boundaries of this region. Therefore, they are defined in terms of specifications [7,93], which determine finite-volume conditional probabilities when the exterior of the volume is known. In one spatial dimension ($d=1$), identifying $\mathbb{Z}$ with a time axis, this corresponds to conditioning both in the past and in the future. In contrast, families of transition probabilities with an exponential continuity rate define the so-called left-interval specifications (LIS) [90,94]. This leads to nonequivalent notions of “Gibbsianness” [95].

In contrast to the potentials studied up to now, the potential (19) is defined from transition probabilities (18), which are not necessarily time-translation invariant. This is the reason why the potential is noted $\varphi (n,x)$, as it depends explicitly on time n and the configuration x. This case is closer to the setting where potential or energy is not necessarily invariant when moving along a lattice in statistical physics, therefore not constrained by the stationarity assumption made up to now. As we discuss in the next section, this is quite helpful in the study of neuronal network dynamics.

3. Thermodynamic Formalism in Neuroscience

In this section, we make the connection between Thermodynamic Formalism and spiking neuronal dynamics. From the standpoint of mathematics, there are at least two ways to consider spiking neuronal networks. First, they can be considered as biological objects whose activity can be experimentally recorded while using Multi-Electrode Arrays (MEA), often requiring sophisticated mathematical methods and algorithms for data analysis [96,97]. Second, neuronal networks are characterised by dynamical models, more or less derived from biophysics [35,36].

Here, we begin considering the statistical analysis of spike trains recorded from neuronal networks. For this case, Thermodynamic Formalism provides a powerful and insightful method to analyse the spatio-temporal statistics from experimental spike trains. We briefly mention that this formalism has afforded us to develop algorithm for spike train analysis [49,89,98] leading to the software PRANAS [99] freely available at https://team.inria.fr/biovision/pranas-software/, although we do not develop along these lines in this paper. We focus then on a specific question. When dealing with a model of spiking neurons, how much of the intrinsic dynamics of neurons, their interaction via synapses, and the influence of stimuli, constrain the collective spatio-temporal spike statistics?

Neurons are (nonlinear) entities that evolve in a concerted way (as they interact via synapses) and responding to external stimuli. The theoretical analysis of this high dimensional systems can be made thanks to mathematical methods (dynamical systems, bifurcations theory, stochastic processes, partial differential equations) or theoretical physics (statistical physics, nonlinear physics). Here, one might be interested in what Thermodynamic Formalism can contribute when considering neuronal models dynamics. In this spirit, we consider two models, the Integrate and Fire model and the Galves Löcherbach model. Most of our presentation focuses on stationary situations that are characterised by equilibrium states. Nevertheless, we consider the extension of Thermodynamic Formalism to non-stationary situations.

At the end of the section, we address a couple of open questions.

• What is the natural alphabet for spiking neuron dynamics? As we shall see, although the binary representation of spikes is a good candidate, it is too naive, as the relevant alphabet can be constructed on time blocks of spikes. A subsidiary question is about the size (time depth) of these blocks.
• Under which conditions can Thermodynamic Formalism machinery be faithfully applied to a spiking neuronal network model?
• What are the limits when the main theorems of Thermodynamic Formalism can and cannot be applied and what are the consequences for neuronal dynamics and spike statistics?

3.1. Statistics of Spike Trains and Gibbs Distribution

The human brain is composed of about a hundred billion neurons that mostly communicate among themselves together using sequences of spikes that are binary events (Sub-threshold oscillations also play an important role [100] and in organs like the retina, where most neurons do not spike.) Although the action potentials can vary in duration, amplitude, and shape, depending, e.g., on the type of neuron, they have a stereotyped shape, so that they can be considered as identical events. The main physiological reason for spike occurrence is that they can propagate information on different scales in the nervous system (centimeters to one meter) essentially without attenuation (active conduction as opposed to passive, Ohmic, conduction). However, from a contemporary point of view, spikes are also considered as events constituting “bits” of information. In this paradigm, it is tempting to consider spike trains as objects containing a “neural code” [37], i.e., a language that neurons use to communicate and that one could decipher. This terminology should not be considered literally, because, as opposed to computer codes, spike trains have a wide variability (e.g., the repetition of the same stimulus, even under controlled experimental conditions does not induce the same sequence of spikes as a response). In addition, nothing guarantees that there is only one code. When considered from the perspective of Thermodynamic Formalism, the notion of neural codes can have several meanings. (1) Spike trains constitute a symbolic coding of voltage dynamics (which depends on neuronal interactions and stimuli); and, (2) the way how neuronal dynamical systems (especially, spike trains) are affected by stimuli, provides a way for downstream networks to infer the stimulus (e.g., the retina encodes a visual scene in spike trains which are decoded by the visual cortex, capable of reconstructing a representation of the visual scene). Here, we essentially want to address the following questions: how to use Thermodynamic Formalism to fit experimental (or numerically generated) spike train data and which Gibbs distribution is produced by a network of neurons whose dynamics is known.

It is useful to consider spikes as instantaneous events (while the duration is about 1ms) and identify the maximum in the action potential course as the “time of the spike” [101]. This implicitly assumes that one considers dynamics on time scales larger than one millisecond. The binary representation is obtained using a window of a constant “binning size” (of order 10–20 ms) over the continuous time course of membrane potentials and count how many spikes there are per neuron within each time bin. Two or more spikes may occur within the same time bin, in that case, the convention is to consider these events equivalent to just one spike. This procedure [102] transforms experimental data into sequences of binary patterns (see Figure 2) leading to the following symbolic description.

Denoting the discrete time index by n, the spike-state of neuron k is denoted by $x_n^k\in \{0,1\}$, depending on whether the k-th neuron emits a spike during the n-th time bin or not. A spike pattern is the spike-state of all the neurons in a network of N neurons at a given time bin, and it is denoted by $x_n:={\left[x_n^k\right]}_{k=1}^N$. A spike block denoted by $x_{n,r}:=x_n{x}_{n+1}\cdots x_r$ is a sequence of spike patterns. The length of the spike block $x_{n,r}$ is $r-n+1$. A spike train denoted by x is the spike block representing the whole sequence of spike patterns. We consider spike trains of finite and infinite length. The set of all possible spike blocks of length R in a network of N neurons is denoted by $A_R^N$.

Thus, in comparison to the previous section, and especially Section 2.1, symbols here are spike blocks of length R. The alphabet, previously denoted A, is denoted here $A_R^N$, making explicit the dependence on the number of neurons, N, and the block depth R. It is important to make this dependence explicit, as we consider, later in this review, the effect of increasing N and R.

3.2. Conditional Probabilities for Spike Trains

The probability that a biological neuron, embedded in a network, emits a spike in a given time bin depends on the history of all variables determining the evolution of the neural network (voltage, conductances, concentrations of ions, neurotransmitters, etc.). Most of these variables are not experimentally accessible. Even if they were, there would be no hope of predicting, from this huge amount of information, the statistics of spikes. Dealing with neuronal models, the situation is simpler as there are fewer variables to control and their dynamics are known explicitly. However, even in this case, it is generally not possible to access spike statistics from dynamics. A simplification is to consider that the probability of a spike pattern only depends on the spikes emitted in the past by the network. This way, one can ignore the hidden dynamics of inaccessible variables and compute the statistics from what can be measured. Still, characterising the probability of a spike pattern given the history of the system, is generally out of reach (with, at least, two exceptions described in the next section).

The idea is to characterise the spike train statistics through a family of transition probabilities of the form:

$$P(x_{n+1}\mid x_{n-R,n})\equiv e^{\varphi \left(x_{n-R,n+1}\right)}>0.$$

(20)

where R is the memory of the spike sequence, i.e., the time horizon on which the present depends on the past. Having these transition probabilities and an initial condition (or initial distribution), one can define a Markov chain (or a chain with infinite memory if $R\to +\infty$). It is possible, for some models, to explicitly write these probabilities. In Equation (20), we have assumed that all transition probabilities are strictly positive. This is necessary to ensure the uniqueness of the corresponding Gibbs distribution. Subsequently, one can associate to (20) a range R potential $\varphi \left(x_{n-R,n}\right) > -\infty$. On experimental grounds, the problem is to estimate these probabilities from data. Since there are $2^{NR}$ possible spike blocks, for N and/or R big (e.g., $NR>20$) it becomes rapidly impossible to estimate these transition probabilities from experimental data while using a frequentist approach, as most of these transitions do not even occur within the finite experimental sample.

However, one can try to guess the form of these transition probabilities. One possibility is to start with an ad-hoc form, capturing the main features of neuronal dynamics. A canonical example is the Generalised Linear Model (GLM) [103], where the transition probabilities take the form:

$$P_n(x_{n+1}^i\mid x_{n-R,n})=f\left(b^i+\sum _{j\in {\mathcal{B}}_i}\left(H_{ij}\ast r^j\right)\left(n\right)\right)$$

(21)

where f is a non-linear function. The term $b^i$ is a constant fixing the baseline firing rate of neuron i. $H_{ij}$ is the memory kernel, ∗ is the convolution product, and $r^j\left(n\right)$ is the spike train of neuron j before time n. In this case, the memory kernel considers the spikes between $n-R$ and n, but R can go arbitrarily to the past. Here we do not consider the influence of a time-dependent external stimulus.

As we show in Section 4.3 this form can be established for discrete-time Integrate and Fire models. Equation (21) gives an Ansatz for the marginal probability that neuron i spikes at time $n+1$, given the history of the network. Equation (20), provides the joint probability of having the spike pattern at time $n+1$, and is obtained assuming that neurons are conditionally independent. This can be justified if one assumes that, due to synaptic transmission and delays, neurons do not have the time to interact within one time bin. This means that time bins must not be too large and/or synapses must not be too fast (like gap junctions [104]).

Remark 3.

The GLM, instead of describing conditional probabilities, characterises the spike rate or conditional intensity of an auto-regressive Poisson process.

A second approach is based on the variational principle (14), maximising entropy under constraints. Both of the approaches can be addressed from the perspective of Thermodynamic Formalism.

3.3. The Hammersley–Clifford Theorem

In our representation spikes take a binary value 0 or 1. Thus, any potential of range R is a function taking a finite set of values. A general theorem from Hammersley and Clifford [105,106] states that any range-R observable, in particular, the potential $\varphi \left(x_{n-R,n}\right)$, can be written in the form

$$\varphi \left(x_{n-R,n}\right)=\sum _l{\varphi }_l{m}_l\left(x_{n-R,n}\right),$$

(22)

where the coefficients ${\varphi }_l$ correspond to the decomposition of $\varphi$ in the space of finite range R-observables. This is analogous to Equation (5), with two main differences. First, the linear combination in (5) is used as an example making a link with Thermodynamics and the Maximum Entropy principle. Here, the decomposition (22) is a systematic expansion of any potential of range R defined over spike sequences. Second, in contrast to (5), the observables, denoted $f_k$ in (5) only consider finitely many values.

Equation (22) is a linear decomposition on a basis of eigenfunctions referred to from now on as monomials [107]. They have the form:

$$m_l\left(x_{n-R,n}\right)=\prod _{k=1}^d{x}_{i_k}^{n_k}.$$

where $n_k=1,\cdots ,N$ is a neuron index, and $i_k=n-R,\cdots ,n$ a time index. Thus, $m_l\left(x_{n-R,n}\right)=1$ if and only if, in the spike sequence $x_{n-R,n}$ neuron $n_k$ spikes at time $i_k$ for all $k=1,\cdots ,d$. Otherwise, $m_l\left(x_{n-R,n}\right)=0$. The number d is the degree of the monomial; degree one monomials have the form $x_{i_1}^{n_1}$, taking the value 1 if and only if neuron $n_1$ spikes at time $i_1$. Degree two monomials have the form $x_{i_1}^{n_1}{x}_{i_2}^{n_2}$, taking the value 1 if and only if neuron $n_1$ spikes at time $i_1$ and neuron $n_2$ spikes at time $i_2$, and so on. Thus, monomials provide a notion of spike interactions, similar to spins interactions in magnetic systems. For example, monomials of degree two correspond to pairwise interactions, like, e.g., in an Ising model. In contrast to the Ising model, the interactions that are considered here may involve time delays between spikes.

There are $2^{NR}$ monomials for N neurons and a given range R. One can index them by an integer l in one-to-one correspondence with the set of pairs $(i_k,n_k)$. The advantage of the monomial representation is that it focuses on spike events, which is natural for spiking neuronal dynamics. Thus, the Hammersley Clifford decomposition gives a canonical way to write any range R potentials as a linear combination of monomials of maximum degree R. This includes the GLM potential which can be embedded in the same framework [104].

As emphasised above, the Hammersley Clifford decomposition is analogous to the expression of thermodynamic potentials as a sum of products of an intensive quantity (e.g., temperature) with an extensive one (e.g., the energy). Depending on the physical constraints of the problem, one defines a thermodynamic ensemble, where the average value of extensive quantities (energy, number of particles, volume, magnetisation) is prescribed. Whereas, the first principles allow for guessing the form of the potential in thermodynamics, there is no such recipe in neuroscience. Moreover, one cannot use the complete expansion (22) on practical grounds, simply because large degree monomials have a vanishing empirical probability. More precisely, the average value of a monomial of degree d decays exponentially fast with d. This leads to two problems. (i) How to determine (from data) the constraints which are necessary to correctly characterise the spike train statistics? (ii) Are there constraints that are equivalent? (i) can be addressed in the context of information geometry [108] while (ii) can be approached using cohomology [107]. We do not further develop these aspects here, but rather refer the reader to the cited articles.

The variational principle (4) (or its finite version (12)) provides a systematic way of inferring Gibbs distributions from empirical average values of spike interactions (monomials). We make the construction explicit in the next subsections.

3.3.1. Finite Memory, Markov Chains and Gibbs Distributions

We now explicitly show how to build a Gibbs measure from a finite set of experimental averages as constraints of the maximum entropy variational problem. We assume that these constraints involve events (monomials) over a memory depth R. We build the corresponding Markov chain while using the material of Section 2.3.2. We associate to each spike block $x_{n,n+R-1}$ an integer $w_n$,

$$w_n=\sum _{r=0}^{R-1}\sum _{k=1}^N{2}^{k-1+Nr}{x}_{n+r}^k,$$

(23)

we write $w_n\sim x_{n,n+R-1}$. In this way a sequence of spike patterns (spike block) can be encoded as sequences of integers, that define the alphabet. Next, we define the incidence matrix ("grammar") between symbols of the alphabet. Not all transitions between symbols are legal or allowed. A transition between the two symbols denoted by $w_n\to w_{n+1}$ or $w_n,w_{n+1}$ is legal if the corresponding blocks overlap according to this pattern, i.e., they have the block $x_{n+1}\cdots x_{n+R-1}$ in common (see Figure 3).

This defines an incidence matrix $I(w^{\prime },w)=1$ if the transition between symbols $w^{\prime }$ and w is legal and 0 otherwise. This incidence matrix defines the grammar of allowed and forbidden words or sequences of symbols. From this incidence matrix, we define the Perron–Frobenius transfer matrix ${\mathcal{L}}_{\psi }$ in the same way, as in (15). In order to obtain the unique Markov transition matrix of maximum entropy, we follow the procedure that is explained in Section 2.3.2.

Thus, for a given choice of monomials, we associate a potential of the form (22), where the ${\lambda }_l$s that do not correspond to a chosen monomial are set to 0. Subsequently, one computes the empirical average of the chosen monomials from data. From the Perron–Frobenius theorem, there is a unique Gibbs measure, the Markov measure ${\mu }_{\lambda }$ of the Markov chain, which solves the variational problem (14), giving a statistical model of data, minimizing the KL divergence between the empirical measure and ${\mu }_{\lambda }$. Here, $\lambda$ is the set of parameters ${\lambda }_l$ that achieve the variational principle. These parameters can be numerically computed, either by using the explicit form of the measure (17) [49] or by while using the MonteCarlo methods [89]. The software PRANAS allows for the handling of spike train statistics by numerically computing the Gibbs distribution, solving (14) for up to 100 neurons [99].

3.3.2. Spectral Gap and Thermodynamic Limit

For a potential of finite range R and a finite number of neurons provided ${\lambda }_l>-\infty$, for all l in (22), the Perron–Frobenius theorem guarantees the uniqueness of the Gibbs measure ${\mu }_{\lambda }$ solving the variational principle (14). Moreover, the pressure being a real analytic function for finite range potentials, is infinitely differentiable with respect to parameters and there is an exponential decay of correlations with respect to time. This last aspect is due to the gap in the spectrum of the transfer matrix (15). These properties may not hold if either $R\to +\infty$ or, if $N\to \infty$ (corresponding to a thermodynamic limit) where the potential may lose regularity. Here, one has to consider Thermodynamic Formalism in infinite dimension, on the functional space of continuous functions. The case $R\to +\infty$ corresponds to a potential with infinite range associated, in our case, to spike statistics with infinite memory. This is discussed in the next section and in Section 4, where we show that neuronal models can have such an unbounded memory. More generally, the limits $N\to \infty$ or $R\to +\infty$ can induce important effects, such as phase transitions, which will be commented upon further in the discussion section.

4. Spiking Neuronal Network Models: The Leaky Integrate-and-Fire and Beyond

There are many models for neuronal dynamics both at the level of individual neurons and neuronal networks [35,36,109]. Here, we consider a canonical example of such a model. The first model proposed to the scientific community was introduced by Lapicque in 1907 [110]. The main interest was to make a nice link between dynamics, coding, and spikes, paving the way to use Thermodynamic Formalism in order to analyse the spike train statistics.

4.1. Dynamics and Spikes

A fundamental equation in neuronal membrane potential dynamics is the conservation of electric charge, written in its most canonical form, as follows:

$$C\frac{dV}{dt}=-\sum _X{g}_X(V-V_X)+I\left(t\right),$$

where C is the membrane capacitance of the neuron and V is its membrane potential. The sum ${\sum }_X$ holds on ionic currents of the form $i_X=-g_X(V-V_X)$ involving specific ionic channels permeable to specific ions (e.g., $N{a}^{+}$, $K^{+}$, $C{l}^{-}$, $C{a}^{2+}$). Here, we include the neuron’s intrinsic currents’ (e.g., sodium and potassium currents triggering a spike [111]) and synaptic currents [109]. The conductance, $g_X$, of channels of type X depends, in general, non-linearly on activation variables, themselves dependent on the voltage. $V_X$ is the Nernst reversal potential, i.e., the value of the membrane potential at which the current $i_X$ reverses its direction. Finally, $I\left(t\right)$ is an external current that can mimic, e.g., an injection by an electrode or an external stimulus.

In its simplest form, for a single isolated neuron, this equation takes the form:

$$C\frac{dV}{dt}=-\frac{1}{R}V+I\left(t\right),$$

(24)

where R is the membrane resistance and the term $g=\frac{1}{R}$ corresponds to a unique passive conductance. In this case, we consider only a leak current where the leak reversal potential is set to 0. This is the equation of an RC circuit, which is quite simple, but quite far from a real neuron, as this equation does not even produce spikes. To circumvent this problem, one introduces a threshold, $\theta$, such that Equation (24) holds whenever $V\left(t\right)<\theta$ (sub-threshold dynamics), corresponding to the “ntegrate” phase. In contrast, for all times $t_{\left(r\right)}$ such that $V\left(t_{\left(r\right)}\right)=\theta$, two effects take place: (i) The membrane potential of the neuron is reset instantaneously to a rest value, here 0, without a loss of generality; (ii) a spike is recorded at times $t_{\left(r\right)}$ called “spike times”. This is the “fire” phase (see Figure 2B, bottom).

While this is a simple artificial way to generate spikes, there is a huge price to pay on mathematical grounds because the threshold introduces a singularity set in the phase space where the dynamic is not differentiable. We develop this aspect below.

The generalisation of (24) to a network of N neurons is straightforward. Adding the contribution of synaptic currents, $I_{\left(syn\right)}^k\left(t\right)$, building the network interactions, we obtain:

$$C^k\frac{d{V}^k}{dt}+\frac{1}{R}{V}^k=I_{\left(syn\right)}^k\left(t\right)+S^k\left(t\right)+{\sigma }_B{\xi }^k\left(t\right),\mathrm{if}{V}^k\left(t\right)<\theta ,$$

(25)

where $C^k$ is the membrane capacitance of neuron k, $V^k$, its membrane potential. The resistance R is assumed to be the same for all neurons. In the synaptic current,

$$I_{\left(syn\right)}^k\left(t\right)=\sum _{j,r}{W}_{kj}\alpha \left(t-t_{\left(r\right)}^j\right),$$

(26)

the parameter $W_{kj}$ represents the synaptic strength (“weight”) from the pre-synaptic neuron j to the post-synaptic neuron k (see Figure 2B). Synaptic weights can be negative (inhibition) or positive (excitation). By convention, $W_{kj}=0$ if there is no connection from j to k. This way, the sum in (26) holds for $j=1\cdots N$.

The function $\alpha$ represents the time profile of the postsynaptic current induced by a pre-synaptic spike [112]. It has been experimentally observed that the tail of this function is exponential. On mathematical grounds this is essential. The sum in Equation (26) considers all the spike times $t_{\left(r\right)}^j$ emitted by all the pre-synaptic neurons j before time t. When considering the asymptotic regime $t\to -\infty$ (to get rid of initial conditions) this sum may contain an infinite number of terms. Thus, in order to ensure the sumability of this series one needs $\alpha$ to decay sufficiently fast (here exponentially fast).

Equation (25) holds in the sub-threshold regime. The term $S^k\left(t\right)$ represents an external stimulus, and ${\xi }^k\left(t\right)$ is white noise whose amplitude is modulated by ${\sigma }_B$. When the membrane potential of neuron k reaches the firing threshold at a firing time $t_{\left(r\right)}^k$, for some r, i.e., $V^k\left(t_{\left(r\right)}^k\right)\ge \theta$, the neuron k fires an action potential and its membrane potential is reset to a fixed reset value instantaneously (see Figure 2B).

While Equations (25) and (26) look rather simple, the right hand side of Equation (25) depends, via (26) on a possibly uncountable set of events (the spike times) corresponding to a possible infinite history of the voltage dynamics of the network. In this sense, these equations do not represent a classical dynamical system, where the knowledge of the variables at a given time allows one to compute the variables’ value at a future time by integrating the flow. Here, we require knowledge of the spike history, back to the last time where the neuron was reset to make the integration. This history can go quite far into the past, with a dependence decaying like the tail of the function $\alpha$.

In order to circumvent these problems, we need to first get rid of the fact that spike times belong a priori to an uncountable set. There are two alternatives. The first one, briefly explored in this section, consists of discretising time as done e.g., in [113]. This leads to important results related with Thermodynamic Formalism (see [23,114] for details). The second alternative, to keep time continuous, is commented on the discussion section.

4.2. A Discrete Time Version of the Leaky-Integrate and Fire Model

The time discretisation of the model (25) reads:

$$\left\{\begin{array}{ccc}{V}_{n+1}^k\hfill & =\gamma {\displaystyle V_n^k}+\sum _{j=1}^N{W}_{kj}{x}_n^j+S_n^k+{\sigma }_B{\xi }_n^k,\hfill & \mathrm{if}{\displaystyle V_n^k}<\theta ,\mathrm{Integrate}\mathrm{phase};\hfill \\ \hfill & \\ V_{n+1}^k\hfill & =0\mathrm{and}{x}_n^k=1,\hfill & \mathrm{if}{V}_n^k\ge \theta ,\mathrm{Firing}\mathrm{phase}.\hfill \end{array}\right.$$

(27)

For simplicity, we have assumed that all neurons have the same capacitance $C^k=C$, and set $\gamma =1-\frac{dt}{\tau }$, where ($\tau$ is the characteristic time scale of the membrane response, $dt$ is an integration time step which has to be much smaller than $\tau$ to preserve the physical relevance, whereas it has to be strictly positive to have a time-discretization scheme.) $\tau =RC$, with $0<\gamma <1$, and then taken $dt=1$. We have assumed that synapses are instantaneous. Subsequently, the synaptic input is ${\sum }_j{W}_{kj}{x}_n^j$, that correspond to the pre-synaptic neuron j that acts on the post-synaptic neuron k whenever j spikes, $x_n^j=1$. If, at some discrete-time n, $V_n^k$ exceeds the threshold $\theta$, the membrane potential is reset at time $n+1$ and a spike is recorded at n for neuron k, i.e., $x_n^k=1$. Below the threshold, the random dynamical system is ruled by (27). $S_n^k$ is the time discretization of the external stimulus. ${\xi }_n^k$ are independent standard Gaussian random variables.

It is easy to integrate Equation (27) conditionally upon a fixed spike sequence x. A trajectory $V=\left\{V_n^k,k=1\cdots N,n\in \mathbb{Z}\right\}$ is compatible with this spike sequence if $\chi \left(V_n^k>\theta \right)=x_n^k$, $\forall k=1\cdots N,n\in \mathbb{Z}$, where $\chi \left(A\right)$ is the indicator function of the logical event A, $\chi \left(A\right)=1$ if A is true, $\chi \left(A\right)=0$ otherwise. We discuss the compatibility condition in more detail in Section 4.4, when dealing with symbolic coding. For the moment, assume that V and x are compatible. We note ${\tau }_k(x,n)=max\left\{l,l<n|x_l^k=1\right\}$ the last time before n where neuron k has spiked, thus whose voltage was reset to 0. Then:

$$V_{n+1}^k=\sum _{j=1}^N{W}_{kj}{\eta }_{kj}(n,x)+\sum _{l={\tau }_k(x,n)}^n{\gamma }^{n-l}{S}_l^k+{\sigma }_B\sum _{l={\tau }_k(x,n)}^n{\gamma }^{n-l}{\xi }^k\left(l\right),$$

(28)

where:

$${\eta }_{kj}(n,x)=\sum _{l={\tau }_k(x,n)}^n{\gamma }^{n-l}{x}_l^j,$$

integrates the influence of pre-synaptic neuron j on the time interval $\left[{\tau }_k(x,n)+1,n\right]$. Each spike emitted by this neuron, at times l in this time interval, contributes with a weight ${\gamma }^{n-l}$ and there is no contribution at times where $x_l^j=0$. The condition $\gamma <1$ implies an exponential decay in the spike history dependence with characteristic time:

$${\tau }_{\gamma }=-\left[\frac{1}{log\gamma }\right].$$

Likewise, ${\sum }_{l={\tau }_k(x,n)}^n{\gamma }^{n-l}{S}_l^k$ integrates the stimulus influence on neuron k and ${\sigma }_B{\sum }_{l={\tau }_k(x,n)}^n{\gamma }^{n-l}{\xi }^k\left(l\right)$ is the integrated noise term. This is a Gaussian random variable, with mean zero and variance ${\sigma }_B^2\frac{1-{\gamma }^{2(n+1-{\tau }_k(x,n))}}{1-{\gamma }^2}$. In (28), the dependence on the initial condition does not appear because we assume the initial time to be $n_0\to -\infty$. So, either neuron k has spiked in the time interval $[-\infty ,n]$ and the voltage is reset to 0, or it has not spiked, but the initial condition dependence decays like ${\gamma }^{n-n_0}$ which vanishes when $n_0\to -\infty$.

4.3. Gibbs Distribution of the Discrete Lif Model

Thanks to the integrated Equation (28) and because the integrated noise is Gaussian it is now easy to compute the probability that neuron k spikes at time $n+1$ given the history x, $P\left(x_{n+1}^k=1|x_{-\infty ,n}\right)=P\left(V_{n+1}^k\ge \theta |x_{-\infty ,n}\right)$:

$$P\left(x_{n+1}^k=1|x_{-\infty ,n}\right)=f\left(\frac{\theta -{\sum }_{j=1}^N{W}_{kj}{\eta }_{kj}(n,x)-{\sum }_{l={\tau }_k(x,n)}^n{\gamma }^{n-l}{S}_l^k}{{\sigma }_B\sqrt{\frac{1-{\gamma }^{2(n+1-{\tau }_k(x,n))}}{1-{\gamma }^2}}}\right)$$

(29)

where $f\left(z\right)={\int }_z^{+\infty }\frac{e^{-\frac{u^2}{2}}}{\sqrt{2\pi }}dz$. Here, we have used a small abuse of notation. The conditioning upon $x_{-\infty ,n}$ means, in fact, the conditioning upon the sequence $x_{n-1,R\left(x\right)}$, where $R\left(x\right)=min\left\{k=1\cdots N|{\tau }_k(x,n)\right\}$. We condition upon the spike history prior to n back to the last time where all neurons had been reset at least once. While, for Equation (29). we just need to consider the history back to ${\tau }_k(x,n)$, the conditioning upon $x_{n-1,R\left(x\right)}$ is necessary when considering the conditional join probability of spiking patterns. The joint probability is conditionally independent given the past:

$$P\left(x_{n+1}|x_{-\infty ,n}\right)=\prod _{k=1}^N\left(x_{n+1}^{k}P\left(x_{n+1}^k=1|x_{-\infty ,n}\right)+\left(1-x_{n+1}^k\right)\left(1-P\left(x_{n+1}^k=1|x_{-\infty ,n}\right)\right)\right).$$

(30)

Let us comment this result. Equation (30) is the transition probability, of the form (20), where the normalised Gibbs potential $\varphi$ can be explicitly written, in terms of the synaptic interactions and the parameter ${\sigma }_B$ controlling the noise amplitude. However, note that, in contrast to (20), where the memory of the spike sequence was fixed independently of x, here the memory depends of x, providing a variable length Markov chain [115,116]. Actually, $R\left(x\right)$ is an unbounded function of x as one can find, for all $r\in \left\{-\infty ,n\right\}$, a sequence x, such that $R\left(x\right)=r$ (take the sequences where all $x_n^k=0$, $k=1\cdots N$, $n>r$ and $x_r^k=1$ for at least one k). Thus, we have to deal with the extension of Markov chains, to chains with unbounded memory introduced in Section 2.7. The existence and uniqueness of a Gibbs distribution compatible with this chain is guaranteed by the exponential decay of the memory controlled by $\gamma <1$ [23]. In this case, the potential fulfills the conditions described in Section 2.7. Finally, in (30), the transition probabilities explicitly depend on time because of the stimulus dependent term. They are, therefore, not translation invariant. While the extension of Gibbs distributions to non time-translation invariant chains can be rigorously done (upon the exponential decay of memory [90]), we restrict ourselves now to the case without stimulus ($S_n^k=0,\forall k=1\cdots N,n\in \mathbb{Z}$) to apply Thermodynamic Formalism, until Section 4.5, where we discuss linear response.

Note that (29) bares a strong resemblance to the GLM Ansatz (21).

4.4. Markov Partition and Symbolic Coding

In this section, we consider the deterministic discrete-time neuronal network model obtained considering (25) with ${\sigma }_B=0$, studied in detail in [114]. The threshold $\theta$ of the voltage in a network of N neurons induces a natural partition of ${\mathbb{R}}^N$, $\mathcal{P}={\prod }_{k=1}^N{\mathcal{P}}_{x^k}$, where $x^k \in \left\{0,1\right\}$, ${\mathcal{P}}_0=[-B,\theta ],{\mathcal{P}}_1=[\theta ,B]$ where the bound B depends on synaptic weights [114,117]. If $V_n^k\in {\mathcal{P}}_0$, it evolves according to the sub-threshold dynamics (25) and it does not emit a spike at time n. In contrast, if $V_n^k\in {\mathcal{P}}_1$, it emits a spike at time n and its trajectory is set back to ${\mathcal{P}}_0$ at time $n+1$. Thus, $\mathcal{P}$ is a natural partition in the sense that it informs about the spikes of each neuron. Therefore, to each trajectory $V\equiv \left\{V_n^k,k=1\cdots N,n\in \mathbb{Z}\right\}$, there is associated an infinite spike sequence x such that $x_n^k=0⟺V_n^k\in {\mathcal{P}}_0$ and $x_n^k=1⟺V_n^k\in {\mathcal{P}}_0$.

This partition can be used to generate a Markov partition [114], but, in general, the Markov partition is not $\mathcal{P}$ but a finite refinement $\mathcal{Q}$ of $\mathcal{P}$. What ensures that a Markov partition exists is that (27) is contracting. More precisely, it contracts, in one step, at speed $\gamma$ for directions (neurons) such that $V^k<\theta$, and it contracts, with an infinite speed, for directions such that $V^k\ge \theta$ (reset). However, this generates discontinuities in the mapping (27), and a singularity set $\mathcal{S}=\left\{V\in {\mathbb{R}}^N|\exists k\in \left\{1\cdots N\right\},V^k=\theta \right\}$, where the map associated with (27), hereafter denoted by $\mathcal{G}$, is discontinuous. Thus, $\mathcal{G}$ is piecewise continuous and piecewise contracting.

Now, recall that $\mathcal{Q}$ is a Markov partition for the dynamics with contracting map $\mathcal{G}$ if its elements satisfy $\mathcal{G}\left({\mathcal{Q}}_n\right)\cap {\mathcal{Q}}_{n^{\prime }}\ne \varnothing \Rightarrow \mathcal{G}\left({\mathcal{Q}}_n\right)\subset {\mathcal{Q}}_{n^{\prime }}$. In other words, the image of ${\mathcal{Q}}_n$ is included in ${\mathcal{Q}}_{n^{\prime }}$ whenever the transition $n\to n^{\prime }$ is legal. Here, in general, the elements of $\mathcal{P}$ do not satisfy this condition. This is because the image of a domain of $\mathcal{P}$ usually intersects in several domains (in this case, the image intersects the singularity set). From the neural network’s perspective this means that, in general, it is not possible to know the spiking pattern at time $n+1$ knowing the spiking pattern at time n. There are several possibilities, depending on the membrane potential values and not only on the firing state of the neurons. Of course, if, say ${\mathcal{P}}_n$ is such that $\mathcal{G}\left({\mathcal{P}}_n\right)$ intersects several domains ${\mathcal{P}}_{n_1},\cdots ,{\mathcal{P}}_{n_l}$ one can take the preimages of these domains ${\mathcal{G}}^{-1}\left({\mathcal{P}}_{n_r}\right)$ to construct a refinement of $\mathcal{P}$, such that the Markov partition requirement is satisfied in one iteration of the map. However, nothing guarantees that, at the second iteration, some elements of this new partition will not intersect the singularity set under ${\mathcal{G}}^2$.

Can we find a finite refinement of $\mathcal{P}$, such that the trajectory of the partition elements never intersects several partition elements? It is shown in [114] that this property is satisfied for generic values of the synaptic weights $W_{ij}$. Essentially, it is based on the fact that the distance between the $\Omega$-limit set of (27) and the singularity set $\mathcal{S}$, is generically positive. In other words, each point in the partition elements ${\mathcal{Q}}_n$ has a local stable manifold with a finite diameter.

As a consequence, the deterministic discrete-time neuronal network model (27) admits a Markov partition, a refinement of the natural partition, providing a symbolic coding of the membrane potential trajectories in terms of spike sequences. In particular, once the initial condition dependence has been removed, the evolution (28) (without noise) is only constrained by the stimulus. Thus, (28) provides a coding scheme of the stimulus in terms of spike sequences. The Markov partition is made of spike blocks, with finite memory depth R, which can be used to apply Thermodynamic Formalism in the presence of noise. However, R—the memory depth of the corresponding Markov chain—depends on the parameters and, in particular, the synaptic weights.

In addition, note that the presence of a singularity set induces a weak form of initial conditions. Although the dynamic is contracting, a small perturbation of a trajectory can induce an evolution drastically different from the unperturbed trajectory, if the perturbation crosses the singularity set. In this case, e.g., there is a neuron, k, which does not spike, at time n in the unperturbed trajectory, and spikes at time n in the perturbed one, inducing a completely different evolution (cascade effect). This phenomenon has been exposed in the context of spiking neurons, where the coexistence of stable and unstable dynamics is investigated [118]. The singularity set also induces the existence of ghost orbits, $\exists k\in \left\{1\cdots N\right\},\forall n>0,V_n^k<\theta$ and ${lim\; sup}_{n\to +\infty }{V}_n^k=\theta$. However, ghost orbits are non-generic in a topological and a measure-theoretic sense. As a corollary, the $\Omega$-set is generically composed of finitely many periodic orbits with a finite period (whose length depends on parameters of the model, in particular, synaptic weights).

4.5. Extensions

4.5.1. Explicit Form of the Potential. GLM vs. MaxEnt

Model (27) makes a link between the dynamics of a neuronal network and the transition probabilities (20), where the dependence on the model parameters (in particular, synaptic weights, and stimuli) is explicit. We have an explicit potential for this model, which, here, takes a GLM-like form (29), but is more general, as in contrast to the GLM, the effective interactions depend on time via powers of the leak term $\gamma$. This potential can also be written in terms of monomials using the Hammersley–Clifford decomposition (Section 3.3), through a series expansion of the function f. This procedure generates a series of monomials with coefficients that can be explicitly computed (using the fact that, from the monomials definition (22) ${\left(x_{i_k}^{n_k}\right)}^m=x_{i_k}^{n_k}$, for any integer $m>0$). These coefficients are proportional to powers of $\gamma <1$, so their strength decays exponentially fast, allowing for truncating the potential to a finite number of terms, which produce canonical Markov approximations of different orders [92]. One obtains, to the lowest order, a Bernoulli potential, then pairwise terms, and so on.

4.5.2. Linear Response

Another interesting consequence of the analysis of this model is that the potential may depend on a time dependent stimulus. When considering that the stimulus is of small amplitude and additive, one can take a Taylor expansion of the potential as powers of the stimulus allowing one to go beyond the stationarity assumption central to equilibrium statistical mechanics and Thermodynamic Formalism. In this case, it is possible to show that the variations in the spike statistics induced by the stimulus, can be described in terms of a linear response theory [119,120,121,122]. The main result is that the variation, in the average of an observable f, resulting from the application of a stimulus reads:

$$\delta \mu \left[f\right]\left(n\right)\equiv {\mu }^S\left[f\right]\left(n\right)-{\mu }^{\left(sp\right)}\left[f\right]=\left(K_f\ast S\right)\left(n\right),$$

where ${\mu }^{\left(sp\right)}$ is the Gibbs distribution in spontaneous activity (without stimulus), and ${\mu }^S$ is the Gibbs distribution in the evoked activity regime (with stimulus), ${\mu }^S\left[f\right]\left(n\right)$ means the average of f with respect to ${\mu }^S$, which depends on time (if the stimulus does), and ${\mu }^{\left(sp\right)}\left[f\right]$ means the average of f with respect to ${\mu }^{\left(sp\right)}$, which does not depends on time. This variation in average is given by a convolution between a kernel $K_f$, depending on f (which can be expressed in terms of time correlation functions between monomials) and of the stimulus. The coefficients in the expansion of $K_f$ depend on the parameters constraining dynamics (e.g., the synaptic weights in (27)). The correlations are computed with respect to the invariant Gibbs measure ${\mu }^{\left(sp\right)}$. In addition, the influence of monomials in the expansion decreases with their order, so that one can obtain a reasonable approximation of the convolution kernel only considering averages of order two monomials (time dependent pairwise correlations). Therefore, this is sa result in the form of a fluctuation-dissipation “theorem” in statistical physics, with the difference that one considers time dependent correlations. This formula has proven to give astonishingly good results when computing the response to a time dependent stimulus in the model (27) [122].

4.6. The Galves-Löcherbach Model

Here, we present a second example, where the Gibbs potential can be computed. This model is known as the Galves–Löcherbach model introduced by Antonio Galves and Eva Löcherbach in [91] (see also [123,124]). This model is a generalization of [22], but considering an infinite (countable) network of neurons interacting in time with memory of variable length.

The model is built when considering a stochastic chain ${\left(X_t\right)}_{t\in \mathbb{Z}}$ taking values in ${\{0,1\}}^I$, where I is a countable set of neurons. The probability of a spike depends on the accumulated activity of the system since the last spike, thus, each spike depends on a variable length history, defining also a non-Markovian stochastic process. Extensions of this model have been made considering the hydrodynamic limit of the interacting neuronal system [125], classifying the collective behavior according to parameter values [126], and the generalization to the continuous time [127,128].

For each neuron $i\in I$ and each time $t\in \mathbb{Z}$, let ${\tau }_i(x,t)$ denote the last time before t at which neuron i fired a spike in the spike train x:

$${\tau }_i(x,t)=sup\left\{s<t:x_s^i=1\right\},$$

and suppose that the synaptic weights $W_{ij}$ have the uniform summability property:

$$\underset{i\in I}{sup}\sum _j\left|W_{ij}\right| < \infty .$$

The joint probabilities are conditionally independent given the past:

$$P\left(x_t\mid x_{-\infty ,t-1}\right)=\prod _{i\in I}P\left(x_t^i\mid x_{{\tau }_i(x,t),t-1}\right),$$

(31)

where the probability of neuron i having a spike at time t is given by:

$$P\left(x_t^i\mid x_{{\tau }_i(x,t),t-1}\right)=h_i\left(\sum _j{W}_{ij}\sum _{s={\tau }_i(x,t)}^{t-1}{g}_j(t-s)x_s^j\right).$$

(32)

where $g_j(t-s)$ plays the role of the exponential $\alpha$-kernel in (26). The transition probabilities (32) have the form of a GLM model.

Under technical conditions of the functions $h_i$ and $g_j$ and $W_{ij}$, there exists a unique probability measure consistent with (31) and (32) (see Theorem 1 of [91]). To prove this claim, they use a Kalikow-type decomposition of the infinite order transition probabilities. This type of decomposition has also been considered in Ref [91,129,130]. The setup considered in this work extends to infinite size and infinite memory.

5. Discussion and Perspectives

In this review, we introduce different ideas and tools from Thermodynamical Formalism and show how they can be applied in theoretical neuroscience. As a summary, we grouped these approaches depending on two main characteristics. The first one is the number of neurons N that affect the cardinality of the alphabet considered and the range R of the potential (memory in transition probabilities) associated to the equilibrium measure characterising the system. This is represented in

Table 1. Types of Gibbs measures potentially found in experimental data analysis or in the analysis of mathematical models of networks of interacting spiking neurons.

Thermodynamic Formalism and Gibbs Measures
Number of neurons	Memory of the potential
	Memoryless	Finite	Infinite
Finite	Boltzmann-Gibbs	Gibbs in the sense of Bowen	Chains with complete connections
Infinite	Countable state Bernoulli	Countable state Markov	Chains with variable length

. The infinite number of neurons and infinite range cases are further discussed in Section 5.2.

While there have been interesting applications of Thermodynamic Formalism to neuroscience, there are interesting ideas and developments still to come. In this concluding section we raise questions, challenges and new avenues for the application of Thermodynamic Formalism to theoretical neurosciences.

5.1. Thermodynamic Formalism for More Biologically Plausible Neuron Models

In this review, we have considered rather academic models of neuronal networks, where, especially, time is considered to be discrete. There are good reasons for that. As we remarked at the beginning of Section 3, we were considering models, like the Integrate-and-Fire, where spikes arise instantaneously, in continuous time, thereby providing a possibly uncountable set of potential spike trains. The question is whether we are dealing here with a realistic property of biological neuronal networks or with an artifact created by the instantaneous reset. Real spikes have a duration (a few ms) and a refractory period (also a few ms), so, for a fixed initial condition, spike trains produced by a continuous-time neuronal network are countable. Now, it might be that the set of spike trains depend continuously on the initial condition, so that we are still left with an uncountable set of spikes.

It is out of the scope of this review to discuss from a biological perspective, whether or not neuronal networks have the cardinality of the continuum (see [101,117,131] for a discussion on this topic). Instead, we have considered a strategy consisting of discretizing time, avoiding the problem of potentially uncountable spikes. Here, we briefly mention another strategy, allowing for associating a countable set of spike trains to continuous time networks with a countable set of spike trains. We first make the remark that the instantaneous reset of voltage is physical and biological nonsense, inducing pathologies in the dynamics [132]. On this basis, we use a convenient mathematical trick that is explained in the next paragraph, which can certainly be criticized on phenomenological grounds [131], especially when the dynamical system representing the neuronal activity is deterministic.

After spiking, a biological neuron stays at rest a certain time (refractory period). Accordingly, the trick is the following. Fix $\delta >0$ and define a spiking variable $x_n^k\in \left\{0,1\right\}$, where n is an integer, where $x_n^k=1$ if neuron k emits a spike in the time interval $[n\delta ,(n+1\left)\delta \right]$ and $x_n^k=0$ otherwise. Recall that $t_{\left(r\right)}^k$ denotes the time at which neuron k emits its r-th spike. This reads:

$$x_n^k=\left\{\begin{array}{cc}1,\hfill & if\exists r,t_{\left(r\right)}^k\in [n\delta ,(n+1)\delta ];\hfill \\ 0,\hfill & \mathrm{otherwise}.\hfill \end{array}\right.$$

Spiking variables are therefore time-discrete events with a time resolution $\delta$. When $V^k$ reaches the threshold at time $t_{\left(r\right)}^k$ it is reset to 0, and stays there until time $(n+1)\delta$. After this, follows the sub-threshold evolution (25) until the next time where $V^k$ reaches the threshold. Note that, in this modelling, $\delta$ can be quite small when compared to the time scales of the dynamics. In this way, the set of spike trains x becomes at most countable.

This trick can be used to generalise the Integrate-and-Fire model into a conductance-based Integrate-and-Fire model, which was introduced by Rudolph and Destexhe in [133] and mathematically studied in [23,24,117], where the synaptic conductance depends on the spike history. One can still show that a unique Gibbs distribution with infinite range potential exist in this case, characterising the spike train statistics. The potential can be explicitly computed as a function of network parameters, even in the presence of a time-dependent stimulus.

Now, would Thermodynamic Formalism apply to more realistic neuronal models, like Hodgkin–Huxley [111], FitzHugh–Nagumo [134,135], Morris–Lecar models [136]? (see [36,109,137] for a complete presentation of canonical neuronal models). In these models, closer to biology, the spikes have a time course and, thus, are not considered as point events. However, we do not know any result establishing, e.g., the existence of Markov partitions and symbolic coding for these models and this seems to be out of reach for the moment. Still, one can bin the time and proceed as done in experiments where voltage is a time-continuous signal. Thus, one can still use the approach used in (20) to characterise the spike train statistics.

A natural question in this context is what is the link between spike train statistics and the underlying dynamical model, with “hidden” dynamical variables, such as membrane potential, but also, e.g., activation/inactivation variables? If we think in terms of spike coding, the alternative is the following. Either spike trains contain all the necessary information to characterise the dynamics, e.g., the spike response to a stimulus, and then characterising (20) is sufficient. Or, there is additional information, not conveyed by spikes (e.g., sub-threshold oscillations [100,138,139]), and the "neural code" is not entirely contained in the spikes, somewhat ruining the hope of encoding neuronal messages purely in terms of spikes. This question is much more general than the validity of the Thermodynamic Formalism approach for these models.

What Thermodynamic Formalism brings to the analysis of these models is twofold: (1) a way to rigorously handle probabilistic representations of spikes (20); and, (2) to provide conceptual and mathematical tools to analyse spiking neuronal network models, like (27), where a dynamical system formulation of biophysical variables can be mathematically related to spike coding and spike train statistics.

5.2. Phase Transitions

Several studies have shown that the population of vertebrate retinal ganglion cells responding to naturalistic stimulus is poised near a “critical state” [73,74]. From the maximum entropy joint distribution (13), a family of Gibbs distributions can be built introducing a parameter $1/\beta$ (analogous to the inverse temperature), which scales all of the Lagrange multipliers of the inferred Hamiltonian. When $\beta \to 0$ (infinite temperature), the uniform distribution is obtained, and when $\beta \to +\infty$ (zero temperature), the Dirac delta supported at the spike configuration(s) of minimal energy is obtained. $1/\beta =1$ corresponds to the inferred maximum entropy distribution from data. These studies have only analysed memoryless Gibbs distributions (13).

From this representation, it is possible to compute the fluctuations (variance) of the energy $U_{\lambda }$ as a function of the “temperature” parameter T. This quantity can be obtained as the second derivative of the pressure, Equation (6), which is, in thermodynamics, related to the heat capacity $C_T$. On numerical grounds, this quantity can be computed while using MonteCarlo simulations and plotted as a function of the “temperature” $T=\frac{1}{\beta }$, for different network sizes, (see Figure 4).

The form of $C_T$ versus T plot, for maximum entropy models of Ising type obtained from the recording of retinal ganglion cells responding to naturalistic stimuli are shown in Figure 4 (redrawn from [74]). It can be observed that there is a clear, increasing peak at $T=1$, which starts to manifest itself when larger and larger groups of neurons are considered. This presumable divergence of the heat capacity (or variance of U) when $N\to \infty$ (thermodynamic limit) is interpreted as a second order phase transition (a so-called “critical regime” [140]).

The behavior of the specific heat that is observed in Figure 4 suggests that the heat capacity of a maximum entropy distribution, fitted over an increasingly large group of neurons in the retina, diverges. This phenomenon has been considered to be a “signature of criticality" (details of this study and a discussion about whether criticality is functional for retinal ganglion cells can be found in Ref. [74]). Some criticism regarding this approach to diagnostic criticality has appeared arguing that the maximum entropy principle is likely to yield models that are close to singular values of parameters, akin to critical points in physics where phase transitions occur. Statistically distinguishable models tend to accumulate close to critical points, where the susceptibility (inverse Fisher Information) diverges in infinite systems [141]. These ideas have also been applied to numerical simulations of a canonical feed-forward population model showing that the specific heat diverges whenever the average correlation strength is independent of the population size [75], as in the random subsampling of correlations used in [74]. Additionally, note that, for spike trains obtained from discrete Markov processes, binning generates a stochastic process with unbounded memory akin to inducing spurious phase transitions [102].

This interesting approach leads, nevertheless, to several questions in the context of Thermodynamic Formalism.

• Does this signature of criticality extend to Gibbs distributions with potentials of range $R>1$, i.e., with memory? How does it depend on R? We are not aware of any experimental results addressing this issue. This question is related to the following:
• What is this signature of criticality from the point of view of Thermodynamic Formalism? The occurrence of a second-order phase transition mathematically means that the pressure is $C^1$ but not $C^2$ when some limit is taken. Here, we have two possible limits: the range of the potential R tends to infinity or the number of neurons N tends to infinity. These two limits could also be addressed simultaneously and they do not necessarily commute. For potentials associated to finite R and N the Perron-Frobenius theorem guarantees the existence and uniqueness of the Gibbs measure and the analyticity of the pressure can also be proved, preventing phase transitions. When R or N are infinite, the properties of the RPF operator (Section 2.3.2) characterises the presence or absence of phase transitions. Indeed, there are conditions ensuring a spectral gap for this operator, ensuring the exponential decay of correlations. Now, Equation (7) characterise the second derivative of the pressure as a time series of correlations, which converge when the correlations decay exponentially. On the opposite side, the non-summability of time correlation function implies the non-existence of the second derivative, and thus, of a second-order phase transition. Therefore, a possibility to have a second-order phase transition is when the spectral gap property for the RPF operator when $R\to +\infty$ or $N\to +\infty$ is absent. In statistical mechanics, second-order phase transitions can be characterised by how the zeros of the partition function, written as a polynomial, pinch the real axis (Lee-Yang phenomenon) [142,143,144]. In our case, when $R>1$, the object of interest is not the partition function, but rather the largest eigenvalue of ${\mathcal{L}}_{\varphi }$, which has to stay analytic in the limit R, or N, $\to +\infty$. The absence of the spectral gap property presents an analogy with the Lee-Yang phenomenon, although we do not know about results establishing a deeper link.
• Can we relate known examples of dynamical systems exhibiting phase transitions to models in neuroscience? Another possible example to be interpreted in neuroscience is the Dyson model [145], in which there exists a phase transition in the sense of spontaneous magnetisation when the temperature goes to zero, due to an infinite range potential whose correlation does not decay exponentially fast. In our case, the range of the potential should be taken in time, keeping (possibly) the number of neurons finite. Other examples exist of rigorous characterisations of phase transitions in the thermodynamic description of Pomeau–Manneville intermittent maps, passing from an integrable density function associated with the measure to heavy-tailed densities [146]. An interesting result may hint at the connection between the topological Markov map of the interval and stochastic chains of infinite order or chains with complete connections. Ref [147] presents how to build a topological Markov map of the interval whose invariant probability measure is the stationary law of a given stochastic chain of infinite order. This is interesting in this context because as we presented in (27), there are mathematical models of spiking neurons whose spike statistics are represented by chains of infinite order. This result or its inverse i.e, how to build a stochastic chain of infinite order from a topological Markov map may hint at conditions in the parameters or conditions of the mathematical models of spiking neurons to exhibit second order phase transitions.
• What could the dynamical or mechanistic origins of a second-order phase transition be in a spiking neuronal network model? Handling experimental data is of course important, but for long experiments with living neuronal tissue, one cannot control the size of the sample, the stationarity of data, and so on. Accordingly, assume that we have been able to find an example of a Gibbs distribution exhibiting a second-order phase transition when $R\to +\infty$ or $N\to +\infty$. Can we build a spiking dynamical system, with finite R and N, which has this Gibbs distribution in the limit R, or N, $\to +\infty$, so that we observe a phase transition in the model? Then, what are the mechanistic origins (in the neuronal dynamics) of second-order phase transitions? It could be an interesting example to study the existence of a second-order phase transition in a simple neuronal model. Returning back to the discrete LIF model, the failure in the second-order differentiability of the pressure means the loss of exponential mixing, which, in the model (27) can arise in, at least, two cases. First, if $\gamma =1-ϵ$, $ϵ\to 0$. This is a way to obtain a potential with increasing range as $ϵ\to 0$ with loss the summability of correlations. The corresponding orbits (reminiscent of the ghost orbits discussed in Section 4.4) are such that it may take a long time for some neurons to be reset. Thereby, the memory to be considered is very long. However, this is a case hardly interpretable from the neuroscience perspective. A second possibility is to analyse how the pressure depends on the spectrum of the synaptic weights matrix and to check whether there are cases (e.g., small world or scale-free lattices), where the spectral gap of the RPF vanishes.

From the perspective of the maximum entropy distributions built from experimental data of spiking neurons, there have been interesting applications of the Gibbs distributions that were obtained to answer questions related to the retinal code that are not related to criticality [148]. From the maximum entropy joint distribution the conditional distributions can be computed, and questions about the redundancy of the neural code can be addressed such as how predictable is the activity of each neuron based on the knowledge of the activity of other neurons in the population. Can we find a subset of neurons J that together predict with high accuracy the spiking behaviour of the neuron i? Mathematically can be written in this way $p(x^i=1\mid {\left\{x^j\right\}}_{j\in J})$, where J is a subset of neurons in the population of spiking neurons. Other questions related to the neural coding and dimensionality reduction can be addressed studying the energy landscape $U_{\lambda }\left(x\right)$ of (13). For example, the the local minima of an energy landscape correspond to metastable states and several configurations may correspond to the same “valley“ near each local minima. Transitions between valleys have be studied in the context of “retinal coding” (see details in [148]). Alternative methodologies using the maximum entropy principle to study network of sensory neurons have been used in order to classify intrinsic interactions from extrinsic correlations [46] and to reveal the excitatory and inhibitory correlations [45].

5.3. What Else Do Thermodynamic Formalism and Gibbs Distributions May Tell Us about Neuroscience?

The relationship between mathematics and physics has been historically symbiotic and Thermodynamic Formalism is an interesting example of how ideas from physics may help to solve problems and introduce ideas into the field of mathematics. The history of Thermodynamic Formalism also shows how purely mathematical results can be obtained as corollaries of physical laws, inverting the frequently assumed relationship between physics and mathematics [149].

However, in the case considered in this review—the link between Thermodynamic Formalism and neuroscience—the mathematical problem is motivated by biology, not by physics. While Eugene Wigner argues in favour of the “The unreasonable effectiveness of mathematics in the natural sciences” [150], Israel Gelfand, after spending several years working in mathematical problems related to biology, replied with “The unreasonable ineffectiveness of mathematics in biology.” [151]. While there are reasons to argue that this is still the case, it is less clear if one can blame the field of mathematics or just the fact that we have not yet used the right tools or frameworks. In the quest for these “right tools”, there is a long tradition of using ideas from statistical physics to study neural networks, and in particular, to represent the emergence of collective behaviour from microscopic interactions, with the hope that statistical aspects of the collective behaviour will be independent of the details in these systems. This gave rise to major branches of theoretical neuroscience, like dynamic mean-field methods [29,152,153,154] or the Maximum Entropy approach [38,49,148,155], mainly coming from physicists. Although there are considerably less articles using mathematical methods to rigorously analyse the collective behaviour of neuronal networks some promising approach have been recently proposed based on large deviations [156,157], Kalikow-type decomposition [91,158], stochastic processes [159,160,161,162,163], dynamical systems [137,164], etc. As we have developed in this review, Thermodynamic Formalism could also be one of these tools, providing interesting connections between mathematics and physics, dynamics and statistics, applied to neuroscience.

Especially, we have described how Thermodynamic Formalism: (1) provides a conceptual and operational (i.e., allowing to develop algorithms and software [99]) framework to analyse experimental spike train data; (2) allows to derive explicit expressions linking spike statistics to neuronal networks dynamics; (3) extends to non stationarity via linear response theory; and, (4) proposes a realm to address questions related to criticality.

Here, we would like to propose some other extensions, not yet explored so staying at the level of ideas, all based on the power of Thermodynamic Formalism to make explicit and operational links between dynamics, statistics and symbolic coding.

• Geometry of the state space. A prominent aspect of Thermodynamic Formalism, that we haven’t discussed yet in this review, is its link to the characterisation of the geometry of attractors and, especially, fractal sets [165,166]. For example, the composition of contracting mappings along symbolic orbits defines the so-called Iterated Function Systems (IFS) [167] generating fractal sets with tunable geometry and structure. Now, it is interesting to remark that Integrate and Fire models are actually piecewise contracting dynamical systems having a structure similar to IFS where the contracting pieces are symbolically encoded by spike blocks [114]. It would be interesting to investigate, along these lines, the structure of attractors in Integrate and Fire models, and how orbits, encoded by spike blocks, are related to the geometry of attractors (the $\Omega$-limit set).
• Transitions between attractors. The concept of attractor is actually central in describing brain dynamics [168,169]. Especially, a current trend in neuroscience is to associate to brain states attractors (or ghost attractors, see [170] and references therein). The transitions between these states corresponds to transition during tasks or spontaneous activity [171,172,173,174]. It is relatively natural to characterise such transitions by Markov chains [175], which is the first step toward the application of Thermodynamic Formalism and analysing these transitions from a statistical and statistical physics perspective.
• Non-stationarity and link with generating functional formalism. As we mentioned, Thermodynamic Formalism is constructed from a variational approach based on entropy and, thus, requiring time translation invariance. We have briefly described how we can depart from this constraint while using linear response theory. It would be interesting to explore beyond this point and consider general types of response to stimuli (not requiring a small perturbation, as in linear response). For this, one would have to construct a Thermodynamic Formalism based on the optimisation of a quantity, which is not the entropy. This is somehow what generating functional approaches like the dynamic mean-field theory does (see introduction), although using other constraining hypotheses (essentially to be able to describe the infinite size limit by a Gaussian process). It would be interesting to try to close the gap between these two approaches (e.g., via large deviations theory).

One of the biggest challenges in science of the XXI century is to understand the brain functions within a conceptual framework that are capable of unifying the multi-scale dynamics that take place in the brain. This framework should also make sense in the light of the overwhelming amount of experimental data capable of predicting macroscopic phenomena, such as motor behaviour or visual experience from the activity of billions of neurons.

Physicists have been able to make a deep connection between mechanics, statistical physics, and thermodynamics. A similar quest is presumably guiding the research of (some) theoretical and experimental neuroscientists. While there is still a long way to go before achieving this goal (as some argue we are still searching for principles [176]), during the last decades, mathematicians have been playing a relevant role in the rigorous description of neural phenomena, clarifying and raising conceptual problems in neuroscience.

We hope that theoretical tools and ideas from Thermodynamic Formalism and its current application in neuronal dynamics and spike train statistics will lead to a better and unified understanding of the neural phenomena. We also hope that the present review may serve as an encouragement for the mathematical community that is interested in applications of Thermodynamic Formalism in order to study these interesting and important problems.