# Ising-like model replicating time-averaged spiking behaviour of in vitro neuronal networks

Cesar I N Sampaio Filho, Lucilla de Arcangelis, Hans J Herrmann, Dietmar Plenz, Patrick Kells, Tiago Lins Ribeiro, José S Andrade

► To cite this version:

Cesar I N Sampaio Filho, Lucilla de Arcangelis, Hans J Herrmann, Dietmar Plenz, Patrick Kells, et al..
Ising-like model replicating time-averaged spiking behaviour of in vitro neuronal networks. Scientific
Reports, 2024, 14 (1), pp.7002. 10.1038/s41598-024-55922-9 . hal-04582213

HAL Id: hal-04582213

https://u-paris.hal.science/hal-04582213

Submitted on 21 May 2024

**HAL** is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Ising-like model replicating time-averaged spiking behaviour of in vitro neuronal networks

Cesar I N Sampaio Filho, Lucilla de Arcangelis, Hans J Herrmann, Dietmar Plenz, Patrick Kells, Tiago Lins Ribeiro, José S Andrade

► To cite this version:

Cesar I N Sampaio Filho, Lucilla de Arcangelis, Hans J Herrmann, Dietmar Plenz, Patrick Kells, et al..
Ising-like model replicating time-averaged spiking behaviour of in vitro neuronal networks. Scientific
Reports, 2024, 14 (1), pp.7002. 10.1038/s41598-024-55922-9 . hal-04582213

HAL Id: hal-04582213

https://u-paris.hal.science/hal-04582213v1

Submitted on 21 May 2024

**HAL** is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

Check for updates

OPEN

# Ising-like model replicating time-averaged spiking behaviour of in vitro neuronal networks

Cesar I. N. Sampaio Filho¹✉, Lucilla de Arcangelis², Hans J. Herrmann¹,³, Dietmar Plenz⁴, Patrick Kells⁴, Tiago Lins Ribeiro⁴ & José S. Andrade Jr.¹

We analyze time-averaged experimental data from in vitro activities of neuronal networks. Through a Pairwise Maximum-Entropy method, we identify through an inverse binary Ising-like model the local fields and interaction couplings which best reproduce the average activities of each neuron as well as the statistical correlations between the activities of each pair of neurons in the system. The specific information about the type of neurons is mainly stored in the local fields, while a symmetric distribution of interaction constants seems generic. Our findings demonstrate that, despite not being directly incorporated into the inference approach, the experimentally observed correlations among groups of three neurons are accurately captured by the derived Ising-like model. Within the context of the thermodynamic analogy inherent to the Ising-like models developed in this study, our findings additionally indicate that these models demonstrate characteristics of second-order phase transitions between ferromagnetic and paramagnetic states at temperatures above, but close to, unity. Considering that the operating temperature utilized in the Maximum-Entropy method is $T_o = 1$, this observation further expands the thermodynamic conceptual parallelism postulated in this work for the manifestation of criticality in neuronal network behavior.

In 1985 Amit et al.¹ for the first time drew a connection between neural networks and Ising spin glasses. Both systems have in common an energy landscape of many valleys and offer the possibility of delocalized storage of patterns. These analogies were explored further since then²,³. Subsequently, 1991 Miranda and Herrmann⁴ suggested that the brain operates at criticality in the sense that it exhibits avalanches of activity following power-law distributions for their size and duration, with exponents around – 1.5 and – 2.0, respectively. Alternative proposals were made later by Chialvo and Bak⁵. These theoretical predictions were experimentally confirmed by the seminal work of Beggs and Plenz⁶. By recording spontaneous local field potentials using a 60 channel multielectrode array on mature organotypic cultures of acute slices of rat cortex, they found power-law distributions in avalanche size and duration with similar exponents values as reported in⁴,⁵. Since then, many authors have confirmed signs of criticality in the brain⁷⁻⁹ overcoming the substantial challenges of studying scaling laws reported at criticality in the face of experimental constraints when assessing brain dynamics¹⁰,¹¹.

Models from Statistical Physics, when used in conjunction with the Maximum Entropy Method (MEM) developed in information theory, provide a conceptual framework to understand, from experimental data, a given natural process in terms of the “interactions” among its many elementary units¹²⁻¹⁴. The principle of maximum entropy, by itself, contains the essence of the so-called Inverse Ising Problem solution, in which a “Hamiltonian” associated with a given complex system can be inferred from observed statistical correlations among its components. This method is frequently referred to as the Boltzmann machine, since the Boltzmann distribution is present at its core. Generally speaking, the MEM has been applied to systems that can be mapped to Ising-like models, namely, models in which the interacting elements are in an active or inactive state, i.e., a lattice of dipole moments with spins states that are either up or down under the action of an external field and their mutual interactions. For instance, when dealing with neuronal networks, the interactions between neurons reacting to some stimuli are inferred from their firing patterns¹⁵⁻²¹, in which the firing of a spike equate to the spin up state. This strategy has also been successful in characterizing protein-protein interactions²²,²³, the genetic interaction networks from gene expression patterns²⁴⁻²⁶, the collective responses exhibited by flocks of birds²⁷,²⁸,

¹Departamento de Física, Universidade Federal do Ceará, Fortaleza 60451-970, Brazil. ²Department of Mathematics and Physics, University of Campania “Luigi Vanvitelli”, 81100 Caserta, Italy. ³PMMH, ESPCI, CNRS UMR 7636, 7 Quai St. Bernard, 75005 Paris, France. ⁴Section on Critical Brain Dynamics, NIMH, Bethesda, MD 20892, USA. ⁵email: cesar@fisica.ufc.br
<figure><img src="image_2.png" /></figure>

Check for updates

OPEN

# Ising-like model replicating time-averaged spiking behaviour of in vitro neuronal networks

Cesar I. N. Sampaio Filho<fnref n="1" />✉, Lucilla de Arcangelis², Hans J. Herrmann<fnref n="1" />,³, Dietmar Plenz⁴, Patrick Kells⁴, Tiago Lins Ribeiro⁴ & José S. Andrade Jr.<fnref n="1" />

We analyze time-averaged experimental data from in vitro activities of neuronal networks. Through a Pairwise Maximum-Entropy method, we identify through an inverse binary Ising-like model the local fields and interaction couplings which best reproduce the average activities of each neuron as well as the statistical correlations between the activities of each pair of neurons in the system. The specific information about the type of neurons is mainly stored in the local fields, while a symmetric distribution of interaction constants seems generic. Our findings demonstrate that, despite not being directly incorporated into the inference approach, the experimentally observed correlations among groups of three neurons are accurately captured by the derived Ising-like model. Within the context of the thermodynamic analogy inherent to the Ising-like models developed in this study, our findings additionally indicate that these models demonstrate characteristics of second-order phase transitions between ferromagnetic and paramagnetic states at temperatures above, but close to, unity. Considering that the operating temperature utilized in the Maximum-Entropy method is $T_o = 1$, this observation further expands the thermodynamic conceptual parallelism postulated in this work for the manifestation of criticality in neuronal network behavior.

In 1985 Amit et al.<fnref n="1" /> for the first time drew a connection between neural networks and Ising spin glasses. Both systems have in common an energy landscape of many valleys and offer the possibility of delocalized storage of patterns. These analogies were explored further since then²,³. Subsequently, 1991 Miranda and Herrmann⁴ suggested that the brain operates at criticality in the sense that it exhibits avalanches of activity following power-law distributions for their size and duration, with exponents around – 1.5 and – 2.0, respectively. Alternative proposals were made later by Chialvo and Bak⁵. These theoretical predictions were experimentally confirmed by the seminal work of Beggs and Plenz⁶. By recording spontaneous local field potentials using a 60 channel multielectrode array on mature organotypic cultures of acute slices of rat cortex, they found power-law distributions in avalanche size and duration with similar exponents values as reported in⁴,⁵. Since then, many authors have confirmed signs of criticality in the brain⁷⁻⁹ overcoming the substantial challenges of studying scaling laws reported at criticality in the face of experimental constraints when assessing brain dynamics¹⁰,¹¹.

Models from Statistical Physics, when used in conjunction with the Maximum Entropy Method (MEM) developed in information theory, provide a conceptual framework to understand, from experimental data, a given natural process in terms of the “interactions” among its many elementary units¹²⁻¹⁴. The principle of maximum entropy, by itself, contains the essence of the so-called Inverse Ising Problem solution, in which a “Hamiltonian” associated with a given complex system can be inferred from observed statistical correlations among its components. This method is frequently referred to as the Boltzmann machine, since the Boltzmann distribution is present at its core. Generally speaking, the MEM has been applied to systems that can be mapped to Ising-like models, namely, models in which the interacting elements are in an active or inactive state, i.e., a lattice of dipole moments with spins states that are either up or down under the action of an external field and their mutual interactions. For instance, when dealing with neuronal networks, the interactions between neurons reacting to some stimuli are inferred from their firing patterns¹⁵⁻²¹, in which the firing of a spike equate to the spin up state. This strategy has also been successful in characterizing protein-protein interactions²²,²³, the genetic interaction networks from gene expression patterns²⁴⁻²⁶, the collective responses exhibited by flocks of birds²⁷,²⁸,

<footnote n="1">¹Departamento de Física, Universidade Federal do Ceará, Fortaleza 60451-970, Brazil. ²Department of Mathematics and Physics, University of Campania “Luigi Vanvitelli”, 81100 Caserta, Italy. ³PMMH, ESPCI, CNRS UMR 7636, 7 Quai St. Bernard, 75005 Paris, France. ⁴Section on Critical Brain Dynamics, NIMH, Bethesda, MD 20892, USA. ⁵email: cesar@fisica.ufc.br</footnote>
and the emergence of collective behavior from the eye movement patterns of a group of people while watching commercial videos²⁹ or reading texts³⁰.

## Methods

Our initial goal here is to investigate the experimental data reported in Ref.⁶ on the transient synchronization of local neuronal groups, i.e., “spikes”, recorded with microelectrode arrays in neuronal cultures, but from the perspective of a pairwise Maximum-Entropy method. As we show next, we follow the approach proposed in Refs.¹⁶,³¹⁻³⁶ to build an Ising-like model that reproduces the same time-averaged spiking rates and pairwise correlations as the experimental data set, and then evaluate its behaviour in the analogous framework of the corresponding thermal equilibrium properties. Accordingly, let the variable $s_i(t)$ be the binarized time series of spikes of neuron i observed at a discretized time $t = 1, 2, \dots, M$, such that $s_i = +1$ if it fires and $s_i = -1$ if it does not fire. From these series, we can calculate the time-averaged activity for each neuron i,

$$ \langle s_i \rangle^{\text{obs}} = \frac{1}{M} \sum_{t=1}^{M} s_i(t), \qquad (1) $$

and the covariances between the spiking sequences of neurons i and j

$$ C_{ij}^{\text{obs}} = \langle s_i s_j \rangle^{\text{obs}} - \langle s_i \rangle^{\text{obs}} \langle s_j \rangle^{\text{obs}}, \qquad (2) $$

where $\langle s_i s_j \rangle = \frac{1}{M} \sum_{t=1}^{M} s_i(t)s_j(t)$. Moreover, to model these experimental sequences, we consider that $s_i$ corresponds to Ising-like variables on a fully connected network of N sites. Therefore, $\{s\} = \{s_1(t), \dots, s_N(t)\}$ can describe the system's state at a given time $t$. The probability distribution $P(\{s\})$ with the smallest number of parameters that represents our system is the one that maximizes the entropy, while reproducing our observations, i.e., $\langle s_i \rangle^{\text{obs}}$ for all N neurons and all $N(N-1)/2$ pairs of $C_{ij}^{\text{obs}}$. Given these constraints, the form of $P(\{s\})$ is the Boltzmann's probability distribution,

$$ P(\{s\}) \sim e^{-\mathcal{H}/T}, \qquad (3) $$

where $T$ is analogous to a temperature, $\mathcal{H}$ to a Hamiltonian³², and we have set $k_B = 1$. This distribution corresponds to the least biased representation for an Ising-like system like ours, with known first and second moments. Specifically, as a first approximation, the energy term has the same form of the Ising model,

$$ \mathcal{H} = - \sum_{i=1}^{N} h_i s_i - \sum_{i>j}^{N} J_{ij} s_i s_j. \qquad (4) $$

This mathematical correspondence naturally leads us to interpret $h_i$ as the action of a local external stimulus on neuron i, analogous to a “random field”, and $J_{ij}$ as a “coupling coefficient” between neurons i and j. Such pairwise couplings or interactions between the neuronal activities give rise to the observed correlations among them. At this point, we compute the local fields $h_i$ and the interactions $J_{ij}$ by directly solving the inverse problem given by Eq. (4). For simplicity, here we arbitrarily set the “operating temperature” to $T_o = 1$. On their turn, the local fields $h_i$ and interaction constants $J_{ij}$ are obtained through the following iteration scheme:

$$ J_{ij}(n+1) = J_{ij}(n) - \eta(n) [C_{ij}^{\text{MC}} - C_{ij}^{\text{obs}}], \qquad (5) $$

$$ h_i(n+1) = h_i(n) - \eta(n) [\langle s_i \rangle^{\text{MC}} - \langle s_i \rangle^{\text{obs}}], \qquad (6) $$

where $n$ is the iteration parameter and we start with $n = 1$ and $h_i(n = 1) = 0$. The covariance $C_{ij}^{\text{MC}}$ between two sites i and j of the Ising network of Eq. (4) is given by $C_{ij}^{\text{MC}} = \langle s_i s_j \rangle^{\text{MC}} - \langle s_i \rangle^{\text{MC}} \langle s_j \rangle^{\text{MC}}$, where the statistical average $\langle \dots \rangle^{\text{MC}}$ is obtained by performing a Monte Carlo simulation of the model Eq. (4) at temperature $T_o = 1$ using $h_i(n)$ and $J_{ij}(n)$. The function $\eta(n)$ is a learning rate which decays like $1/n^{0.437}$. Typically, we iterate till $n = 80,000$. Once we infer the values of $h_i$ and $J_{ij}$ that better reproduce the experimentally observed time-averaged activities $\langle s_i \rangle^{\text{obs}}$ and covariances $C_{ij}$, while maximizing the entropy, the Boltzmann probability distribution of Eq. (3) characterizes the statistics of the in vitro datasets.

We analyzed six in vitro samples (called “1” to “6”) consisting of 60 time series of binarized electrode spikes of 1 s divided in 20, 000 time bin each⁶ (see Figs. 1a-c, S1a, S1b, and S1c, for samples 1, 2, 3, 4, 5, and 6, respectively). They were obtained from coronal slices from rat dorsolateral cortex grown at 35.5 °C for 4-6 weeks before recording, as reported in Ref.⁶.

## Results and discussion

### In vitro sample characterization

In Fig. 2a-c we show the probability $P(K)$ that K electrodes fire simultaneously for in vitro sample 1, 2, and 3, respectively. The data follow an exponential decay $P(K) = \alpha e^{-\beta K}$, whose parameters depend on the sample considered. From the least-squares fit to the data points we obtained $\alpha = 0.50$ and $\beta = 0.35$ for sample 1, $\alpha = 0.21$ and $\beta = 0.25$ for sample 2, and $\alpha = 0.31$ and $\beta = 0.30$ for sample 3. The exponential fittings were confirmed by the Kolmogorov-Smirnov (KS) test, yielding a $p$ value > 0.05 in all cases. Therefore, the data exhibit a Poisson-like distribution, which would suggest uncorrelated events. We calculate next for every pair of sites {$ij$} the
and the emergence of collective behavior from the eye movement patterns of a group of people while watching commercial videos²⁹ or reading texts³⁰.

## Methods

Our initial goal here is to investigate the experimental data reported in Ref.⁶ on the transient synchronization of local neuronal groups, i.e., “spikes”, recorded with microelectrode arrays in neuronal cultures, but from the perspective of a pairwise Maximum-Entropy method. As we show next, we follow the approach proposed in Refs.¹⁶,³¹⁻³⁶ to build an Ising-like model that reproduces the same time-averaged spiking rates and pairwise correlations as the experimental data set, and then evaluate its behaviour in the analogous framework of the corresponding thermal equilibrium properties. Accordingly, let the variable $s_i(t)$ be the binarized time series of spikes of neuron $i$ observed at a discretized time $t = 1, 2, \dots, M$, such that $s_i = +1$ if it fires and $s_i = -1$ if it does not fire. From these series, we can calculate the time-averaged activity for each neuron $i$,

$$ \langle s_i \rangle^{\text{obs}} = \frac{1}{M} \sum_{t=1}^{M} s_i(t), \qquad (1) $$

and the covariances between the spiking sequences of neurons $i$ and $j$

$$ C_{ij}^{\text{obs}} = \langle s_i s_j \rangle^{\text{obs}} - \langle s_i \rangle^{\text{obs}} \langle s_j \rangle^{\text{obs}}, \qquad (2) $$

where $\langle s_i s_j \rangle = \frac{1}{M} \sum_{t=1}^{M} s_i(t)s_j(t)$. Moreover, to model these experimental sequences, we consider that $s_i$ corresponds to Ising-like variables on a fully connected network of $N$ sites. Therefore, $\{s\} = \{s_1(t), \dots, s_N(t)\}$ can describe the system's state at a given time $t$. The probability distribution $P(\{s\})$ with the smallest number of parameters that represents our system is the one that maximizes the entropy, while reproducing our observations, i.e., $\langle s_i \rangle^{\text{obs}}$ for all $N$ neurons and all $N(N-1)/2$ pairs of $C_{ij}^{\text{obs}}$. Given these constraints, the form of $P(\{s\})$ is the Boltzmann's probability distribution,

$$ P(\{s\}) \sim e^{-\mathcal{H}/T}, \qquad (3) $$

where $T$ is analogous to a temperature, $\mathcal{H}$ to a Hamiltonian³², and we have set $k_B = 1$. This distribution corresponds to the least biased representation for an Ising-like system like ours, with known first and second moments. Specifically, as a first approximation, the energy term has the same form of the Ising model,

$$ \mathcal{H} = -\sum_{i=1}^{N} h_i s_i - \sum_{i>j}^{N} J_{ij} s_i s_j. \qquad (4) $$

This mathematical correspondence naturally leads us to interpret $h_i$ as the action of a local external stimulus on neuron $i$, analogous to a “random field”, and $J_{ij}$ as a “coupling coefficient” between neurons $i$ and $j$. Such pairwise couplings or interactions between the neuronal activities give rise to the observed correlations among them. At this point, we compute the local fields $h_i$ and the interactions $J_{ij}$ by directly solving the inverse problem given by Eq. (4). For simplicity, here we arbitrarily set the “operating temperature” to $T_o = 1$. On their turn, the local fields $h_i$ and interaction constants $J_{ij}$ are obtained through the following iteration scheme:

$$ J_{ij}(n+1) = J_{ij}(n) - \eta(n) [C_{ij}^{\text{MC}} - C_{ij}^{\text{obs}}], \qquad (5) $$

$$ h_i(n+1) = h_i(n) - \eta(n) [\langle s_i \rangle^{\text{MC}} - \langle s_i \rangle^{\text{obs}}], \qquad (6) $$

where $n$ is the iteration parameter and we start with $n = 1$ and $h_i(n = 1) = 0$. The covariance $C_{ij}^{\text{MC}}$ between two sites $i$ and $j$ of the Ising network of Eq. (4) is given by $C_{ij}^{\text{MC}} = \langle s_i s_j \rangle^{\text{MC}} - \langle s_i \rangle^{\text{MC}} \langle s_j \rangle^{\text{MC}}$, where the statistical average $\langle \dots \rangle^{\text{MC}}$ is obtained by performing a Monte Carlo simulation of the model Eq. (4) at temperature $T_o = 1$ using $h_i(n)$ and $J_{ij}(n)$. The function $\eta(n)$ is a learning rate which decays like $1/n^{0.437}$. Typically, we iterate till $n = 80,000$. Once we infer the values of $h_i$ and $J_{ij}$ that better reproduce the experimentally observed time-averaged activities $\langle s_i \rangle^{\text{obs}}$ and covariances $C_{ij}$, while maximizing the entropy, the Boltzmann probability distribution of Eq. (3) characterizes the statistics of the in vitro datasets.

We analyzed six in vitro samples (called “1” to “6”) consisting of 60 time series of binarized electrode spikes of 1 s divided in 20, 000 time bin each⁶ (see Figs. 1a-c, S1a, S1b, and S1c, for samples 1, 2, 3, 4, 5, and 6, respectively). They were obtained from coronal slices from rat dorsolateral cortex grown at 35.5 °C for 4-6 weeks before recording, as reported in Ref.⁶.

## Results and discussion
### In vitro sample characterization

In Fig. 2a-c we show the probability $P(K)$ that $K$ electrodes fire simultaneously for in vitro sample 1, 2, and 3, respectively. The data follow an exponential decay $P(K) = \alpha e^{-\beta K}$, whose parameters depend on the sample considered. From the least-squares fit to the data points we obtained $\alpha = 0.50$ and $\beta = 0.35$ for sample 1, $\alpha = 0.21$ and $\beta = 0.25$ for sample 2, and $\alpha = 0.31$ and $\beta = 0.30$ for sample 3. The exponential fittings were confirmed by the Kolmogorov-Smirnov (KS) test, yielding a $p$ value > 0.05 in all cases. Therefore, the data exhibit a Poisson-like distribution, which would suggest uncorrelated events. We calculate next for every pair of sites {$ij$} the
<figure><img src="image_3.png" /><figcaption>**Figure 1.** The ensemble of investigated time series of binarized spikes spanning a total of 1 s for 60 electrodes of the in vitro samples 1, 2 and 3 are shown in (a–c)⁶, respectively. Each horizontal line represents the temporal spiking series of a local group at an electrode or single neuron, respectively, whereas each vertical line represents the state of the system at time *t*. A blue point corresponds to a spike (+1) and an empty place to no spike (−1).</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>**Figure 2.** The distribution $P(K)$ of the measured $K$ simultaneous spikes for the in vitro samples 1, 2, and 3 are shown in (a–c), respectively. The solid black line in each figure is the least-squares fit to the data points of the exponential function $p(K) = \alpha e^{-\beta K}$, with (a) $\alpha = 0.50$ and $\beta = 0.35$ for sample 1, (b) $\alpha = 0.21$ and $\beta = 0.25$ for sample 2, and (c) $\alpha = 0.31$ and $\beta = 0.30$ for sample 3. The exponential fittings were confirmed by the Kolmogorov-Smirnov (KS) test, yielding a $p$ value > 0.05 in all cases. In (d–f) are the distributions of correlation coefficients $C_{ij}^{\text{obs}}$ between pairs of electrodes for the in vitro samples 1, 2 and 3, respectively.</figcaption></figure>

corresponding time-averaged correlation $C_{ij}^{\text{obs}}$. Their distributions are shown in Fig. 2d–f for the in vitro samples 1, 2, and 3, respectively. All three in vitro samples have their maximum around zero, but with skewness towards positive values. Moreover, the in vitro distributions are qualitatively different from each other. The question that naturally arises is how this difference will be expressed by the Ising-like models that we construct next.

### Distributions of the learned parameters of the model in Eq. (4)

In order to apply the Boltzmann machine, we solved Eqs. (5) and (6) simultaneously to calculate for each sample the local fields $h_i$ and the coupling constants $J_{ij}$. Their distributions are shown in Fig. 3. We see that, for all samples, the fields $h_i$ are mostly negative but the most likely values depend on the sample considered. For samples 1 (Fig. 3a) and 2 (Fig. 3b) these are close to zero on the negative side while for sample 3 (Fig. 3c) the most likely value is close to zero on the negative side. Furthermore, all distributions exhibit a skewness towards negative fields. Since our system has only 60 electrodes, the sample to sample variations are considerable. As also shown in Fig. 3d–f the distributions of the interaction constants $J_{ij}$ are in all cases symmetrically centered around zero. Moreover, from the Kolmogorov-Smirnov (KS) test, Gaussian fittings appropriately describe the data for all samples.
<figure><img src="image_3.png" /><figcaption>**Figure 1.** The ensemble of investigated time series of binarized spikes spanning a total of 1 s for 60 electrodes of the in vitro samples 1, 2 and 3 are shown in (a–c)⁶, respectively. Each horizontal line represents the temporal spiking series of a local group at an electrode or single neuron, respectively, whereas each vertical line represents the state of the system at time *t*. A blue point corresponds to a spike (+1) and an empty place to no spike (−1).</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>**Figure 2.** The distribution $P(K)$ of the measured $K$ simultaneous spikes for the in vitro samples 1, 2, and 3 are shown in (a–c), respectively. The solid black line in each figure is the least-squares fit to the data points of the exponential function $p(K) = \alpha e^{-\beta K}$, with (a) $\alpha = 0.50$ and $\beta = 0.35$ for sample 1, (b) $\alpha = 0.21$ and $\beta = 0.25$ for sample 2, and (c) $\alpha = 0.31$ and $\beta = 0.30$ for sample 3. The exponential fittings were confirmed by the Kolmogorov-Smirnov (KS) test, yielding a $p$ value > 0.05 in all cases. In (d–f) are the distributions of correlation coefficients $C_{ij}^{\text{obs}}$ between pairs of electrodes for the in vitro samples 1, 2 and 3, respectively.</figcaption></figure>

corresponding time-averaged correlation $C_{ij}^{\text{obs}}$. Their distributions are shown in Fig. 2d–f for the in vitro samples 1, 2, and 3, respectively. All three in vitro samples have their maximum around zero, but with skewness towards positive values. Moreover, the in vitro distributions are qualitatively different from each other. The question that naturally arises is how this difference will be expressed by the Ising-like models that we construct next.

### Distributions of the learned parameters of the model in Eq. (4)

In order to apply the Boltzmann machine, we solved Eqs. (5) and (6) simultaneously to calculate for each sample the local fields $h_i$ and the coupling constants $J_{ij}$. Their distributions are shown in Fig. 3. We see that, for all samples, the fields $h_i$ are mostly negative but the most likely values depend on the sample considered. For samples 1 (Fig. 3a) and 2 (Fig. 3b) these are close to zero on the negative side while for sample 3 (Fig. 3c) the most likely value is close to zero on the negative side. Furthermore, all distributions exhibit a skewness towards negative fields. Since our system has only 60 electrodes, the sample to sample variations are considerable. As also shown in Fig. 3d–f the distributions of the interaction constants $J_{ij}$ are in all cases symmetrically centered around zero. Moreover, from the Kolmogorov-Smirnov (KS) test, Gaussian fittings appropriately describe the data for all samples.
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Figures (a–c) show the distributions of local fields $h_i$ after learning for in vitro samples 1, 2, and 3, respectively. All distributions display skewness towards negative values. Figures (d–f) show the distributions of the coupling coefficients $J_{ij}$ between $N = 60$ electrodes for in vitro samples. The black continuous lines in (d–f) correspond to a Gaussian fitting and the adequacy of this model was confirmed by the Kolmogorov–Smirnov (KS) test.</figcaption></figure>

## Test to verify how well the Ising-like model of Eq. (4) can reproduce the measured time averages

To verify how well the Ising-like model Eq. (4) is able to reproduce the measured time averages of the experimental data, we show in Fig. 4a–c $\langle s_i \rangle^{MC}$ against $\langle s_i \rangle^{obs}$ for every site $i$ of the in vitro sample 1, 2, and 3, respectively, and in Figs. S2a, S2b, and S2c for the samples 4, 5 and 6, respectively. The same is shown in Fig. 4d–f but for the two-point correlations $C_{ij}^{MC}$ against $C_{ij}^{obs}$ for every pair of sites $i$ and $j$ of the samples 1, 2, and 3, respectively, and in Figs. S3a, S3b, and S3c for the samples 4, 5 and 6, respectively. Clearly, the agreement is excellent in all cases. For samples 1, 2, and 3, as illustrated in Fig. 5a–c, respectively, the model demonstrates remarkable accuracy in predicting the experimentally observed correlations among triplets of neurons. The precision of these predictions is quantitatively supported by the relative errors, $\langle (T_{ijk}^{MC} - T_{ijk}^{obs})/T_{ijk}^{obs} \rangle$, which amount to 0.095%, 0.091%, and 0.120% for samples 1, 2, and 3, respectively. These deviations are unusually smaller than those observed in previous studies¹⁶, highlighting the adequacy of the Ising-like model to describe the time-averaged spiking behavior of these in vitro neuronal networks. Moreover, since the three-point correlations were not explicitly utilized in the correlation function matching procedure, these findings strongly indicate that the machine demonstrates a robust level of generalization for these samples, providing us with a solid basis to confidently infer thermodynamic properties from the learned parameters. In the case of samples 4, 5, and 6, as depicted in Figs. S4a, S4b, and S4c, respectively, the models accurately reflect the rankings of the triplets, in the sense that triplets observed experimentally to be more strongly correlated are also predicted to be relatively more strongly correlated within the set of model predictions. Nonetheless, the models' inability to quantitatively replicate the observed three-point correlations prevents our examination of their thermodynamic properties.

## Critical properties of the Ising-like model

After showing that the Ising-like model describes to a satisfactory degree the time-averaged properties of our experimental samples, we now proceed with a more detailed investigation of its properties in terms of the implicit thermodynamic analogy brought about naturally³³. The model adopted here is in fact an extension of the famous Sherrington-Kirkpatrick (SK) spin glass¹² without averaging over quenched disorder, but including a random field. Since the coupling constants of our model are distributed around zero, we expect frustration effects and, due to the long-range interactions, the occurrence of mean-field behaviour. In addition to the SK spin glass, however, we have a predominantly negative quenched random field. This extension of the SK model was in fact studied by Hadjiagapiou¹³ for the case in which both the interaction constants and the local fields are Gaussian distributed. We assume that, despite some differences in the distributions, our model should have a similar phase diagram as in Ref.¹³, namely, as a function of temperature, either a spin glass to paramagnetic transition for weak negative fields (below the de Almeida-Thouless (AT) line¹⁴) or a ferromagnetic-paramagnetic phase
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Figures (a–c) show the distributions of local fields $h_i$ after learning for in vitro samples 1, 2, and 3, respectively. All distributions display skewness towards negative values. Figures (d–f) show the distributions of the coupling coefficients $J_{ij}$ between $N = 60$ electrodes for in vitro samples. The black continuous lines in (d–f) correspond to a Gaussian fitting and the adequacy of this model was confirmed by the Kolmogorov–Smirnov (KS) test.</figcaption></figure>

## Test to verify how well the Ising-like model of Eq. (4) can reproduce the measured time averages

To verify how well the Ising-like model Eq. (4) is able to reproduce the measured time averages of the experimental data, we show in Fig. 4a–c $\langle s_i \rangle^{MC}$ against $\langle s_i \rangle^{obs}$ for every site $i$ of the in vitro sample 1, 2, and 3, respectively, and in Figs. S2a, S2b, and S2c for the samples 4, 5 and 6, respectively. The same is shown in Fig. 4d–f but for the two-point correlations $C_{ij}^{MC}$ against $C_{ij}^{obs}$ for every pair of sites $i$ and $j$ of the samples 1, 2, and 3, respectively, and in Figs. S3a, S3b, and S3c for the samples 4, 5 and 6, respectively. Clearly, the agreement is excellent in all cases. For samples 1, 2, and 3, as illustrated in Fig. 5a–c, respectively, the model demonstrates remarkable accuracy in predicting the experimentally observed correlations among triplets of neurons. The precision of these predictions is quantitatively supported by the relative errors, $\langle (T_{ijk}^{MC} - T_{ijk}^{obs})/T_{ijk}^{obs} \rangle$, which amount to 0.095%, 0.091%, and 0.120% for samples 1, 2, and 3, respectively. These deviations are unusually smaller than those observed in previous studies¹⁶, highlighting the adequacy of the Ising-like model to describe the time-averaged spiking behavior of these in vitro neuronal networks. Moreover, since the three-point correlations were not explicitly utilized in the correlation function matching procedure, these findings strongly indicate that the machine demonstrates a robust level of generalization for these samples, providing us with a solid basis to confidently infer thermodynamic properties from the learned parameters. In the case of samples 4, 5, and 6, as depicted in Figs. S4a, S4b, and S4c, respectively, the models accurately reflect the rankings of the triplets, in the sense that triplets observed experimentally to be more strongly correlated are also predicted to be relatively more strongly correlated within the set of model predictions. Nonetheless, the models' inability to quantitatively replicate the observed three-point correlations prevents our examination of their thermodynamic properties.

## Critical properties of the Ising-like model

After showing that the Ising-like model describes to a satisfactory degree the time-averaged properties of our experimental samples, we now proceed with a more detailed investigation of its properties in terms of the implicit thermodynamic analogy brought about naturally³³. The model adopted here is in fact an extension of the famous Sherrington-Kirkpatrick (SK) spin glass¹² without averaging over quenched disorder, but including a random field. Since the coupling constants of our model are distributed around zero, we expect frustration effects and, due to the long-range interactions, the occurrence of mean-field behaviour. In addition to the SK spin glass, however, we have a predominantly negative quenched random field. This extension of the SK model was in fact studied by Hadjiagapiou¹³ for the case in which both the interaction constants and the local fields are Gaussian distributed. We assume that, despite some differences in the distributions, our model should have a similar phase diagram as in Ref.¹³, namely, as a function of temperature, either a spin glass to paramagnetic transition for weak negative fields (below the de Almeida-Thouless (AT) line¹⁴) or a ferromagnetic-paramagnetic phase
<figure><img src="image_6.png" /><figcaption>Figure 4. In (a–c) are shown the magnetizations $\langle s_i \rangle^{MC}$ versus $\langle s_i \rangle^{obs}$ for the in vitro samples 1, 2, and 3, respectively. The same in (d–f), but for the correlations $C_{ij}^{MC}$ versus $C_{ij}^{obs}$ obtained from the in vitro samples 1, 2, and 3, respectively. The error bars represent the standard deviations calculated from $10^5$ samples generated through repeated Monte Carlo runs. The solid black lines correspond to the function $y = x$.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Figure 5. Test to verify if the three-point correlations $T_{ijk}^{MC}$ from the Monte Carlo simulations of the Ising-like model of Eq. (4) are able to reproduce the three-point correlations $T_{ijk}^{obs}$ of the experimental in vitro samples 1, 2 and 3, in (a–c), respectively. The solid black lines correspond to the function $y = x$. The triplets are binned into 100 populated bins and the error bars are the standard deviations across the bins.</figcaption></figure>

transition for stronger negative fields (above the AT line). Moreover, the upper critical dimension should be that of the random field Ising model³⁸.

Under this analogous thermodynamic framework, for every spin configuration {$s_i$}, one can define its magnetization $M(T) = (1/N) \sum_i s_i$ and its energy $E(T) = (1/N) \sum_{ij} J_{ij} s_i s_j + (1/N) \sum_i h_i s_i$. The order parameter for a ferromagnetic phase transition is given by the average magnetization $\langle M(T) \rangle^{MC}$, and the fluctuation-dissipation theorem provides the following expression for the specific heat:

$$C(T) = \frac{1}{T^2} \left( \langle E(T)^2 \rangle^{MC} - \langle (E(T))^{MC} \rangle^2 \right), \qquad (7)$$

and the susceptibility

$$\chi(T) = \frac{1}{T} \left( \langle M(T)^2 \rangle^{MC} - \langle (M(T))^{MC} \rangle^2 \right). \qquad (8)$$

In Fig. 6a the ferromagnetic order parameter $\langle M(T) \rangle^{MC}$ is plotted against temperature T for the three in vitro samples. We have verified that these curves do not depend on the initial configuration used in the Monte Carlo
<figure><img src="image_6.png" /><figcaption>Figure 4. In (a–c) are shown the magnetizations $\langle s_i \rangle^{MC}$ versus $\langle s_i \rangle^{obs}$ for the in vitro samples 1, 2, and 3, respectively. The same in (d–f), but for the correlations $C_{ij}^{MC}$ versus $C_{ij}^{obs}$ obtained from the in vitro samples 1, 2, and 3, respectively. The error bars represent the standard deviations calculated from $10^5$ samples generated through repeated Monte Carlo runs. The solid black lines correspond to the function $y = x$.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Figure 5. Test to verify if the three-point correlations $T_{ijk}^{MC}$ from the Monte Carlo simulations of the Ising-like model of Eq. (4) are able to reproduce the three-point correlations $T_{ijk}^{obs}$ of the experimental in vitro samples 1, 2 and 3, in (a–c), respectively. The solid black lines correspond to the function $y = x$. The triplets are binned into 100 populated bins and the error bars are the standard deviations across the bins.</figcaption></figure>

transition for stronger negative fields (above the AT line). Moreover, the upper critical dimension should be that of the random field Ising model³⁸.

Under this analogous thermodynamic framework, for every spin configuration $\{s_i\}$, one can define its magnetization $M(T) = (1/N) \sum_i s_i$ and its energy $E(T) = (1/N) \sum_{ij} J_{ij} s_i s_j + (1/N) \sum_i h_i s_i$. The order parameter for a ferromagnetic phase transition is given by the average magnetization $\langle M(T) \rangle^{MC}$, and the fluctuation-dissipation theorem provides the following expression for the specific heat:

$$C(T) = \frac{1}{T^2} \left( \langle E(T)^2 \rangle^{MC} - \langle (E(T))^{MC} \rangle^2 \right), \quad (7)$$

and the susceptibility

$$\chi(T) = \frac{1}{T} \left( \langle M(T)^2 \rangle^{MC} - \langle (M(T))^{MC} \rangle^2 \right). \quad (8)$$

In Fig. 6a the ferromagnetic order parameter $\langle M(T) \rangle^{MC}$ is plotted against temperature $T$ for the three in vitro samples. We have verified that these curves do not depend on the initial configuration used in the Monte Carlo
<figure><img src="image_8.png" /><figcaption>**Figure 6.** Thermodynamic analysis for the in vitro samples 1, 2, and 3. In (a) is shown the magnetization $M(T)$ as function of temperature $T$. In (b) is shown the specific heat $C(T)$ as a function of temperature $T$. We also show the specific heat for sample 1 after reshuffling the data (solid black line). In (c) is shown the susceptibility $\chi(T)$ against temperature $T$. In all cases, the dashed line corresponds to the operation temperature $T_o = 1$.</figcaption></figure>

simulation. The fact that we do not observe hysteresis indicates that the transition might be second order and not first order as observed, for example, in retinal ganglion cell spike trains<sup>21</sup>. Figure 6b shows the specific heat $C(T)$ as a function of temperature for the three in vitro samples. The curves attain maxima at a temperature $T_c^*$ that depends on the sample considered, with $T_c^* = 1.18$ for sample 1, $T_c^* = 1.20$ for sample 2, and $T_c^* = 1.27$ for sample 3. Also shown in Fig. 6b is the dependence of the specific heat on temperature obtained by randomly shuffling the spikes in the time series of the in vitro sample 1, i.e., by suppressing the intrinsic correlations present in the spike sequence. We see that, in this case, the sharp maximum of the specific heat is dramatically attenuated, which means that there is no ferromagnetic phase transition any more as a function of temperature. Possibly, this case falls below the AT line and might exhibit a spin glass phase at low temperatures. Higher order correlations seem important for neurons that are close to each other<sup>36</sup>. In order to take those into account, the method can be generalized to multi-spin interactions<sup>35</sup>. We also studied the susceptibilities $\chi(T)$ of the Ising-like models for samples 1,2, and 3, shown in Fig. 6c. They also exhibit a maximum at a given temperature $T_c^*$, but unlike specific heat curves, the maxima of the susceptibility occur at the same value $T_c^* = 1.34$. Furthermore, the maxima of the susceptibility $\chi_{max}$ are quite different between the samples, despite having the same number of spins, with $\chi_{max}^{(3)} > \chi_{max}^{(2)} > \chi_{max}^{(1)}$ for sample 3, 2, and 1, respectively.

## Conclusion

In summary, we analyzed time-averaged experimental data from the activities of in vitro neuronal networks. By implementing a Pairwise Maximum-Entropy method alongside an inverse binary Ising-like model, we successfully identified the local fields and interaction couplings that best represent the average neuronal activities and the statistical correlations between neuronal pairs. This approach confirmed the model's effectiveness in capturing the inherent complexities of neuronal network dynamics, where specific information about neuron types is reflected in the local fields, and a symmetric distribution in interaction constants emerges as a characteristic trait of the network's interaction framework.

The investigation into the thermodynamic properties of these models revealed typical features of second-order phase transitions between ferromagnetic and paramagnetic states at temperatures $T_c^*$ above, but not far from unity, which corresponds to the operating temperature, $T_o = 1$, utilized in the Maximum-Entropy method. This finding extends the thermodynamic analogy within our Ising-like models, giving support to the hypothesis of criticality in neuronal network behavior. It points out the relevance of critical phenomena in understanding the functionality and complexity of neuronal networks.

The exceptional alignment between our model predictions and experimental observations, particularly in the accurate rendering of three-point correlations, demonstrates a case of success of the methodology. This precision emphasizes the Ising-like model's adequacy in representing the nuanced spiking behavior of neuronal networks, surpassing previous modeling attempts in terms of accuracy and reliability. These findings not only validate the criticality hypothesis in neuronal networks, but also highlight the potential of using thermodynamic and statistical physics concepts to explore neural computation and information processing.

## Data availability

The source data for all figures in this study are provided for this paper. Source data are provided with this paper. The in vitro data is available at https://github.com/cesampaiof/Ising-like-model-neuronal-networks.

Received: 28 September 2023; Accepted: 28 February 2024

Published online: 25 March 2024

## References

1. Amit, D. J., Gutfreund, H. & Sompolinky, H. Spin-glass models of neural networks. *Phys. Rev. A* **32**, 1007 (1985).

2. Hopfield, J. J. Neural networks and physical systems with emergent collective computational abilities. *Proc. Natl. Acad. Sci. U.S.A.* **79**, 2554 (1982).

3. Amit, D. J. *Modeling Brain Function: The World of Attractor Neural Networks* (Cambridge Press, Cambridge, 1989).

4. Miranda, E. N. & Herrmann, H. J. Self-organized criticality with disorder and frustration. *Physica A* **175**, 339 (1991).
<figure><img src="image_8.png" /><figcaption>**Figure 6.** Thermodynamic analysis for the in vitro samples 1, 2, and 3. In (a) is shown the magnetization $M(T)$ as function of temperature $T$. In (b) is shown the specific heat $C(T)$ as a function of temperature $T$. We also show the specific heat for sample 1 after reshuffling the data (solid black line). In (c) is shown the susceptibility $\chi(T)$ against temperature $T$. In all cases, the dashed line corresponds to the operation temperature $T_o = 1$.</figcaption></figure>

simulation. The fact that we do not observe hysteresis indicates that the transition might be second order and not first order as observed, for example, in retinal ganglion cell spike trains<sup>21</sup>. Figure 6b shows the specific heat $C(T)$ as a function of temperature for the three in vitro samples. The curves attain maxima at a temperature $T_c^*$ that depends on the sample considered, with $T_c^* = 1.18$ for sample 1, $T_c^* = 1.20$ for sample 2, and $T_c^* = 1.27$ for sample 3. Also shown in Fig. 6b is the dependence of the specific heat on temperature obtained by randomly shuffling the spikes in the time series of the in vitro sample 1, i.e., by suppressing the intrinsic correlations present in the spike sequence. We see that, in this case, the sharp maximum of the specific heat is dramatically attenuated, which means that there is no ferromagnetic phase transition any more as a function of temperature. Possibly, this case falls below the AT line and might exhibit a spin glass phase at low temperatures. Higher order correlations seem important for neurons that are close to each other<sup>36</sup>. In order to take those into account, the method can be generalized to multi-spin interactions<sup>35</sup>. We also studied the susceptibilities $\chi(T)$ of the Ising-like models for samples 1,2, and 3, shown in Fig. 6c. They also exhibit a maximum at a given temperature $T_c^*$, but unlike specific heat curves, the maxima of the susceptibility occur at the same value $T_c^* = 1.34$. Furthermore, the maxima of the susceptibility $\chi_{max}$ are quite different between the samples, despite having the same number of spins, with $\chi_{max}^{(3)} > \chi_{max}^{(2)} > \chi_{max}^{(1)}$ for sample 3, 2, and 1, respectively.

## Conclusion

In summary, we analyzed time-averaged experimental data from the activities of in vitro neuronal networks. By implementing a Pairwise Maximum-Entropy method alongside an inverse binary Ising-like model, we successfully identified the local fields and interaction couplings that best represent the average neuronal activities and the statistical correlations between neuronal pairs. This approach confirmed the model's effectiveness in capturing the inherent complexities of neuronal network dynamics, where specific information about neuron types is reflected in the local fields, and a symmetric distribution in interaction constants emerges as a characteristic trait of the network's interaction framework.

The investigation into the thermodynamic properties of these models revealed typical features of second-order phase transitions between ferromagnetic and paramagnetic states at temperatures $T_c^*$ above, but not far from unity, which corresponds to the operating temperature, $T_o = 1$, utilized in the Maximum-Entropy method. This finding extends the thermodynamic analogy within our Ising-like models, giving support to the hypothesis of criticality in neuronal network behavior. It points out the relevance of critical phenomena in understanding the functionality and complexity of neuronal networks.

The exceptional alignment between our model predictions and experimental observations, particularly in the accurate rendering of three-point correlations, demonstrates a case of success of the methodology. This precision emphasizes the Ising-like model's adequacy in representing the nuanced spiking behavior of neuronal networks, surpassing previous modeling attempts in terms of accuracy and reliability. These findings not only validate the criticality hypothesis in neuronal networks, but also highlight the potential of using thermodynamic and statistical physics concepts to explore neural computation and information processing.

## Data availability

The source data for all figures in this study are provided for this paper. Source data are provided with this paper. The in vitro data is available at https://github.com/cesampaiof/Ising-like-model-neuronal-networks.

Received: 28 September 2023; Accepted: 28 February 2024

Published online: 25 March 2024

## References

1. Amit, D. J., Gutfreund, H. & Sompolinky, H. Spin-glass models of neural networks. *Phys. Rev. A* **32**, 1007 (1985).

2. Hopfield, J. J. Neural networks and physical systems with emergent collective computational abilities. *Proc. Natl. Acad. Sci. U.S.A.* **79**, 2554 (1982).

3. Amit, D. J. *Modeling Brain Function: The World of Attractor Neural Networks* (Cambridge Press, Cambridge, 1989).

4. Miranda, E. N. & Herrmann, H. J. Self-organized criticality with disorder and frustration. *Physica A* **175**, 339 (1991).
5. Chialvo, D. R. & Bak, P. Learning from mistakes. *Neuroscience* **90**, 1137 (1999).

6. Beggs, J. M. & Plenz, D. Neuronal avalanches in neocortical circuits. *J. Neurosci.* **23**, 11167 (2003).

7. de Arcangelis, L., Perrone Capano, C. & Herrmann, H. J. Self-organized criticality model for brain plasticity. *Phys. Rev. Lett.* **96**, 028107 (2006).

8. Levina, A., Herrmann, J. M. & Geisel, T. Dynamical synapses causing self-organized criticality in neural networks. *Nat. Phys.* **3**, 857 (2007).

9. Plenz, D. et al. Self-organized criticality in the brain. *Front. Phys.* **9**, 639389 (2021).

10. Beggs, J. M. & Timme, N. Being critical of criticality in the brain. *Front. Physiol.* **3**, 163 (2012).

11. Nandi, M. K., Sarracino, A., Herrmann, H. J. & de Arcangelis, L. Scaling of avalanche shape and activity power spectrum in neuronal networks. *Phys. Rev. E* **106**, 024304 (2022).

12. Sherrington, D. & Kirkpatrick, S. Solvable model of a spin-glass. *Phys. Rev. Lett.* **35**, 1792 (1975).

13. Hadjiagapiou, I. A. The Sherrington–Kirkpatrick spin glass model in the presence of a random field with a joint Gaussian probability density function for the exchange interactions and random fields. *Physica A* **397**, 1 (2014).

14. De Almeida, J. R. L. & Thouless, D. J. Stability of the Sherrington–Kirkpatrick solution of a spin glass model. *J. Phys. A: Math. Theor.* **11**, 983 (1978).

15. Cocco, S., Leibler, S. & Monasson, R. Neuronal couplings between retinal ganglion cells inferred by efficient inverse statistical physics methods. *Proc. Natl. Acad. Sci. U.S.A.* **106**, 14058–14062 (2009).

16. Tkacik, G., Schneidman, E., Berry II, M. J. & Bialek, W. Spin glass models for a network of real neurons (2009). arXiv:0912.5409v1.

17. Shlens, J. et al. The structure of multi-neuron firing patterns in primate retina. *J. Neurosci.* **26**, 8254–8266 (2006).

18. Tang, A. et al. A maximum entropy model applied to spatial and temporal correlations from cortical networks in vitro. *J. Neurosci.* **28**, 505–518 (2008).

19. Mora, T., Deny, S. & Marre, O. Dynamical criticality in the collective activity of a population of retinal neurons. *Phys. Rev. Lett.* **114**, 078105 (2015).

20. Lotfi, N. et al. Signatures of brain criticality unveiled by maximum entropy analysis across cortical states. *Phys. Rev. E* **102**, 012408 (2020).

21. Ioffe, M. L. & Berry, M. J. II. The structured ‘low temperature’ phase of the retinal population code. *PLoS Comput. Biol.* **13**, e1005792 (2017).

22. Morcos, F. et al. Direct-coupling analysis of residue coevolution captures native contacts across many protein families. *Proc. Natl. Acad. Sci. U.S.A.* **108**, E1293–E1301 (2011).

23. Weigt, M., White, R. A., Szurmant, H., Hoch, J. A. & Hwa, T. Identification of direct residue contacts in protein-protein interaction by message passing. *Proc. Natl. Acad. Sci. U.S.A.* **106**, 67–72 (2009).

24. Stein, R. R., Marks, D. S. & Sander, C. Inferring pairwise interactions from biological data using maximum-entropy probability models. *PLoS Comput. Biol.* **11**, e1004182 (2015).

25. Lezon, T. R., Banavar, J. R., Cieplak, M., Maritan, A. & Fedoroff, N. V. Using the principle of entropy maximization to infer genetic interaction networks from gene expression patterns. *Proc. Natl. Acad. Sci. U.S.A.* **103**, 19033–19038 (2006).

26. Locasale, J. W. & Wolf-Yadlin, A. Maximum entropy reconstructions of dynamic signaling networks from quantitative proteomics data. *PLoS ONE* **4**, e6522 (2009).

27. Bialek, W. et al. Statistical mechanics for natural flocks of birds. *Proc. Natl. Acad. Sci. U.S.A.* **109**, 4786–4791 (2012).

28. Bialek, W. et al. Social interactions dominate speed control in poising natural flocks near criticality. *Proc. Natl. Acad. Sci. U.S.A.* **111**, 7212–7217 (2014).

29. Burleson-Lesser, K., Morone, F., DeGuzman, P., Parra, L. C. & Makse, H. A. Collective behaviour in video viewing: A thermodynamic analysis of gaze position. *PLoS ONE* **12**, e0168995 (2017).

30. Torres, D. et al. Eye-tracking as a proxy for coherence and complexity of texts. *PLoS ONE* **16**, e0260236 (2021).

31. Schneidman, E., Berry, M. J. II., Segev, R. & Bialek, W. Weak pairwise correlations imply strongly correlated network states in a neural population. *Nature* **440**, 1007 (2006).

32. Tkacik, G. et al. Thermodynamics and signatures of criticality in a network of neurons. *Proc. Natl. Acad. Sci. U.S.A.* **112**, 11513 (2015).

33. Tkacik, G. et al. Searching for collective behavior in a large network of sensory neurons. *PLoS Comput. Biol.* **10**, e1003408 (2014).

34. Tkacik, G., Schneidman, E., Berry II, M. J. & Bialek, W. Ising models for networks of real neurons (2006). arXiv:q-bio/0611072.

35. Albert, J. & Swendsen, R. H. Detecting multi-spin interactions in the inverse Ising problem. *Physica A* **483**, 293 (2017).

36. Ohiorhenuan, I. E. et al. Sparse coding and high-order correlations in fine-scale cortical networks. *Nature* **466**, 617 (2010).

37. Nguyen, H. C., Zecchina, R. & Berg, J. Inverse statistical problems: From the inverse Ising problem to data science. *Adv. Phys.* **66**, 197–261 (2017).

38. Ahrens, B. & Hartmann, A. Critical behavior of the random-field Ising model at and beyond the upper critical dimension. *Phys. Rev. E* **83**, 014205 (2011).

## Acknowledgements

We thank the Brazilian agencies CNPq, CAPES, FUNCAP and the National Institute of Science and Technology for Complex Systems (INCT-SC) in Brazil for financial support. LdA acknowledges financial support from the MIUR PRIN 2017WZFTZP. LdA also thanks the support by #NEXTGENERATIONEU (NGEU) and funded by the Ministry of University and Research (MUR), National Recovery and Resilience Plan (NRRP), project MNESYS (PE0000006)—A Multiscale integrated approach to the study of the nervous system in health and disease (DN. 1553 11.10.2022). This research was supported by the Division of the Intramural Research Program (DIRP) of the National Institute of Mental Health (NIMH), USA, ZIAMH002797, ZIAMH002971, and the BRAIN initiative Grant U19 NS107464-01 (DP, PK, TLR).

## Author contributions

Conceptualization: C.I.N.S.F., L.A., H.J.H., D.P., and J.S.A.Jr.; Data Curation: D.P., P.K., T.L.R.; Formal Analysis: C.I.N.S.F., L.A., D.P., H.J.H., J.S.A.Jr.; Funding Acquisition: D.P. and J.S.A.Jr.; Project Administration: D.P. and J.S.A.Jr.; Writing—Original Draft: H.J.H. and J.S.A.Jr.; Writing—Review Draft: C.I.N.S.F., L.A., H.J.H., D.P., P.K., T.L.R., J.S.A.Jr.

## Competing interests

The authors declare no competing interests.
5. Chialvo, D. R. & Bak, P. Learning from mistakes. *Neuroscience* **90**, 1137 (1999).

6. Beggs, J. M. & Plenz, D. Neuronal avalanches in neocortical circuits. *J. Neurosci.* **23**, 11167 (2003).

7. de Arcangelis, L., Perrone Capano, C. & Herrmann, H. J. Self-organized criticality model for brain plasticity. *Phys. Rev. Lett.* **96**, 028107 (2006).

8. Levina, A., Herrmann, J. M. & Geisel, T. Dynamical synapses causing self-organized criticality in neural networks. *Nat. Phys.* **3**, 857 (2007).

9. Plenz, D. et al. Self-organized criticality in the brain. *Front. Phys.* **9**, 639389 (2021).

10. Beggs, J. M. & Timme, N. Being critical of criticality in the brain. *Front. Physiol.* **3**, 163 (2012).

11. Nandi, M. K., Sarracino, A., Herrmann, H. J. & de Arcangelis, L. Scaling of avalanche shape and activity power spectrum in neuronal networks. *Phys. Rev. E* **106**, 024304 (2022).

12. Sherrington, D. & Kirkpatrick, S. Solvable model of a spin-glass. *Phys. Rev. Lett.* **35**, 1792 (1975).

13. Hadjiagapiou, I. A. The Sherrington–Kirkpatrick spin glass model in the presence of a random field with a joint Gaussian probability density function for the exchange interactions and random fields. *Physica A* **397**, 1 (2014).

14. De Almeida, J. R. L. & Thouless, D. J. Stability of the Sherrington–Kirkpatrick solution of a spin glass model. *J. Phys. A: Math. Theor.* **11**, 983 (1978).

15. Cocco, S., Leibler, S. & Monasson, R. Neuronal couplings between retinal ganglion cells inferred by efficient inverse statistical physics methods. *Proc. Natl. Acad. Sci. U.S.A.* **106**, 14058–14062 (2009).

16. Tkacik, G., Schneidman, E., Berry II, M. J. & Bialek, W. Spin glass models for a network of real neurons (2009). arXiv:0912.5409v1.

17. Shlens, J. et al. The structure of multi-neuron firing patterns in primate retina. *J. Neurosci.* **26**, 8254–8266 (2006).

18. Tang, A. et al. A maximum entropy model applied to spatial and temporal correlations from cortical networks in vitro. *J. Neurosci.* **28**, 505–518 (2008).

19. Mora, T., Deny, S. & Marre, O. Dynamical criticality in the collective activity of a population of retinal neurons. *Phys. Rev. Lett.* **114**, 078105 (2015).

20. Lotfi, N. et al. Signatures of brain criticality unveiled by maximum entropy analysis across cortical states. *Phys. Rev. E* **102**, 012408 (2020).

21. Ioffe, M. L. & Berry, M. J. II. The structured ‘low temperature’ phase of the retinal population code. *PLoS Comput. Biol.* **13**, e1005792 (2017).

22. Morcos, F. et al. Direct-coupling analysis of residue coevolution captures native contacts across many protein families. *Proc. Natl. Acad. Sci. U.S.A.* **108**, E1293–E1301 (2011).

23. Weigt, M., White, R. A., Szurmant, H., Hoch, J. A. & Hwa, T. Identification of direct residue contacts in protein-protein interaction by message passing. *Proc. Natl. Acad. Sci. U.S.A.* **106**, 67–72 (2009).

24. Stein, R. R., Marks, D. S. & Sander, C. Inferring pairwise interactions from biological data using maximum-entropy probability models. *PLoS Comput. Biol.* **11**, e1004182 (2015).

25. Lezon, T. R., Banavar, J. R., Cieplak, M., Maritan, A. & Fedoroff, N. V. Using the principle of entropy maximization to infer genetic interaction networks from gene expression patterns. *Proc. Natl. Acad. Sci. U.S.A.* **103**, 19033–19038 (2006).

26. Locasale, J. W. & Wolf-Yadlin, A. Maximum entropy reconstructions of dynamic signaling networks from quantitative proteomics data. *PLoS ONE* **4**, e6522 (2009).

27. Bialek, W. et al. Statistical mechanics for natural flocks of birds. *Proc. Natl. Acad. Sci. U.S.A.* **109**, 4786–4791 (2012).

28. Bialek, W. et al. Social interactions dominate speed control in poising natural flocks near criticality. *Proc. Natl. Acad. Sci. U.S.A.* **111**, 7212–7217 (2014).

29. Burleson-Lesser, K., Morone, F., DeGuzman, P., Parra, L. C. & Makse, H. A. Collective behaviour in video viewing: A thermodynamic analysis of gaze position. *PLoS ONE* **12**, e0168995 (2017).

30. Torres, D. et al. Eye-tracking as a proxy for coherence and complexity of texts. *PLoS ONE* **16**, e0260236 (2021).

31. Schneidman, E., Berry, M. J. II., Segev, R. & Bialek, W. Weak pairwise correlations imply strongly correlated network states in a neural population. *Nature* **440**, 1007 (2006).

32. Tkacik, G. et al. Thermodynamics and signatures of criticality in a network of neurons. *Proc. Natl. Acad. Sci. U.S.A.* **112**, 11513 (2015).

33. Tkacik, G. et al. Searching for collective behavior in a large network of sensory neurons. *PLoS Comput. Biol.* **10**, e1003408 (2014).

34. Tkacik, G., Schneidman, E., Berry II, M. J. & Bialek, W. Ising models for networks of real neurons (2006). arXiv:q-bio/0611072.

35. Albert, J. & Swendsen, R. H. Detecting multi-spin interactions in the inverse Ising problem. *Physica A* **483**, 293 (2017).

36. Ohiorhenuan, I. E. et al. Sparse coding and high-order correlations in fine-scale cortical networks. *Nature* **466**, 617 (2010).

37. Nguyen, H. C., Zecchina, R. & Berg, J. Inverse statistical problems: From the inverse Ising problem to data science. *Adv. Phys.* **66**, 197–261 (2017).

38. Ahrens, B. & Hartmann, A. Critical behavior of the random-field Ising model at and beyond the upper critical dimension. *Phys. Rev. E* **83**, 014205 (2011).

## Acknowledgements

We thank the Brazilian agencies CNPq, CAPES, FUNCAP and the National Institute of Science and Technology for Complex Systems (INCT-SC) in Brazil for financial support. LdA acknowledges financial support from the MIUR PRIN 2017WZFTZP. LdA also thanks the support by #NEXTGENERATIONEU (NGEU) and funded by the Ministry of University and Research (MUR), National Recovery and Resilience Plan (NRRP), project MNESYS (PE0000006)—A Multiscale integrated approach to the study of the nervous system in health and disease (DN. 1553 11.10.2022). This research was supported by the Division of the Intramural Research Program (DIRP) of the National Institute of Mental Health (NIMH), USA, ZIAMH002797, ZIAMH002971, and the BRAIN initiative Grant U19 NS107464-01 (DP, PK, TLR).

## Author contributions

Conceptualization: C.I.N.S.F., L.A., H.J.H., D.P., and J.S.A.Jr.; Data Curation: D.P., P.K., T.L.R.; Formal Analysis: C.I.N.S.F., L.A., D.P., H.J.H., J.S.A.Jr.; Funding Acquisition: D.P. and J.S.A.Jr.; Project Administration: D.P. and J.S.A.Jr.; Writing—Original Draft: H.J.H. and J.S.A.Jr.; Writing—Review Draft: C.I.N.S.F., L.A., H.J.H., D.P., P.K., T.L.R., J.S.A.Jr.

## Competing interests

The authors declare no competing interests.
**Additional information**

**Supplementary Information** The online version contains supplementary material available at https://doi.org/10.1038/s41598-024-55922-9.

**Correspondence** and requests for materials should be addressed to C.I.N.S.F.

**Reprints and permissions information** is available at www.nature.com/reprints.

**Publisher's note** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

<figure><img src="image_9.png" /></figure>

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International

License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

© The Author(s) 2024
**Additional information**

**Supplementary Information** The online version contains supplementary material available at https://doi.org/10.1038/s41598-024-55922-9.

**Correspondence** and requests for materials should be addressed to C.I.N.S.F.

**Reprints and permissions information** is available at www.nature.com/reprints.

**Publisher's note** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

<figure><img src="image_9.png" /></figure>

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International

License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

© The Author(s) 2024