# Thermodynamic analog of integrate-and-fire neuronal networks by maximum entropy modelling

T S a N Simões, C I N Sampaio Filho, H J Herrmann, J S Andrade, L. de Arcangelis

► To cite this version:

T S a N Simões, C I N Sampaio Filho, H J Herrmann, J S Andrade, L. de Arcangelis. Thermodynamic analog of integrate-and-fire neuronal networks by maximum entropy modelling. Scientific Reports, 2024, 14 (1), pp.9480. 10.1038/s41598-024-60117-3. hal-04583164

HAL Id: hal-04583164

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

Submitted on 22 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
# Thermodynamic analog of integrate-and-fire neuronal networks by maximum entropy modelling

T S a N Simões, C I N Sampaio Filho, H J Herrmann, J S Andrade, L. de Arcangelis

► To cite this version:

T S a N Simões, C I N Sampaio Filho, H J Herrmann, J S Andrade, L. de Arcangelis. Thermodynamic analog of integrate-and-fire neuronal networks by maximum entropy modelling. Scientific Reports, 2024, 14 (1), pp.9480. 10.1038/s41598-024-60117-3. hal-04583164

HAL Id: hal-04583164

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

Submitted on 22 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>

OPEN

# Thermodynamic analog of integrate-and-fire neuronal networks by maximum entropy modelling

T. S. A. N. Simões<fnref n="1" />✉, C. I. N. Sampaio Filho², H. J. Herrmann²,³, J. S. Andrade Jr.² & L. de Arcangelis<fnref n="1" />

Recent results have evidenced that spontaneous brain activity signals are organized in bursts with scale free features and long-range spatio-temporal correlations. These observations have stimulated a theoretical interpretation of results inspired in critical phenomena. In particular, relying on maximum entropy arguments, certain aspects of time-averaged experimental neuronal data have been recently described using Ising-like models, allowing the study of neuronal networks under an analogous thermodynamical framework. This method has been so far applied to a variety of experimental datasets, but never to a biologically inspired neuronal network with short and long-term plasticity. Here, we apply for the first time the Maximum Entropy method to an Integrate-and-fire (IF) model that can be tuned at criticality, offering a controlled setting for a systematic study of criticality and finite-size effects in spontaneous neuronal activity, as opposed to experiments. We consider generalized Ising Hamiltonians whose local magnetic fields and interaction parameters are assigned according to the average activity of single neurons and correlation functions between neurons of the IF networks in the critical state. We show that these Hamiltonians exhibit a spin glass phase for low temperatures, having mostly negative intrinsic fields and a bimodal distribution of interaction constants that tends to become unimodal for larger networks. Results evidence that the magnetization and the response functions exhibit the expected singular behavior near the critical point. Furthermore, we also found that networks with higher percentage of inhibitory neurons lead to Ising-like systems with reduced thermal fluctuations. Finally, considering only neuronal pairs associated with the largest correlation functions allows the study of larger system sizes.

Biological neural networks are highly complex systems, due to the large number of interacting degrees of freedom and the connectivity properties of its constituent elements. Neurons interact by generating action potentials ("firing" or "spiking"), with a biophysical mechanism first explained quantitatively by Hodgkin and Huxley in 1952¹. Their dynamics can be characterized by discretized time series of these firing patterns in binary notation². This description suggests an analogy with the binary spin states in Ising models. Indeed, Ising-like descriptions of brain activity can be traced back to fifty years ago, to the work by Little³ and Hopfield⁴, followed by the work of Amit et al.⁵, where the authors used a spin glass model to describe neural networks. More recently⁶⁻¹³, generalized Ising models have been used to describe the dynamics of recordings from stimulated neuronal activity, using the so-called Maximum Entropy Modelling (MEM) method¹⁴. The method consists in finding the least biased (or maximum entropy) probability distribution that is consistent with a given set of statistical measurements from the system under consideration but otherwise imposes no further constraints¹⁵. When this method is applied to describe the individual firing rates and correlation functions between neurons, the resulting statistical description is equivalent to the one of a specific Ising model with frustrated spins⁶,⁷. Therefore, the MEM approach can be interpreted as an effective mapping of the dynamics from an out-of-equilibrium, multi-component system into a "Hamiltonian" one⁹, specified by relatively few parameters when compared to the full set of possible states of the system in question, allowing to apply the framework of thermodynamics. Outside the scope of neuroscience,

<footnote n="1">¹Department of Mathematics and Physics, University of Campania "Luigi Vanvitelli", Viale Lincoln, 5, 81100 Caserta, Italy. ²Departamento de Física, Fortaleza, Universidade Federal do Ceará, Ceará 60451-970, Brazil. ³ESPCI, PMMH, Paris, 7 quai St., 75005 Bernard, France. ⁴email: tiago.dasilvaalvesdenogueirasimoes@unicampania.it</footnote>
<figure><img src="image_2.png" /></figure>

# Thermodynamic analog of integrate-and-fire neuronal networks by maximum entropy modelling

T. S. A. N. Simões<fnref n="1" />✉, C. I. N. Sampaio Filho², H. J. Herrmann²,³, J. S. Andrade Jr.² & L. de Arcangelis<fnref n="1" />

Recent results have evidenced that spontaneous brain activity signals are organized in bursts with scale free features and long-range spatio-temporal correlations. These observations have stimulated a theoretical interpretation of results inspired in critical phenomena. In particular, relying on maximum entropy arguments, certain aspects of time-averaged experimental neuronal data have been recently described using Ising-like models, allowing the study of neuronal networks under an analogous thermodynamical framework. This method has been so far applied to a variety of experimental datasets, but never to a biologically inspired neuronal network with short and long-term plasticity. Here, we apply for the first time the Maximum Entropy method to an Integrate-and-fire (IF) model that can be tuned at criticality, offering a controlled setting for a systematic study of criticality and finite-size effects in spontaneous neuronal activity, as opposed to experiments. We consider generalized Ising Hamiltonians whose local magnetic fields and interaction parameters are assigned according to the average activity of single neurons and correlation functions between neurons of the IF networks in the critical state. We show that these Hamiltonians exhibit a spin glass phase for low temperatures, having mostly negative intrinsic fields and a bimodal distribution of interaction constants that tends to become unimodal for larger networks. Results evidence that the magnetization and the response functions exhibit the expected singular behavior near the critical point. Furthermore, we also found that networks with higher percentage of inhibitory neurons lead to Ising-like systems with reduced thermal fluctuations. Finally, considering only neuronal pairs associated with the largest correlation functions allows the study of larger system sizes.

Biological neural networks are highly complex systems, due to the large number of interacting degrees of freedom and the connectivity properties of its constituent elements. Neurons interact by generating action potentials ("firing" or "spiking"), with a biophysical mechanism first explained quantitatively by Hodgkin and Huxley in 1952¹. Their dynamics can be characterized by discretized time series of these firing patterns in binary notation². This description suggests an analogy with the binary spin states in Ising models. Indeed, Ising-like descriptions of brain activity can be traced back to fifty years ago, to the work by Little³ and Hopfield⁴, followed by the work of Amit et al.⁵, where the authors used a spin glass model to describe neural networks. More recently⁶⁻¹³, generalized Ising models have been used to describe the dynamics of recordings from stimulated neuronal activity, using the so-called Maximum Entropy Modelling (MEM) method¹⁴. The method consists in finding the least biased (or maximum entropy) probability distribution that is consistent with a given set of statistical measurements from the system under consideration but otherwise imposes no further constraints¹⁵. When this method is applied to describe the individual firing rates and correlation functions between neurons, the resulting statistical description is equivalent to the one of a specific Ising model with frustrated spins⁶,⁷. Therefore, the MEM approach can be interpreted as an effective mapping of the dynamics from an out-of-equilibrium, multi-component system into a "Hamiltonian" one⁹, specified by relatively few parameters when compared to the full set of possible states of the system in question, allowing to apply the framework of thermodynamics. Outside the scope of neuroscience,

<footnote n="1">¹Department of Mathematics and Physics, University of Campania "Luigi Vanvitelli", Viale Lincoln, 5, 81100 Caserta, Italy. ²Departamento de Física, Fortaleza, Universidade Federal do Ceará, Ceará 60451-970, Brazil. ³ESPCI, PMMH, Paris, 7 quai St., 75005 Bernard, France. ⁴email: tiago.dasilvaalvesdenogueirasimoes@unicampania.it</footnote>
this method has been used to construct novel statistical models¹⁶, and to study geographic distributions of species¹⁷, protein structures¹⁸,¹⁹, gene mutation effects²⁰, the collective behavior of flocks of birds²¹ and their diversity distribution²², correlations in eye movements while watching videos²³ or reading texts²⁴, and to conduct urban-oriented studies such as flood risk assessment²⁵, analyzing urban mobility patterns²⁶ and property valuations in real estate markets²⁷. In turn, within the context of neuroscience, MEM has been mostly used to analyse experimental data of neuronal recordings, such as visual inputs from cells of a salamander⁶,⁹,¹⁰,²⁸,²⁹ and rat retina³⁰, numerical simulations from a phenomenological model of retinal ganglion cells³¹, responses from hippocampus place cells in rodent brains¹¹,¹³, synchronized and desynchronized neuronal activity from the primary visual cortex of anaesthetised cats and awake monkeys³², in vivo and in vitro neuronal activity from cortical tissue of rodents³³ and the nervous system of the nematod *C. Elegans*¹², an organism whose pattern of connectivity between all its 302 neurons is well known³⁴.

Overall, these studies have given insights regarding collective⁶ and functional characteristics²⁸ of biological neural networks, as well as helping unraveling the information content of neuronal responses⁹,¹⁰ and its analogous thermodynamic properties¹⁰. Particularly, it was shown that the thermal fluctuations present in these Ising systems tend to diverge with system size,⁷,¹⁰,¹²,³²,³³ a sign of critical behaviour. However, experimental measurements in real networks can pose challenges that may impact the relevance and range of applicability of the conclusions drawn using the MEM method. For instance, contemporary neuronal recording techniques are restricted by the duration and sampling rate of the recordings²⁹,³⁵ and a precise association between which neuron generated which spike is an open problem, known as “spike sorting”³⁶, possibly affecting the accuracy of the measured correlation functions³⁷. Furthermore, precise estimations of the fraction of excitatory and inhibitory neuronal populations in real systems are also often difficult to obtain³⁸. To mitigate these challenges, numerical models can offer a more controlled environment for studying neuronal activity, allowing one to systematically change many parameters of choice, tuning the activity state, make systematic size studies which do not rely on subsampling, and generate many equivalent independent samples, thus producing statistically relevant data. However, this method has never been applied to data generated by biologically inspired neuronal network models reproducing the fundamental features of neuronal activity in the resting state.

In this context, we apply the MEM method to spontaneous neuronal activity generated by an Integrate-and-fire (IF) model implemented on a scale-free network with short- and long-term plasticity³⁹. This IF model implements the main features of neuronal activity, as firing at threshold, refractory period, as well as long- and short-term plasticity for synaptic connections. In particular, short-term plasticity models the recovery of synaptic resources and relies on a tuning parameter that regulates the dynamical state of the system³⁹. At appropriate values of the tuning parameter, the model generates bursts of firing activity with characteristic spatio-temporal statistics, known as neuronal avalanches, first observed in 2003, in acute slices of rat cortex⁴⁰, displaying sizes S and durations *D* that are power-law distributed, as *P(S)* ~ *S*⁻¹.⁵ and *P(D)* ~ *D*⁻², suggesting that biological neural networks operate near a critical point, as proposed in previous numerical and analytical studies⁴¹,⁴², and is consistent with the divergence of thermal energy fluctuations observed in relatively small Ising-like models constructed from experimental neuronal data using the MEM method⁷,¹⁰,¹²,³². Results suggest that the brain might act between a quiescence-like state and a state of hyperactivity, possibly offering several biological advantages such as optimal information transmission and storage⁴³. Signs of criticality have been observed in a wide variety of animals, including humans⁴⁴⁻⁴⁶, monkeys⁴⁷, cats³², salamanders⁴⁸, turtles⁴⁹, worms¹², and fish⁵⁰, a hallmark of universality⁴¹.

Our aim here is to apply the MEM approach to neuronal activity data generated by the IF neuronal network in the critical state. In particular, we construct, using the MEM method mentioned above, fully-connected Ising models of frustrated spins with local fields and interaction constants which allow the average spin state to reproduce the average local neuronal activity, as well as the two-point correlation functions in the two models. We will investigate how the parameters of the Ising model and its thermodynamics properties change with the size and connectivity of the neuronal network. Additionally, we examine how different fractions of inhibitory neurons may influence the associated Ising models. We also consider partially-connected spin networks with only a subset of the total pairwise interactions, associated with the strongest neuronal pairwise correlations, to be able to analyse networks of sizes larger than the ones typically considered so far in the literature³⁵. The aim is to implement the mapping of the neuronal model into a thermodynamic framework, which may open the way to novel insights into brain functions.

The manuscript is organized in four sections. The first one gives a general presentation of Results, in particular the firing dynamics generated by the IF model, the mapping of the IF model into an Ising-like model using the MEM method and a study of its associated properties. In the following sections, we present a general Discussion followed by the Conclusions. In the final section, we describe the details of all methods, namely the IF model implementation and the MEM method.

## Results

We consider scale-free neural networks of different size $N \in [20, 500]$ with short- and long-term plasticity, as well as a refractory mechanism, where neurons remain inactive for a single timestep immediately after firing (see “Methods”). Together with plastic adaptation, the refractory mechanism has been shown to play a crucial role in the appearance of critical behaviour and in shaping the topology of the network⁵¹. All measurements are performed on networks in the critical state. This is done by appropriately tuning the neurotransmitter recovery parameter $\delta u_{\text{rec}}$ (see “Integrate-and-fire model” in Methods) and verifying that the distributions of neuronal avalanche sizes $P(S) \propto S^{-\tau_S}$ and durations $P(D) \propto D^{-\tau_D}$ decay as power laws with exponents $\tau_S = 1.5$ and $\tau_D = 2$ (see Fig. S1 in Supplementary Information), characteristic of the mean field self-organized branching process⁵², and consistent with both numerical and experimental observations³⁹,⁴⁰,⁴⁶. Initially, we consider only
this method has been used to construct novel statistical models¹⁶, and to study geographic distributions of species¹⁷, protein structures¹⁸,¹⁹, gene mutation effects²⁰, the collective behavior of flocks of birds²¹ and their diversity distribution²², correlations in eye movements while watching videos²³ or reading texts²⁴, and to conduct urban-oriented studies such as flood risk assessment²⁵, analyzing urban mobility patterns²⁶ and property valuations in real estate markets²⁷. In turn, within the context of neuroscience, MEM has been mostly used to analyse experimental data of neuronal recordings, such as visual inputs from cells of a salamander⁶,⁹,¹⁰,²⁸,²⁹ and rat retina³⁰, numerical simulations from a phenomenological model of retinal ganglion cells³¹, responses from hippocampus place cells in rodent brains¹¹,¹³, synchronized and desynchronized neuronal activity from the primary visual cortex of anaesthetised cats and awake monkeys³², in vivo and in vitro neuronal activity from cortical tissue of rodents³³ and the nervous system of the nematod *C. Elegans*¹², an organism whose pattern of connectivity between all its 302 neurons is well known³⁴.

Overall, these studies have given insights regarding collective⁶ and functional characteristics²⁸ of biological neural networks, as well as helping unraveling the information content of neuronal responses⁹,¹⁰ and its analogous thermodynamic properties¹⁰. Particularly, it was shown that the thermal fluctuations present in these Ising systems tend to diverge with system size,⁷,¹⁰,¹²,³²,³³ a sign of critical behaviour. However, experimental measurements in real networks can pose challenges that may impact the relevance and range of applicability of the conclusions drawn using the MEM method. For instance, contemporary neuronal recording techniques are restricted by the duration and sampling rate of the recordings²⁹,³⁵ and a precise association between which neuron generated which spike is an open problem, known as “spike sorting”³⁶, possibly affecting the accuracy of the measured correlation functions³⁷. Furthermore, precise estimations of the fraction of excitatory and inhibitory neuronal populations in real systems are also often difficult to obtain³⁸. To mitigate these challenges, numerical models can offer a more controlled environment for studying neuronal activity, allowing one to systematically change many parameters of choice, tuning the activity state, make systematic size studies which do not rely on subsampling, and generate many equivalent independent samples, thus producing statistically relevant data. However, this method has never been applied to data generated by biologically inspired neuronal network models reproducing the fundamental features of neuronal activity in the resting state.

In this context, we apply the MEM method to spontaneous neuronal activity generated by an Integrate-and-fire (IF) model implemented on a scale-free network with short- and long-term plasticity³⁹. This IF model implements the main features of neuronal activity, as firing at threshold, refractory period, as well as long- and short-term plasticity for synaptic connections. In particular, short-term plasticity models the recovery of synaptic resources and relies on a tuning parameter that regulates the dynamical state of the system³⁹. At appropriate values of the tuning parameter, the model generates bursts of firing activity with characteristic spatio-temporal statistics, known as neuronal avalanches, first observed in 2003, in acute slices of rat cortex⁴⁰, displaying sizes S and durations *D* that are power-law distributed, as *P(S)* ~ *S*⁻¹.⁵ and *P(D)* ~ *D*⁻², suggesting that biological neural networks operate near a critical point, as proposed in previous numerical and analytical studies⁴¹,⁴², and is consistent with the divergence of thermal energy fluctuations observed in relatively small Ising-like models constructed from experimental neuronal data using the MEM method⁷,¹⁰,¹²,³². Results suggest that the brain might act between a quiescence-like state and a state of hyperactivity, possibly offering several biological advantages such as optimal information transmission and storage⁴³. Signs of criticality have been observed in a wide variety of animals, including humans⁴⁴⁻⁴⁶, monkeys⁴⁷, cats³², salamanders⁴⁸, turtles⁴⁹, worms¹², and fish⁵⁰, a hallmark of universality⁴¹.

Our aim here is to apply the MEM approach to neuronal activity data generated by the IF neuronal network in the critical state. In particular, we construct, using the MEM method mentioned above, fully-connected Ising models of frustrated spins with local fields and interaction constants which allow the average spin state to reproduce the average local neuronal activity, as well as the two-point correlation functions in the two models. We will investigate how the parameters of the Ising model and its thermodynamics properties change with the size and connectivity of the neuronal network. Additionally, we examine how different fractions of inhibitory neurons may influence the associated Ising models. We also consider partially-connected spin networks with only a subset of the total pairwise interactions, associated with the strongest neuronal pairwise correlations, to be able to analyse networks of sizes larger than the ones typically considered so far in the literature³⁵. The aim is to implement the mapping of the neuronal model into a thermodynamic framework, which may open the way to novel insights into brain functions.

The manuscript is organized in four sections. The first one gives a general presentation of Results, in particular the firing dynamics generated by the IF model, the mapping of the IF model into an Ising-like model using the MEM method and a study of its associated properties. In the following sections, we present a general Discussion followed by the Conclusions. In the final section, we describe the details of all methods, namely the IF model implementation and the MEM method.

## Results

We consider scale-free neural networks of different size $N \in [20, 500]$ with short- and long-term plasticity, as well as a refractory mechanism, where neurons remain inactive for a single timestep immediately after firing (see “Methods”). Together with plastic adaptation, the refractory mechanism has been shown to play a crucial role in the appearance of critical behaviour and in shaping the topology of the network⁵¹. All measurements are performed on networks in the critical state. This is done by appropriately tuning the neurotransmitter recovery parameter $\delta u_{\text{rec}}$ (see “Integrate-and-fire model” in Methods) and verifying that the distributions of neuronal avalanche sizes $P(S) \propto S^{-\tau_S}$ and durations $P(D) \propto D^{-\tau_D}$ decay as power laws with exponents $\tau_S = 1.5$ and $\tau_D = 2$ (see Fig. S1 in Supplementary Information), characteristic of the mean field self-organized branching process⁵², and consistent with both numerical and experimental observations³⁹,⁴⁰,⁴⁶. Initially, we consider only
fully-excitatory networks. Later, we will also consider networks with a fraction $p_{in} > 0$ of inhibitory neurons. The time evolution of a system with $N = 120$ neurons in the critical state is shown in the raster plot of Fig. 1.

First, we analyse the firing statistics of the IF model networks. Then, we apply the MEM method to map the dynamics of the IF model networks of size $N \in [20, 120]$ into a pairwise Ising model, defined by $N$ “local fields” $h_i$ and $N \cdot (N - 1)/2$ “interaction constants” $J_{ij}$, and study its thermodynamical properties. More specifically, this mapping is achieved by using the so-called Boltzmann Machine (BM) algorithm. This learning algorithm searches for the set of parameters $h_i$ and $J_{ij}$ that best fit the pairwise Ising model to the data of the IF model by comparison with the average local activities $\langle \sigma_i \rangle$ and correlation functions $C_{ij}$ generated by each model (see “Maximum entropy modelling” in Methods).

## Firing statistics

To measure the IF network firing statistics, we consider a time bin $\Delta t_b = 5$ timesteps. We assign a binary value $\sigma_i^k \in \{-1, 1\}$ to each neuron at each time bin according to the following rule: $\sigma_i^k = 1$ if neuron $i$ fired at least once during the $k$-th time bin, $\sigma_i^k = -1$ if the neuron $i$ was inactive during the whole time bin $k$. We then calculate the average local activity $\langle \sigma_i \rangle$ of the $N$ neurons over $N_b$ time bins, defined as

$$ \langle \sigma_i \rangle = \frac{1}{N_b} \sum_{k=1}^{N_b} \sigma_i^k, \qquad (1) $$

as well as the average two-point activity $\langle \sigma_i \sigma_j \rangle$ between all $N \cdot (N - 1)/2$ pairs of neurons $i$ and $j$

$$ \langle \sigma_i \sigma_j \rangle = \frac{1}{N_b} \sum_{k=1}^{N_b} \sigma_i^k \sigma_j^k. \qquad (2) $$

From these quantities we can also define the two-point correlation functions $C_{ij}$,

$$ C_{ij} = \langle (\sigma_i - \langle \sigma_i \rangle) \cdot (\sigma_j - \langle \sigma_j \rangle) \rangle = \langle \sigma_i \sigma_j \rangle - \langle \sigma_i \rangle \langle \sigma_j \rangle, \qquad (3) $$

<figure><img src="image_3.png" /><figcaption>**Figure 1.** Time series for a system at criticality with $N = 120$ neurons. The top raster plot presents the time evolution of the firing states of $N = 120$ neurons over a certain number of timesteps, from a fully-excitatory IF network in the critical state. A coloured dot indicates that the respective neuron fired during that timestep, while an absence of a dot means it was inactive. The alternating green and blue colours are just a visual aid to distinguish between different avalanches. The bottom raster plot is a zoom-in of the top one, where the coloured shaded areas indicate avalanches of size $S > 1$ and duration $D$.</figcaption></figure>
fully-excitatory networks. Later, we will also consider networks with a fraction $p_{in} > 0$ of inhibitory neurons. The time evolution of a system with $N = 120$ neurons in the critical state is shown in the raster plot of Fig. 1.

First, we analyse the firing statistics of the IF model networks. Then, we apply the MEM method to map the dynamics of the IF model networks of size $N \in [20, 120]$ into a pairwise Ising model, defined by $N$ “local fields” $h_i$ and $N \cdot (N - 1)/2$ “interaction constants” $J_{ij}$, and study its thermodynamical properties. More specifically, this mapping is achieved by using the so-called Boltzmann Machine (BM) algorithm. This learning algorithm searches for the set of parameters $h_i$ and $J_{ij}$ that best fit the pairwise Ising model to the data of the IF model by comparison with the average local activities $\langle \sigma_i \rangle$ and correlation functions $C_{ij}$ generated by each model (see “Maximum entropy modelling” in Methods).

## Firing statistics

To measure the IF network firing statistics, we consider a time bin $\Delta t_b = 5$ timesteps. We assign a binary value $\sigma_i^k \in \{-1, 1\}$ to each neuron at each time bin according to the following rule: $\sigma_i^k = 1$ if neuron $i$ fired at least once during the $k$-th time bin, $\sigma_i^k = -1$ if the neuron $i$ was inactive during the whole time bin $k$. We then calculate the average local activity $\langle \sigma_i \rangle$ of the $N$ neurons over $N_b$ time bins, defined as

$$ \langle \sigma_i \rangle = \frac{1}{N_b} \sum_{k=1}^{N_b} \sigma_i^k, \qquad (1) $$

as well as the average two-point activity $\langle \sigma_i \sigma_j \rangle$ between all $N \cdot (N - 1)/2$ pairs of neurons $i$ and $j$

$$ \langle \sigma_i \sigma_j \rangle = \frac{1}{N_b} \sum_{k=1}^{N_b} \sigma_i^k \sigma_j^k. \qquad (2) $$

From these quantities we can also define the two-point correlation functions $C_{ij}$,

$$ C_{ij} = \langle (\sigma_i - \langle \sigma_i \rangle) \cdot (\sigma_j - \langle \sigma_j \rangle) \rangle = \langle \sigma_i \sigma_j \rangle - \langle \sigma_i \rangle \langle \sigma_j \rangle, \qquad (3) $$

<figure><img src="image_3.png" /><figcaption>**Figure 1.** Time series for a system at criticality with $N = 120$ neurons. The top raster plot presents the time evolution of the firing states of $N = 120$ neurons over a certain number of timesteps, from a fully-excitatory IF network in the critical state. A coloured dot indicates that the respective neuron fired during that timestep, while an absence of a dot means it was inactive. The alternating green and blue colours are just a visual aid to distinguish between different avalanches. The bottom raster plot is a zoom-in of the top one, where the coloured shaded areas indicate avalanches of size $S > 1$ and duration $D$.</figcaption></figure>
where $\langle\sigma_i\rangle$ is closely related to the firing rate $r_i = (\langle\sigma_i\rangle + 1)/2\Delta t_b$ of neuron $i^9$, whereas $C_{ij}$ quantifies the tendency for neurons $i$ and $j$ to fire together in the time interval $\Delta t_b$. Notice that, although likely related, the $C_{ij}$ are completely distinct from the synaptic strengths $w_{ij}$. For instance, the $C_{ij}$ are defined for each pair $i$ and $j$ and are symmetric by definition ($C_{ij} = C_{ji}$), whereas $w_{ij}$ are not. For all simulations we use $N_b = 10^7$.

In Fig. 2 we show the distributions of $\langle\sigma_i\rangle$ and $C_{ij}$ obtained from $N_c$ different neural network configurations (see Fig. S1 in Supplementary Information for the respective avalanche statistics), for systems with different number of neurons $N$. We vary $N_c$ according to $N$ so that the total number of neurons considered $NN_c \sim 10,000$. The average local activity $\langle\sigma_i\rangle$ is negative for all neurons (Fig. 2a-c), implying that firing is a relatively rare event when considering individual neurons. Moreover, the distribution becomes narrower and shifts towards more negative values as the system size increases. In Fig. 2d-f, we show that the correlation functions are small, but mostly non-zero for all neuron pairs, with their distributions peaked near zero, becoming sharper as the system size $N$ increases. It is also noteworthy that, in experimental studies of the vertebrate retina$^6$, even weak pairwise correlations have been shown to likely become statistically significant in large networks, influencing the dynamics at the scale of the whole network.

Another useful quantity that characterizes the firing dynamics is the probability $P(K)$ to observe $K \in [0, N]$ neurons firing simultaneously within a time bin $\Delta t_b$,

$$P(K) = \frac{1}{N_b} \sum_{k=1}^{N_b} \delta_{K,K^k}, \qquad (4)$$

where $\delta_{K,K^k}$ is the Kronecker delta function and $K^k = \sum_i^N (\sigma_i^k + 1)/2$ counts the number of neurons firing during the $k$-th time bin. In Fig. 3 we plot the probability (4) for several networks of different size $N$. For all $N$, $P(K)$ is well fitted by an exponential distribution in an intermediate regime. The exponential factor of the distributions consistently decreases as $N$ increases, indicating that it is more likely to observe concurrent firing between neurons in larger networks.

## Results of the Boltzmann machine learning

Next, we consider Ising models with pairwise interactions for different system sizes $N = \{20, 40, 80, 120\}$. The interaction constants and local magnetic fields are fitted to reproduce the temporal averages of the correlation functions and local activities generated by the IF model, using the Boltzmann Machine (BM) algorithm (see “Maximum entropy modelling” in Methods). Each spin $i$ can either be in the up-state ($\sigma_i = +1$) or down-state ($\sigma_i = -1$), defining a particular spin configuration $\sigma = \{\sigma_1, \sigma_2, ..., \sigma_N\}$. The Ising Hamiltonian $H(\sigma)$ is given by

<figure><img src="image_4.png" /><figcaption>**Figure 2.** Distributions of the average local activity $\langle\sigma_i\rangle$ and of the correlation functions $C_{ij}$ in IF networks. Average local activity $\langle\sigma_i\rangle$ (a-c) and correlation functions $C_{ij}$ (d-f), calculated as averages over $N_b = 10^7$ time bins. Distributions are obtained for sizes $N = \{40, 120, 500\}$, for $N_c \sim 10,000/N$ different fully-excitatory networks.</figcaption></figure>
where $\langle \sigma_i \rangle$ is closely related to the firing rate $r_i = (\langle \sigma_i \rangle + 1)/2\Delta t_b$ of neuron $i^9$, whereas $C_{ij}$ quantifies the tendency for neurons $i$ and $j$ to fire together in the time interval $\Delta t_b$. Notice that, although likely related, the $C_{ij}$ are completely distinct from the synaptic strengths $w_{ij}$. For instance, the $C_{ij}$ are defined for each pair $i$ and $j$ and are symmetric by definition ($C_{ij} = C_{ji}$), whereas $w_{ij}$ are not. For all simulations we use $N_b = 10^7$.

In Fig. 2 we show the distributions of $\langle \sigma_i \rangle$ and $C_{ij}$ obtained from $N_c$ different neural network configurations (see Fig. S1 in Supplementary Information for the respective avalanche statistics), for systems with different number of neurons $N$. We vary $N_c$ according to $N$ so that the total number of neurons considered $NN_c \sim 10,000$. The average local activity $\langle \sigma_i \rangle$ is negative for all neurons (Fig. 2a-c), implying that firing is a relatively rare event when considering individual neurons. Moreover, the distribution becomes narrower and shifts towards more negative values as the system size increases. In Fig. 2d-f, we show that the correlation functions are small, but mostly non-zero for all neuron pairs, with their distributions peaked near zero, becoming sharper as the system size $N$ increases. It is also noteworthy that, in experimental studies of the vertebrate retina$^6$, even weak pairwise correlations have been shown to likely become statistically significant in large networks, influencing the dynamics at the scale of the whole network.

Another useful quantity that characterizes the firing dynamics is the probability $P(K)$ to observe $K \in [0, N]$ neurons firing simultaneously within a time bin $\Delta t_b$,

$$P(K) = \frac{1}{N_b} \sum_{k=1}^{N_b} \delta_{K,K^k}, \qquad (4)$$

where $\delta_{K,K^k}$ is the Kronecker delta function and $K^k = \sum_i^N (\sigma_i^k + 1)/2$ counts the number of neurons firing during the $k$-th time bin. In Fig. 3 we plot the probability (4) for several networks of different size $N$. For all $N$, $P(K)$ is well fitted by an exponential distribution in an intermediate regime. The exponential factor of the distributions consistently decreases as $N$ increases, indicating that it is more likely to observe concurrent firing between neurons in larger networks.

## Results of the Boltzmann machine learning

Next, we consider Ising models with pairwise interactions for different system sizes $N = \{20, 40, 80, 120\}$. The interaction constants and local magnetic fields are fitted to reproduce the temporal averages of the correlation functions and local activities generated by the IF model, using the Boltzmann Machine (BM) algorithm (see “Maximum entropy modelling” in Methods). Each spin $i$ can either be in the up-state ($\sigma_i = +1$) or down-state ($\sigma_i = -1$), defining a particular spin configuration $\sigma = \{\sigma_1, \sigma_2, ..., \sigma_N\}$. The Ising Hamiltonian $H(\sigma)$ is given by

<figure><img src="image_4.png" /><figcaption>**Figure 2.** Distributions of the average local activity $\langle \sigma_i \rangle$ and of the correlation functions $C_{ij}$ in IF networks. Average local activity $\langle \sigma_i \rangle$ (a-c) and correlation functions $C_{ij}$ (d-f), calculated as averages over $N_b = 10^7$ time bins. Distributions are obtained for sizes $N = \{40, 120, 500\}$, for $N_c \sim 10,000/N$ different fully-excitatory networks.</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Probability $P(K)$ of observing $K$ neurons firing simultaneously during a time bin $\Delta t_b = 5$ timesteps. Results are averaged over $N_c$ different fully-excitatory IF networks configurations of size $N$. $N_c$ varies with $N$ as in Fig. 2. For each configuration, $P(K)$ is estimated by averaging over $N_b = 10^7$ bins. Error bars are given by the standard error, and most of them are smaller than the symbol size.</figcaption></figure>

$$H(\sigma) = -\sum_{i}^{N} h_i \sigma_i - \sum_{i}^{N} \sum_{j<i}^{N} J_{ij} \sigma_i \sigma_j, \qquad (5)$$

with $h_i$ and $J_{ij}$ the fitting parameters, acting as local fields on spin $i$ and interaction constants between spins $i$ and $j$, respectively. When appropriate, we use the superscript (IF) and (BM) to distinguish between measurements using the IF model and those obtained with the BM using Monte Carlo. More specifically, $\langle \dots \rangle^{(\text{IF})}$ indicates an average over time bins in the IF model while $\langle \dots \rangle^{(\text{BM})}$ is an average over spin configurations of the Ising model.

In Fig. 4 we present the values of {$h_i$} and the distributions of the interaction constants {$J_{ij}$}, fitted to reproduce the temporal averages {$\langle \sigma_i \rangle^{(\text{IF})}$} and {$\langle \sigma_i \sigma_j \rangle^{(\text{IF})}$} of fully-excitatory IF networks of size $N$, tuned to the critical state.

More specifically, Fig. 4a–d show ranked plots of the fields $h_i$, where $h_1$ ($h_N$) corresponds to the spin associated with the neuron that fired the most (least) in the IF model. The fields $h_i$ are mostly negative for all $N$. This is

<figure><img src="image_6.png" /><figcaption>**Figure 4.** Sets of learned fields $h_i$ and interaction constants $J_{ij}$ of the Ising model reproducing the time averages of the IF model. Plots of the fields $h_i$ (a–d), sorted by the average local activity $\langle \sigma_i \rangle^{(\text{IF})}$ of the associated neuron $i$, in order of decreasing $\langle \sigma_i \rangle^{(\text{IF})}$, and distributions of the interaction constants $J_{ij}$ (e–h), learned by the BM to reproduce the average local activity $\langle \sigma_i \rangle$ and two-point correlation functions $C_{ij}$ of a fully-excitatory IF neural network at criticality, with $N = \{20, 40, 80, 120\}$ neurons.</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Probability $P(K)$ of observing $K$ neurons firing simultaneously during a time bin $\Delta t_b = 5$ timesteps. Results are averaged over $N_c$ different fully-excitatory IF networks configurations of size $N$. $N_c$ varies with $N$ as in Fig. 2. For each configuration, $P(K)$ is estimated by averaging over $N_b = 10^7$ bins. Error bars are given by the standard error, and most of them are smaller than the symbol size.</figcaption></figure>

$$H(\sigma) = -\sum_{i}^{N} h_i \sigma_i - \sum_{i}^{N} \sum_{j<i}^{N} J_{ij} \sigma_i \sigma_j, \qquad (5)$$

with $h_i$ and $J_{ij}$ the fitting parameters, acting as local fields on spin $i$ and interaction constants between spins $i$ and $j$, respectively. When appropriate, we use the superscript (IF) and (BM) to distinguish between measurements using the IF model and those obtained with the BM using Monte Carlo. More specifically, $\langle \dots \rangle^{(IF)}$ indicates an average over time bins in the IF model while $\langle \dots \rangle^{(BM)}$ is an average over spin configurations of the Ising model.

In Fig. 4 we present the values of {$h_i$} and the distributions of the interaction constants {$J_{ij}$}, fitted to reproduce the temporal averages {$\langle \sigma_i \rangle^{(IF)}$} and {$\langle \sigma_i \sigma_j \rangle^{(IF)}$} of fully-excitatory IF networks of size $N$, tuned to the critical state.

More specifically, Fig. 4a–d show ranked plots of the fields $h_i$, where $h_1$ ($h_N$) corresponds to the spin associated with the neuron that fired the most (least) in the IF model. The fields $h_i$ are mostly negative for all $N$. This is

<figure><img src="image_6.png" /><figcaption>**Figure 4.** Sets of learned fields $h_i$ and interaction constants $J_{ij}$ of the Ising model reproducing the time averages of the IF model. Plots of the fields $h_i$ (a–d), sorted by the average local activity $\langle \sigma_i \rangle^{(IF)}$ of the associated neuron $i$, in order of decreasing $\langle \sigma_i \rangle^{(IF)}$, and distributions of the interaction constants $J_{ij}$ (e–h), learned by the BM to reproduce the average local activity $\langle \sigma_i \rangle$ and two-point correlation functions $C_{ij}$ of a fully-excitatory IF neural network at criticality, with $N = \{20, 40, 80, 120\}$ neurons.</figcaption></figure>
not surprising considering that a neuron in the IF model fires rarely, as discussed in previous sections, and the negative fields are a consequence of this sparse activity. We note, however, that, as the system size N increases, an increasingly larger fraction of spins have positive $h_i$, corresponding to the spins associated with the most active neurons. In Fig. 4e–h we plot the distributions of the learned interaction constants $J_{ij}$. The distributions are bimodal, with the absolute maximum near zero and the second one at positive $J_{ij}$. However, the height of the second peak appears to decrease with the system size N, suggesting that the limiting distribution when $N \to \infty$ becomes consistent with a normal distribution centered around zero. It is interesting to recall that random, normally distributed $J_{ij}$ characterize the so-called Sherrington-Kirkpatrick (SK) model$^{53,54}$, an Ising-like model which has been shown to exhibit both ferromagnetic and spin glass phases$^{53}$. Unlike the SK model, however, the learned Ising models exhibit high heterogeneity in the distribution of the fields $h_i$, and, since the learning technique fits simultaneously both the fields and interaction constants, the inference of the {$J_{ij}$} should not be considered numerically decoupled from that of the {$h_i$}.

To test the quality of the learning process, in Fig. 5 we compare the {$\langle\sigma_i\rangle^{(BM)}$} and {$C_{ij}^{(BM)}$} of the Ising-like model with the neural network data of the IF model {$\langle\sigma_i\rangle^{(IF)}$} and {$C_{ij}^{(IF)}$}.

As described above (by virtue of Eqs. (16) and (17), in the “Methods” section), if the learning is successful, these quantities should be identical since these are precisely the constraints imposed in the learning procedure. The validity of this method is confirmed for both $\langle\sigma_i\rangle$ and $C_{ij}$, for all N, with all points following closely the bisector $y = x$.

An interesting quantity to assess the predictive power of the model are the three-point correlation functions $T_{ijk}$ between all $N \cdot (N-1) \cdot (N-2)/6$ triplets of neurons i, j and k,

$$ 
\begin{aligned}
T_{ijk} &= \langle (\sigma_i - \langle\sigma_i\rangle) \cdot (\sigma_j - \langle\sigma_j\rangle) \cdot (\sigma_k - \langle\sigma_k\rangle) \rangle \\
&= \langle \sigma_i \sigma_j \sigma_k \rangle - \langle \sigma_i \rangle \langle \sigma_j \sigma_k \rangle - \langle \sigma_j \rangle \langle \sigma_i \sigma_k \rangle - \langle \sigma_k \rangle \langle \sigma_i \sigma_j \rangle + 2 \langle \sigma_i \rangle \langle \sigma_j \rangle \langle \sigma_k \rangle,
\end{aligned} 
\quad (6) $$

where the averages $\langle...\rangle$ are defined analogously as in Eqs. (1) or (3). In Fig. 6 we compare the triplets between the IF and the Ising models, for systems with $N = \{20, 40, 80, 120\}$ considered previously.

Even though the $T_{ijk}$ are not constrained by the learning procedure, the generalized Ising model still captures their systematic qualitative behaviour. Unlike {$\langle\sigma_i\rangle$} and {$C_{ij}$}, however, the triplets appear to be consistently over-estimated by the Ising model, when compared to their IF model counterpart. We note that this behaviour was also seen in experimental recordings from stimulated neurons in a salamander retina$^7$, albeit to a less significant degree.

Another useful measure to assess how well the Ising-like model describes the IF model data is the probability for simultaneous firing $P(K)$ within a small time window, as defined in Eq. (4). In the Ising model, the sum over time bins in Eq. (4) is replaced by a sum over spin configurations generated by Monte Carlo simulation, where

<figure><img src="image_7.png" /><figcaption>**Figure 5.** Quality test for the BM learning process. Comparison between the average local activity $\langle\sigma_i\rangle$ (a–d) and correlation functions $C_{ij}$ (e–h) of the Ising-like model (y-axes) and IF model network (x-axes), for system sizes $N = \{20, 40, 80, 120\}$. The blue dashed lines are the bisector $y = x$. Results are averages over $N_b = 10^7$ time bins for the IF model (IF), and over $M_c = 3 \cdot 10^6$ spin configurations for the Ising model (BM). Error bars are given by the standard error, and are smaller or equal to the symbol size.</figcaption></figure>
not surprising considering that a neuron in the IF model fires rarely, as discussed in previous sections, and the negative fields are a consequence of this sparse activity. We note, however, that, as the system size N increases, an increasingly larger fraction of spins have positive $h_i$, corresponding to the spins associated with the most active neurons. In Fig. 4e–h we plot the distributions of the learned interaction constants $J_{ij}$. The distributions are bimodal, with the absolute maximum near zero and the second one at positive $J_{ij}$. However, the height of the second peak appears to decrease with the system size N, suggesting that the limiting distribution when $N \to \infty$ becomes consistent with a normal distribution centered around zero. It is interesting to recall that random, normally distributed $J_{ij}$ characterize the so-called Sherrington-Kirkpatrick (SK) model$^{53,54}$, an Ising-like model which has been shown to exhibit both ferromagnetic and spin glass phases$^{53}$. Unlike the SK model, however, the learned Ising models exhibit high heterogeneity in the distribution of the fields $h_i$, and, since the learning technique fits simultaneously both the fields and interaction constants, the inference of the {$J_{ij}$} should not be considered numerically decoupled from that of the {$h_i$}.

To test the quality of the learning process, in Fig. 5 we compare the {$\langle\sigma_i\rangle^{(BM)}$} and {$C_{ij}^{(BM)}$} of the Ising-like model with the neural network data of the IF model {$\langle\sigma_i\rangle^{(IF)}$} and {$C_{ij}^{(IF)}$}.

As described above (by virtue of Eqs. (16) and (17), in the “Methods” section), if the learning is successful, these quantities should be identical since these are precisely the constraints imposed in the learning procedure. The validity of this method is confirmed for both $\langle\sigma_i\rangle$ and $C_{ij}$, for all N, with all points following closely the bisector $y = x$.

An interesting quantity to assess the predictive power of the model are the three-point correlation functions $T_{ijk}$ between all $N \cdot (N-1) \cdot (N-2)/6$ triplets of neurons i, j and k,

$$ T_{ijk} = \langle (\sigma_i - \langle\sigma_i\rangle) \cdot (\sigma_j - \langle\sigma_j\rangle) \cdot (\sigma_k - \langle\sigma_k\rangle) \rangle \\ = \langle \sigma_i \sigma_j \sigma_k \rangle - \langle \sigma_i \rangle \langle \sigma_j \sigma_k \rangle - \langle \sigma_j \rangle \langle \sigma_i \sigma_k \rangle - \langle \sigma_k \rangle \langle \sigma_i \sigma_j \rangle + 2 \langle \sigma_i \rangle \langle \sigma_j \rangle \langle \sigma_k \rangle, \quad (6) $$

where the averages $\langle...\rangle$ are defined analogously as in Eqs. (1) or (3). In Fig. 6 we compare the triplets between the IF and the Ising models, for systems with $N = \{20, 40, 80, 120\}$ considered previously.

Even though the $T_{ijk}$ are not constrained by the learning procedure, the generalized Ising model still captures their systematic qualitative behaviour. Unlike $\langle\sigma_i\rangle$ and $C_{ij}$, however, the triplets appear to be consistently over-estimated by the Ising model, when compared to their IF model counterpart. We note that this behaviour was also seen in experimental recordings from stimulated neurons in a salamander retina$^7$, albeit to a less significant degree.

Another useful measure to assess how well the Ising-like model describes the IF model data is the probability for simultaneous firing $P(K)$ within a small time window, as defined in Eq. (4). In the Ising model, the sum over time bins in Eq. (4) is replaced by a sum over spin configurations generated by Monte Carlo simulation, where

<figure><img src="image_7.png" /><figcaption>**Figure 5.** Quality test for the BM learning process. Comparison between the average local activity $\langle\sigma_i\rangle$ (a–d) and correlation functions $C_{ij}$ (e–h) of the Ising-like model (y-axes) and IF model network (x-axes), for system sizes $N = \{20, 40, 80, 120\}$. The blue dashed lines are the bisector $y = x$. Results are averages over $N_b = 10^7$ time bins for the IF model (IF), and over $M_c = 3 \cdot 10^6$ spin configurations for the Ising model (BM). Error bars are given by the standard error, and are smaller or equal to the symbol size.</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>**Figure 6.** Predictive capability of the Ising model for the three-point correlation functions $T_{ijk}$. Comparison between the three-point correlation functions $T_{ijk}$ of the Ising-like model (y-axes) and the IF model network (x-axes), for system sizes $N = \{20, 40, 80, 120\}$. The blue dashed lines are the bisector $y = x$. Results are averages over $N_b = 10^7$ time bins for the IF model (IF), and over $M_c = 3 \cdot 10^6$ spin configurations for the Ising model (BM). Error bars are given by the standard error, and are smaller or equal to the symbol size.</figcaption></figure>

$K$ now counts the number of up-spins in a single configuration in the Monte Carlo evolution. We compare this quantity between the IF and Ising models in Fig. 7 for systems of size $N = \{20, 40, 80, 120\}$.

The Ising model seems to predict the $P(K)$ correctly only when $K$ is very small, of the order $K/N \lesssim 0.10$, with noticeable discrepancies seen also for this range of $K$, particularly for the probability of no activity ($K = 0$), in systems of size $N > 20$. For larger $K$, there is an intermediate range where the Ising model underestimates $P(K)$, followed by a region in which $P(K)$ is greatly overestimated. This effect seems to increase with the system size $N$. This particular behaviour was also observed when training a pairwise Ising model to reproduce the time averages of stimulated neuronal activity from a salamander retina⁹.

## Thermodynamics of the pairwise Ising models

Having mapped the neural network to a pairwise Ising model, we can start to analyse its properties. Given a system with $N$ spins, for each spin configuration we can measure the magnetization $M(\sigma) = \sum_i^N \sigma_i$ and the energy $E(\sigma) = -H(\sigma) = \sum_i^N h_i \sigma_i + \sum_i^N \sum_{j<i}^N J_{ij} \sigma_i \sigma_j$. From the fluctuations of these quantities, according to the fluctuation-dissipation theorem, we can calculate the susceptibility $\chi$ and the specific heat $C_v$ of the Ising system as a function of temperature $T$,

$$ \chi = \frac{1}{T} \cdot (\langle M^2 \rangle - \langle M \rangle^2), \qquad (7) $$

$$ C_v = \frac{1}{T^2} \cdot (\langle E^2 \rangle - \langle E \rangle^2), \qquad (8) $$

where $\langle ... \rangle$ indicates an average over spin configurations generated by Monte Carlo simulations and the Boltzmann constant $k_B$ is set equal to one. By changing $T$, one can probe the thermodynamic properties of these spin

<figure><img src="image_9.png" /><figcaption>**Figure 7.** Predictive capability of the Ising model for the probability of simultaneous firing $P(K)$. Comparison between the simultaneous firing/up-state probability $P(K)$ of the Ising-like model (red) and the IF model network (blue), for system sizes $N = \{20, 40, 80, 120\}$. Results are averaged over $N_b = 10^7$ time bins for the IF model (IF), and over $M_c = 3 \cdot 10^6$ spin configurations for the Ising model (BM). Error bars are given by the standard error, and most are smaller or equal to the symbol size.</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>Figure 6. Predictive capability of the Ising model for the three-point correlation functions $T_{ijk}$. Comparison between the three-point correlation functions $T_{ijk}$ of the Ising-like model (y-axes) and the IF model network (x-axes), for system sizes $N = \{20, 40, 80, 120\}$. The blue dashed lines are the bisector $y = x$. Results are averages over $N_b = 10^7$ time bins for the IF model (IF), and over $M_c = 3 \cdot 10^6$ spin configurations for the Ising model (BM). Error bars are given by the standard error, and are smaller or equal to the symbol size.</figcaption></figure>

$K$ now counts the number of up-spins in a single configuration in the Monte Carlo evolution. We compare this quantity between the IF and Ising models in Fig. 7 for systems of size $N = \{20, 40, 80, 120\}$.

The Ising model seems to predict the $P(K)$ correctly only when $K$ is very small, of the order $K/N \lesssim 0.10$, with noticeable discrepancies seen also for this range of $K$, particularly for the probability of no activity ($K = 0$), in systems of size $N > 20$. For larger $K$, there is an intermediate range where the Ising model underestimates $P(K)$, followed by a region in which $P(K)$ is greatly overestimated. This effect seems to increase with the system size $N$. This particular behaviour was also observed when training a pairwise Ising model to reproduce the time averages of stimulated neuronal activity from a salamander retina⁹.

## Thermodynamics of the pairwise Ising models

Having mapped the neural network to a pairwise Ising model, we can start to analyse its properties. Given a system with $N$ spins, for each spin configuration we can measure the magnetization $M(\sigma) = \sum_i^N \sigma_i$ and the energy $E(\sigma) = -H(\sigma) = \sum_i^N h_i \sigma_i + \sum_i^N \sum_{j<i}^N J_{ij} \sigma_i \sigma_j$. From the fluctuations of these quantities, according to the fluctuation-dissipation theorem, we can calculate the susceptibility $\chi$ and the specific heat $C_v$ of the Ising system as a function of temperature $T$,

$$ \chi = \frac{1}{T} \cdot (\langle M^2 \rangle - \langle M \rangle^2), \qquad (7) $$

$$ C_v = \frac{1}{T^2} \cdot (\langle E^2 \rangle - \langle E \rangle^2), \qquad (8) $$

where $\langle ... \rangle$ indicates an average over spin configurations generated by Monte Carlo simulations and the Boltzmann constant $k_B$ is set equal to one. By changing $T$, one can probe the thermodynamic properties of these spin

<figure><img src="image_9.png" /><figcaption>Figure 7. Predictive capability of the Ising model for the probability of simultaneous firing $P(K)$. Comparison between the simultaneous firing/up-state probability $P(K)$ of the Ising-like model (red) and the IF model network (blue), for system sizes $N = \{20, 40, 80, 120\}$. Results are averaged over $N_b = 10^7$ time bins for the IF model (IF), and over $M_c = 3 \cdot 10^6$ spin configurations for the Ising model (BM). Error bars are given by the standard error, and most are smaller or equal to the symbol size.</figcaption></figure>
systems. We start from a random spin configuration and thermalize it as described in the "Methods" section. Notice, however, that the MEM method only claims that the Ising model is representative of the IF model for $T = 1 \equiv T_0$, since that is the temperature used during the BM learning procedure. Since there is no evident analogy between the control parameter of the learned Ising models, $T$, and that of the IF model, $\delta u_{rec}$, temperatures $T \neq T_0$ have no obvious physical meaning in these learned Ising models, besides being a useful parameter to check if $T_0$ has a particular thermodynamic role¹⁰.

In Fig. 8 we plot the average magnetization per spin $m = \langle M \rangle/N$, the susceptibility $\chi$ and the specific heat $C_v$ as a function of the temperature $T$ for systems of size $N = \{20, 40, 80, 120\}$. For each $N$, we consider five Ising systems with distinct parameters learned from different IF networks of identical size at criticality, to assess the sample-to-sample variations for IF networks with different specific connectivities between the neurons.

The striking result is that all quantities exhibit the behavior expected close to a critical transition. The magnetization takes negative values at low temperatures, because of the majority of negative local fields, and tends to zero for high temperatures with a behavior that becomes steeper for increasing system sizes. The maxima of the susceptibility and specific heat appear near the “default” temperature $T_0 = 1$ associated with the IF model neuronal data, and their height increases with $N$, indicating that the neural networks tend to maximize thermodynamic fluctuations. This phenomenon, where fluctuations diverge with increasing system size, is a hallmark of a system operating near a critical point⁵⁵. Looking now at different network configurations of the same size $N$, we can clearly detect large variations across configurations. Particularly, for different IF networks with the same $N$, the maximum value of the susceptibility $\chi$ varies by a factor ~2 for sizes $N = \{20, 40\}$ and ~1.5 for $N = \{80, 120\}$.

Another interesting result is the cloud of random values observed for $T < 1$ for all thermodynamic quantities. This may be an indication that, for low temperatures, the Ising model transitions into a spin-glass phase. Since this phase usually exhibits complex energy landscapes⁵⁴, thermal fluctuations in this temperature regime might not be sufficient to drive the system away from the randomly chosen initial configuration and into the ground state. Indeed, starting with all $\sigma_i = -1$ as initial configuration for the Monte Carlo simulations removes this effect entirely (see Fig. S3 Supplementary Information).

## Subnetworks

In experiments, neuronal recordings are usually performed over a subset of neurons out of the total population⁹,¹³. As such, it is of practical interest to study how network subsampling might affect the analysis. A very recent study¹³ shows that the properties of firing statistics change with the spatial distribution of the subsampled neuronal patch. Taking this into account, we will consider two types of subnetworks, also studied in the aforementioned work: a subset of neurons that are spatially close together (Fig. 9a), or picked at random positions (Fig. 9b). For the former case, we select the neurons that are closest to the center of the cubic lattice.

In Fig. 10 we plot the results of the simultaneous firing probability $P(K)$ and triplets $T_{ijk}$ measured in subnetworks of $n = 40$ neurons in systems with $N = 500$ total neurons, for both spatial distributions described previously. We also compare these quantities with the predictions of the corresponding Ising models. The learning efficiency for the parameters of these models is similar to the cases using full networks, as assessed by inspecting the constrained temporal averages, $\langle \sigma_i \rangle$ and $C_{ij}$ (see Fig. S4 in Supplementary Information).

As in the case of entire networks, the Ising models fail to predict $P(K)$ for large $K$, with even larger discrepancies (Fig. 10). For the IF network, the $P(K)$ distribution decays exponentially for both subnetworks, as for the entire network case. Analogously, the probability $P(K)$ for the learned Ising models shows an initial exponential decay for low $K$, for both spatial distributions. Conversely, at $K/n \approx 0.35$, the behaviour of $P(K)$ for the Ising model is slightly different between the two spatial distributions. For the local case, the probability for observing

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

**Figure 8.** Thermodynamic functions of Ising models associated with fully-excitatory IF networks of different system sizes $N$. Average magnetization per spin $m$ (a), susceptibility $\chi$ (b) and specific heat $C_v$ (c) as a function of the temperature $T \in [0.1, 3.0]$, for systems with $N = \{20, 40, 80, 120\}$ spins. Different curves with the same colour and symbol correspond to Ising systems with parameters fitted to different IF network configurations with the same size $N$ but different specific connectivities. The vertical dashed lines indicate $T = T_0 = 1$, the “default” temperature used in the BM learning procedure to fit the respective Ising parameters to each IF network. The cloud of random values observed for $T < 1$ suggests the presence of a spin-glass phase, where the thermal energy is insufficient to drive the system away from the initial random spin configuration. Results are averages over $M_c = 3 \cdot 10^6$ spin configurations. Error bars are given by the standard error and are overall smaller or equal to the symbol size.
systems. We start from a random spin configuration and thermalize it as described in the "Methods" section. Notice, however, that the MEM method only claims that the Ising model is representative of the IF model for $T = 1 \equiv T_0$, since that is the temperature used during the BM learning procedure. Since there is no evident analogy between the control parameter of the learned Ising models, $T$, and that of the IF model, $\delta u_{rec}$, temperatures $T \neq T_0$ have no obvious physical meaning in these learned Ising models, besides being a useful parameter to check if $T_0$ has a particular thermodynamic role¹⁰.

In Fig. 8 we plot the average magnetization per spin $m = \langle M \rangle/N$, the susceptibility $\chi$ and the specific heat $C_v$ as a function of the temperature $T$ for systems of size $N = \{20, 40, 80, 120\}$. For each $N$, we consider five Ising systems with distinct parameters learned from different IF networks of identical size at criticality, to assess the sample-to-sample variations for IF networks with different specific connectivities between the neurons.

The striking result is that all quantities exhibit the behavior expected close to a critical transition. The magnetization takes negative values at low temperatures, because of the majority of negative local fields, and tends to zero for high temperatures with a behavior that becomes steeper for increasing system sizes. The maxima of the susceptibility and specific heat appear near the “default” temperature $T_0 = 1$ associated with the IF model neuronal data, and their height increases with $N$, indicating that the neural networks tend to maximize thermodynamic fluctuations. This phenomenon, where fluctuations diverge with increasing system size, is a hallmark of a system operating near a critical point⁵⁵. Looking now at different network configurations of the same size $N$, we can clearly detect large variations across configurations. Particularly, for different IF networks with the same $N$, the maximum value of the susceptibility $\chi$ varies by a factor ~2 for sizes $N = \{20, 40\}$ and ~1.5 for $N = \{80, 120\}$.

Another interesting result is the cloud of random values observed for $T < 1$ for all thermodynamic quantities. This may be an indication that, for low temperatures, the Ising model transitions into a spin-glass phase. Since this phase usually exhibits complex energy landscapes⁵⁴, thermal fluctuations in this temperature regime might not be sufficient to drive the system away from the randomly chosen initial configuration and into the ground state. Indeed, starting with all $\sigma_i = -1$ as initial configuration for the Monte Carlo simulations removes this effect entirely (see Fig. S3 Supplementary Information).

## Subnetworks

In experiments, neuronal recordings are usually performed over a subset of neurons out of the total population⁹,¹³. As such, it is of practical interest to study how network subsampling might affect the analysis. A very recent study¹³ shows that the properties of firing statistics change with the spatial distribution of the subsampled neuronal patch. Taking this into account, we will consider two types of subnetworks, also studied in the aforementioned work: a subset of neurons that are spatially close together (Fig. 9a), or picked at random positions (Fig. 9b). For the former case, we select the neurons that are closest to the center of the cubic lattice.

In Fig. 10 we plot the results of the simultaneous firing probability $P(K)$ and triplets $T_{ijk}$ measured in subnetworks of $n = 40$ neurons in systems with $N = 500$ total neurons, for both spatial distributions described previously. We also compare these quantities with the predictions of the corresponding Ising models. The learning efficiency for the parameters of these models is similar to the cases using full networks, as assessed by inspecting the constrained temporal averages, $\langle \sigma_i \rangle$ and $C_{ij}$ (see Fig. S4 in Supplementary Information).

As in the case of entire networks, the Ising models fail to predict $P(K)$ for large $K$, with even larger discrepancies (Fig. 10). For the IF network, the $P(K)$ distribution decays exponentially for both subnetworks, as for the entire network case. Analogously, the probability $P(K)$ for the learned Ising models shows an initial exponential decay for low $K$, for both spatial distributions. Conversely, at $K/n \approx 0.35$, the behaviour of $P(K)$ for the Ising model is slightly different between the two spatial distributions. For the local case, the probability for observing

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

**Figure 8.** Thermodynamic functions of Ising models associated with fully-excitatory IF networks of different system sizes $N$. Average magnetization per spin $m$ (a), susceptibility $\chi$ (b) and specific heat $C_v$ (c) as a function of the temperature $T \in [0.1, 3.0]$, for systems with $N = \{20, 40, 80, 120\}$ spins. Different curves with the same colour and symbol correspond to Ising systems with parameters fitted to different IF network configurations with the same size $N$ but different specific connectivities. The vertical dashed lines indicate $T = T_0 = 1$, the “default” temperature used in the BM learning procedure to fit the respective Ising parameters to each IF network. The cloud of random values observed for $T < 1$ suggests the presence of a spin-glass phase, where the thermal energy is insufficient to drive the system away from the initial random spin configuration. Results are averages over $M_c = 3 \cdot 10^6$ spin configurations. Error bars are given by the standard error and are overall smaller or equal to the symbol size.
<figure><img src="image_11.png" /><figcaption>(b) Random</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Figure 9.** Schematic representations of IF model subnetworks with different spatial distribution. Subnetworks (green) with neurons closely packed together (a) or picked at random positions (b). Subnetworks have $n = 40$ neurons in a system of size $N = 500$.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>Random</figcaption></figure>

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

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

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

<figure><img src="image_17.png" /><figcaption>**Figure 10.** Predictive capability of the Ising model for IF model subnetworks. Comparison between Monte Carlo sampling of the Ising model (BM) and the neural network data (IF) of the probability $P(K)$ (a,b) and triplets $T_{ijk}$ (c,d), for subnetworks with $n = 40$ neurons, with a local (a,c) and random (b,d) spatial distribution, in a system with $N = 500$ neurons. The blue dashed lines in the bottom plots are the bisector $y = x$. Results are averaged over $N_b = 10^7$ time bins for the IF model (IF), and over $M_c = 3 \cdot 10^6$ spin configurations for the Ising model (BM). Error bars are given by the standard error, and are overall smaller or equal to the symbol size.</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>(b) Random</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Figure 9.** Schematic representations of IF model subnetworks with different spatial distribution. Subnetworks (green) with neurons closely packed together (a) or picked at random positions (b). Subnetworks have $n = 40$ neurons in a system of size $N = 500$.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>Random</figcaption></figure>

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

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

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

<figure><img src="image_17.png" /><figcaption>**Figure 10.** Predictive capability of the Ising model for IF model subnetworks. Comparison between Monte Carlo sampling of the Ising model (BM) and the neural network data (IF) of the probability $P(K)$ (a,b) and triplets $T_{ijk}$ (c,d), for subnetworks with $n = 40$ neurons, with a local (a,c) and random (b,d) spatial distribution, in a system with $N = 500$ neurons. The blue dashed lines in the bottom plots are the bisector $y = x$. Results are averaged over $N_b = 10^7$ time bins for the IF model (IF), and over $M_c = 3 \cdot 10^6$ spin configurations for the Ising model (BM). Error bars are given by the standard error, and are overall smaller or equal to the symbol size.</figcaption></figure>
simultaneously up-spins reaches a second local maximum at $K/n \approx 0.75$. Notice, however, that this is a small effect, magnified by the logarithmic scale. On the other hand, for the random case, a plateau in $P(K)$ of the Ising model is observed in the intermediate regime of $K$, followed by a subsequent decay. These differences between the IF and Ising model, compared to the full network case (Fig. 7), could be due to the fact that we are neglecting influences from neurons that are not included in the subnetworks, and therefore not encoded in the Ising model parameters.

The behavior of the triplets is similar to the case of entire networks (Fig. 6), wherein the triplets are systematically overestimated by the Ising model. We note also that overall the triplets have much smaller values, when compared to the ones found for the full network with $N = 40$. Similar results were found for both $P(K)$ and the $T_{ijk}$ for subnetworks of larger size $n = 80$ in the same IF network with $N = 500$ (see Fig. S5 in Supplementary Information). The fact that both cases of spatial distributions yield Ising models with similar accuracy in reproducing the original IF model data is unexpected, as this is seemingly in contrast to recent experimental studies on mouse brains¹³, where the generalized Ising model predicted better the results from subgroups of ~100 neurons that were spatially clustered together. It should be noted, however, that the experimental study pertains to stimulated neuronal activity due to visual stimuli, whereas our model simulates spontaneous neuronal activity³⁹, in the absence of any external stimulation.

## Networks with inhibitory neurons

Inhibitory neurons hamper propagation of neuronal activity, consequently affecting the dynamics of the network⁵⁶,⁵⁷. Increasing the percentage $p_{in}$ of inhibitory neurons in a neural network moves the system into a subcritical regime⁵⁷. A way to keep the system close to the critical state when $p_{in} > 0$ is by increasing the value of the tuning parameter $\delta u_{rec}$. Here we consider systems of size $N = 80$ and $p_{in} > 0$, in the critical state (see Fig. S2 in Supplementary Information), and study how the presence of inhibition might affect the properties of the associated Ising-like models. The agreement of the average local activities and correlation functions between the Ising model and the IF networks with inhibitory neurons is as good as for fully-excitatory ones (see Fig. S6 in Supplementary Information).

In Fig. 11 we show the distributions of fields $h_i$ and interaction constants $J_{ij}$ learned from data of IF networks with $p_{in} = \{0\%, 10\%, 20\%\}$ inhibitory neurons.

The distribution of fields $h_i$ (Fig. 11a) for $p_{in} > 0$ exhibits a tail towards negative values, which becomes more pronounced as $p_{in}$ increases. Interestingly, these $h_i$ are associated with excitatory neurons (see Fig. S7 in Supplementary Information). A possible explanation is that, in the IF network, the neurons with the smallest average local activity ($\sigma_i$) are not necessarily inhibitory neurons, but rather neurons with incoming connections from inhibitory ones, which in turn are more likely to be excitatory since $p_{in} < 50\%$. This indicates that the {$h_i$} do not encode the information about whether a neuron is excitatory or inhibitory. In Fig. 11b we see that the peak of the distributions of the interaction constants $J_{ij}$, decreases with $p_{in}$, with a more pronounced tail towards negative $J_{ij}$ observed for $p_{in} = 20\%$. Interestingly, the presence of the second peak of $P(J_{ij})$ at positive $J_{ij}$ seems to be robust with respect to changes in $p_{in}$.

In Fig. 12 we plot as a function of the temperature $T$ the thermodynamic functions $m$, $\chi$ and $C_v$ of the Ising systems with the parameters of Fig. 11. For increasing $p_{in}$, the magnetization does not go to zero at high $T$ in the observed range of temperatures (Fig. 12a), indicating a higher tendency for spins to be in the down state ($\sigma_i = -1$) than for the fully excitatory case in the paramagnetic phase, which could stem from the inhibitory neurons. This is consistent with the observation of stronger negative fields $h_i$ for $p_{in} = 20\%$, as seen in Fig. 11a. In Fig. 12b,c, we see that the values of the maxima of $\chi$ and $C_v$ decrease with $p_{in}$, while their position with respect

<figure><img src="image_18.png" /><figcaption>**Figure 11.** Distributions of the learned fields $h_i$ and interaction constants $J_{ij}$ of a Ising model for neural networks with different fractions of inhibitory neurons $p_{in}$. Fields $h_i$ (a) and interaction constants $J_{ij}$ (b) associated with IF networks with $p_{in} = \{0\%, 10\%, 20\%\}$ and $N = 80$, obtained after $N_{BM} = 60,000$ iterations of the BM.</figcaption></figure>
simultaneously up-spins reaches a second local maximum at $K/n \approx 0.75$. Notice, however, that this is a small effect, magnified by the logarithmic scale. On the other hand, for the random case, a plateau in $P(K)$ of the Ising model is observed in the intermediate regime of $K$, followed by a subsequent decay. These differences between the IF and Ising model, compared to the full network case (Fig. 7), could be due to the fact that we are neglecting influences from neurons that are not included in the subnetworks, and therefore not encoded in the Ising model parameters.

The behavior of the triplets is similar to the case of entire networks (Fig. 6), wherein the triplets are systematically overestimated by the Ising model. We note also that overall the triplets have much smaller values, when compared to the ones found for the full network with $N = 40$. Similar results were found for both $P(K)$ and the $T_{ijk}$ for subnetworks of larger size $n = 80$ in the same IF network with $N = 500$ (see Fig. S5 in Supplementary Information). The fact that both cases of spatial distributions yield Ising models with similar accuracy in reproducing the original IF model data is unexpected, as this is seemingly in contrast to recent experimental studies on mouse brains¹³, where the generalized Ising model predicted better the results from subgroups of ~100 neurons that were spatially clustered together. It should be noted, however, that the experimental study pertains to stimulated neuronal activity due to visual stimuli, whereas our model simulates spontaneous neuronal activity³⁹, in the absence of any external stimulation.

## Networks with inhibitory neurons

Inhibitory neurons hamper propagation of neuronal activity, consequently affecting the dynamics of the network⁵⁶,⁵⁷. Increasing the percentage $p_{in}$ of inhibitory neurons in a neural network moves the system into a subcritical regime⁵⁷. A way to keep the system close to the critical state when $p_{in} > 0$ is by increasing the value of the tuning parameter $\delta u_{rec}$. Here we consider systems of size $N = 80$ and $p_{in} > 0$, in the critical state (see Fig. S2 in Supplementary Information), and study how the presence of inhibition might affect the properties of the associated Ising-like models. The agreement of the average local activities and correlation functions between the Ising model and the IF networks with inhibitory neurons is as good as for fully-excitatory ones (see Fig. S6 in Supplementary Information).

In Fig. 11 we show the distributions of fields $h_i$ and interaction constants $J_{ij}$ learned from data of IF networks with $p_{in} = \{0\%, 10\%, 20\%\}$ inhibitory neurons.

The distribution of fields $h_i$ (Fig. 11a) for $p_{in} > 0$ exhibits a tail towards negative values, which becomes more pronounced as $p_{in}$ increases. Interestingly, these $h_i$ are associated with excitatory neurons (see Fig. S7 in Supplementary Information). A possible explanation is that, in the IF network, the neurons with the smallest average local activity ($\sigma_i$) are not necessarily inhibitory neurons, but rather neurons with incoming connections from inhibitory ones, which in turn are more likely to be excitatory since $p_{in} < 50\%$. This indicates that the {$h_i$} do not encode the information about whether a neuron is excitatory or inhibitory. In Fig. 11b we see that the peak of the distributions of the interaction constants $J_{ij}$, decreases with $p_{in}$, with a more pronounced tail towards negative $J_{ij}$ observed for $p_{in} = 20\%$. Interestingly, the presence of the second peak of $P(J_{ij})$ at positive $J_{ij}$ seems to be robust with respect to changes in $p_{in}$.

In Fig. 12 we plot as a function of the temperature $T$ the thermodynamic functions $m$, $\chi$ and $C_v$ of the Ising systems with the parameters of Fig. 11. For increasing $p_{in}$, the magnetization does not go to zero at high $T$ in the observed range of temperatures (Fig. 12a), indicating a higher tendency for spins to be in the down state ($\sigma_i = -1$) than for the fully excitatory case in the paramagnetic phase, which could stem from the inhibitory neurons. This is consistent with the observation of stronger negative fields $h_i$ for $p_{in} = 20\%$, as seen in Fig. 11a. In Fig. 12b,c, we see that the values of the maxima of $\chi$ and $C_v$ decrease with $p_{in}$, while their position with respect

<figure><img src="image_18.png" /><figcaption>**Figure 11.** Distributions of the learned fields $h_i$ and interaction constants $J_{ij}$ of a Ising model for neural networks with different fractions of inhibitory neurons $p_{in}$. Fields $h_i$ (a) and interaction constants $J_{ij}$ (b) associated with IF networks with $p_{in} = \{0\%, 10\%, 20\%\}$ and $N = 80$, obtained after $N_{BM} = 60,000$ iterations of the BM.</figcaption></figure>
<figure><img src="image_19.png" /><figcaption>**Figure 12.** Thermodynamic functions of Ising-like models associated with IF model networks with different fractions of inhibitory neurons $p_{in}$. Average magnetization per spin $m$ (a), susceptibility $χ$ (b) and specific heat $C_v$ (c) as a function of the temperature $T \in [0.1, 3.0]$, simulated using the learned parameters shown in Fig. 11 for different fractions of inhibitory neurons $p_{in} = \{0\%, 10\%, 20\%\}$ and $N = 80$. As in Fig. 8, the cloud of random values for $T < 1$ suggests the presence of a spin-glass phase. Results are averaged over $M_c = 3 \cdot 10^6$ spin configurations. Error bars are given by the standard error and are overall smaller or equal to the symbol size.</figcaption></figure>

to the temperature remains unchanged. This suggests that introducing inhibition in the IF networks reduces thermal fluctuations in the associated Ising model.

## Partially-connected pairwise Ising models

One of the main disadvantages of BM learning is its intense CPU time demand³⁵,⁵⁸, restricting the network sizes one could potentially analyse. To try to circumvent this problem, we will consider Ising models with pruned links²⁸, whereby we remove the couplings associated with the weakest correlated pairs of neurons. Specifically, we set $J_{ij} = 0$ if the corresponding $C_{ij}^{(\text{IF})}$ is below a certain threshold, defined as a fraction of the largest measured correlation function, i.e. if $C_{ij}^{(\text{IF})} < \eta \max(C_{ij}^{(\text{IF})})$, where $\eta \in [0, 1]$ sets the threshold. Thus, the BM only needs to learn a subset of the total number of couplings $J_{ij}$, possibly accelerating the convergence process and consequently allowing the study of networks larger than the ones considered so far. Furthermore, an obvious speed-up is also achieved by making use of the fact that a subset of the $\{J_{ij}\}$ are zero, and therefore can be disregarded in the double sum of Eq. (5) during the Monte Carlo simulations. We will consider three thresholds $\eta \in \{0.10, 0.15, 0.20\}$, using a fully-excitatory system of size $N = 180$ in the critical state. Computing only a fraction of the interaction terms allowed to obtain the BM results for the $N = 180$ system in a less or comparable CPU time with respect to the one required by the $N = 120$ systems where we considered the full set $\{J_{ij}\}$, with a greater CPU speed-up achieved the larger the threshold $\eta$ is.

<figure><img src="image_20.png" /><figcaption>**Figure 13.** Quality test for the BM learning process with a partially-connected Ising model. Comparison between the correlation functions $C_{ij}$ of the partially-connected Ising-like model (y-axes) and IF model network (x-axes), for a fully-excitatory system at criticality of size $N = 180$ and three different thresholds $\eta = \{0.10, 0.15, 0.20\}$ for the removal of the $J_{ij}$. If $C_{ij}^{(\text{IF})} < \eta \max(C_{ij}^{(\text{IF})})$ (grey dots), the corresponding interaction constant is set to $J_{ij} = 0$. The vertical dashed lines indicate the value $\eta \max(C_{ij}^{(\text{IF})})$, with the corresponding approximate percentage of removed $J_{ij}$ reported at the right of this line. The blue dashed lines are the bisector $y = x$. Results are averaged over $N_b = 10^7$ time bins for the IF model (IF), and over $M_c = 3 \cdot 10^7$ spin configurations for the Ising model (BM). Error bars are given by the standard error, and are smaller or equal to the symbol size.</figcaption></figure>
<figure><img src="image_19.png" /><figcaption>**Figure 12.** Thermodynamic functions of Ising-like models associated with IF model networks with different fractions of inhibitory neurons $p_{in}$. Average magnetization per spin $m$ (a), susceptibility $χ$ (b) and specific heat $C_v$ (c) as a function of the temperature $T \in [0.1, 3.0]$, simulated using the learned parameters shown in Fig. 11 for different fractions of inhibitory neurons $p_{in} = \{0\%, 10\%, 20\%\}$ and $N = 80$. As in Fig. 8, the cloud of random values for $T < 1$ suggests the presence of a spin-glass phase. Results are averaged over $M_c = 3 \cdot 10^6$ spin configurations. Error bars are given by the standard error and are overall smaller or equal to the symbol size.</figcaption></figure>

to the temperature remains unchanged. This suggests that introducing inhibition in the IF networks reduces thermal fluctuations in the associated Ising model.

## Partially-connected pairwise Ising models

One of the main disadvantages of BM learning is its intense CPU time demand³⁵,⁵⁸, restricting the network sizes one could potentially analyse. To try to circumvent this problem, we will consider Ising models with pruned links²⁸, whereby we remove the couplings associated with the weakest correlated pairs of neurons. Specifically, we set $J_{ij} = 0$ if the corresponding $C_{ij}^{(\text{IF})}$ is below a certain threshold, defined as a fraction of the largest measured correlation function, i.e. if $C_{ij}^{(\text{IF})} < \eta \max(C_{ij}^{(\text{IF})})$, where $\eta \in [0, 1]$ sets the threshold. Thus, the BM only needs to learn a subset of the total number of couplings $J_{ij}$, possibly accelerating the convergence process and consequently allowing the study of networks larger than the ones considered so far. Furthermore, an obvious speed-up is also achieved by making use of the fact that a subset of the $\{J_{ij}\}$ are zero, and therefore can be disregarded in the double sum of Eq. (5) during the Monte Carlo simulations. We will consider three thresholds $\eta \in \{0.10, 0.15, 0.20\}$, using a fully-excitatory system of size $N = 180$ in the critical state. Computing only a fraction of the interaction terms allowed to obtain the BM results for the $N = 180$ system in a less or comparable CPU time with respect to the one required by the $N = 120$ systems where we considered the full set $\{J_{ij}\}$, with a greater CPU speed-up achieved the larger the threshold $\eta$ is.

<figure><img src="image_20.png" /><figcaption>**Figure 13.** Quality test for the BM learning process with a partially-connected Ising model. Comparison between the correlation functions $C_{ij}$ of the partially-connected Ising-like model (y-axes) and IF model network (x-axes), for a fully-excitatory system at criticality of size $N = 180$ and three different thresholds $\eta = \{0.10, 0.15, 0.20\}$ for the removal of the $J_{ij}$. If $C_{ij}^{(\text{IF})} < \eta \max(C_{ij}^{(\text{IF})})$ (grey dots), the corresponding interaction constant is set to $J_{ij} = 0$. The vertical dashed lines indicate the value $\eta \max(C_{ij}^{(\text{IF})})$, with the corresponding approximate percentage of removed $J_{ij}$ reported at the right of this line. The blue dashed lines are the bisector $y = x$. Results are averaged over $N_b = 10^7$ time bins for the IF model (IF), and over $M_c = 3 \cdot 10^7$ spin configurations for the Ising model (BM). Error bars are given by the standard error, and are smaller or equal to the symbol size.</figcaption></figure>
In Fig. 13 we present the comparison between the correlation functions of the IF model and the partially-connected Ising model for the three different thresholds. We see that, even though more than 40% couplings $J_{ij}$ have been removed in the Ising model for the lowest threshold considered, $\eta = 0.10$, and more than 70% for the largest one, $\eta = 0.20$, the overall data are still reconstructed by the partially-connected Ising model, with a decrease in quality of the fit for the lowest $C_{ij}$ (grey dots in the plots) since they were not considered in the learning process of the BM. On the other hand, all $\langle\sigma_i\rangle$ are well predicted by the Ising model (see Fig. S8 in Supplementary Information) since we fit all the $N = 180$ fields $\{h_i\}$.

In Fig. 14 we plot the learned fields and the distributions of the learned non-zero interaction constants of the partially-connected Ising model for the three different thresholds $\eta$. While the set of fields $\{h_i\}$ remains qualitatively similar for the three thresholds $\eta$, with a majority of negative fields as in the case of the fully-connected Ising models, the distributions of the non-zero $J_{ij}$ change significantly with $\eta$, with the peak at positive $J_{ij}$ increasing and the one at $J_{ij} \approx 0$ decreasing as we remove progressively more $J_{ij}$. Since we remove only the $J_{ij}$ associated with the smallest correlation functions $C_{ij}$, this seemingly indicates that the second peak at $J_{ij} > 0$ also seen in the distributions for the fully-connected Ising models (Fig. 4) is associated with the subset of the largest $C_{ij}$.

As usual, we can also analyze how the Ising model might predict quantities that are not being constrained by the BM algorithm. In Fig. 15 we present the comparison between the tree-point correlation functions $T_{ijk}$ measured in the partially-connected Ising model and the IF model. As in the fully-connected case, the Ising model consistently overestimates the values of $T_{ijk}$. The quality of the fit for the smallest threshold considered, $\eta = 0.10$, is comparable to the fully-connected Ising models (Fig. 6). However, as $\eta$ increases, the quality of the fit worsens, particularly for the smallest $T_{ijk}^{(IF)}$.

Next, we present in Fig. 16 the results for the probability $P(K)$ of simultaneous firing. The predictive capability of the partially-connected Ising model for small $K$ is similar to that of fully-connected ones (Fig. 7), while significant differences can be seen for large $K$. As the threshold $\eta$ is increased, the Ising model seems to predict larger simultaneous activity at the scale of the whole network, where $K \approx N$, indicated by a local maximum in $P(K)$. This seems to indicate that the small $J_{ij}$ that we are disregarding encode the information concerning the large $K$ regime of $P(K)$.

Finally, in Fig. 17 we plot the thermodynamic functions $m$, $\chi$ and $C_v$ versus temperature $T$ for the partially-connected Ising systems, considering the three cases $\eta = \{0.10, 0.15, 0.20\}$. For all quantities results are quite robust with respect to $\eta$, taking into account that there is a difference of $\approx 31.9\%$ of removed $J_{ij}$ between the smallest and largest thresholds $\eta$ considered. The tendency for the susceptibility $\chi$ and specific heat $C_v$ to diverge

<figure><img src="image_21.png" /><figcaption>**Figure 14.** Sets of learned fields $h_i$ and non-zero interaction constants $J_{ij}$ of the partially-connected Ising model. Plots of the fields $h_i$ (a–c), sorted by the average local activity $\langle\sigma_i\rangle^{(IF)}$ of the associated neuron $i$, in order of decreasing $\langle\sigma_i\rangle^{(IF)}$, and distributions of the non-zero interaction constants $J_{ij}$ (d–f), that reproduce the respective data of the Ising model presented in Fig. 13 for the three different thresholds $\eta = \{0.10, 0.15, 0.20\}$, for a fully-excitatory system at criticality of size $N = 180$.</figcaption></figure>
In Fig. 13 we present the comparison between the correlation functions of the IF model and the partially-connected Ising model for the three different thresholds. We see that, even though more than 40% couplings $J_{ij}$ have been removed in the Ising model for the lowest threshold considered, $\eta = 0.10$, and more than 70% for the largest one, $\eta = 0.20$, the overall data are still reconstructed by the partially-connected Ising model, with a decrease in quality of the fit for the lowest $C_{ij}$ (grey dots in the plots) since they were not considered in the learning process of the BM. On the other hand, all $\langle\sigma_i\rangle$ are well predicted by the Ising model (see Fig. S8 in Supplementary Information) since we fit all the $N = 180$ fields $\{h_i\}$.

In Fig. 14 we plot the learned fields and the distributions of the learned non-zero interaction constants of the partially-connected Ising model for the three different thresholds $\eta$. While the set of fields $\{h_i\}$ remains qualitatively similar for the three thresholds $\eta$, with a majority of negative fields as in the case of the fully-connected Ising models, the distributions of the non-zero $J_{ij}$ change significantly with $\eta$, with the peak at positive $J_{ij}$ increasing and the one at $J_{ij} \approx 0$ decreasing as we remove progressively more $J_{ij}$. Since we remove only the $J_{ij}$ associated with the smallest correlation functions $C_{ij}$, this seemingly indicates that the second peak at $J_{ij} > 0$ also seen in the distributions for the fully-connected Ising models (Fig. 4) is associated with the subset of the largest $C_{ij}$.

As usual, we can also analyze how the Ising model might predict quantities that are not being constrained by the BM algorithm. In Fig. 15 we present the comparison between the tree-point correlation functions $T_{ijk}$ measured in the partially-connected Ising model and the IF model. As in the fully-connected case, the Ising model consistently overestimates the values of $T_{ijk}$. The quality of the fit for the smallest threshold considered, $\eta = 0.10$, is comparable to the fully-connected Ising models (Fig. 6). However, as $\eta$ increases, the quality of the fit worsens, particularly for the smallest $T_{ijk}^{(IF)}$.

Next, we present in Fig. 16 the results for the probability $P(K)$ of simultaneous firing. The predictive capability of the partially-connected Ising model for small $K$ is similar to that of fully-connected ones (Fig. 7), while significant differences can be seen for large $K$. As the threshold $\eta$ is increased, the Ising model seems to predict larger simultaneous activity at the scale of the whole network, where $K \approx N$, indicated by a local maximum in $P(K)$. This seems to indicate that the small $J_{ij}$ that we are disregarding encode the information concerning the large $K$ regime of $P(K)$.

Finally, in Fig. 17 we plot the thermodynamic functions $m$, $\chi$ and $C_v$ versus temperature $T$ for the partially-connected Ising systems, considering the three cases $\eta = \{0.10, 0.15, 0.20\}$. For all quantities results are quite robust with respect to $\eta$, taking into account that there is a difference of $\approx 31.9\%$ of removed $J_{ij}$ between the smallest and largest thresholds $\eta$ considered. The tendency for the susceptibility $\chi$ and specific heat $C_v$ to diverge

<figure><img src="image_21.png" /><figcaption>**Figure 14.** Sets of learned fields $h_i$ and non-zero interaction constants $J_{ij}$ of the partially-connected Ising model. Plots of the fields $h_i$ (a–c), sorted by the average local activity $\langle\sigma_i\rangle^{(IF)}$ of the associated neuron $i$, in order of decreasing $\langle\sigma_i\rangle^{(IF)}$, and distributions of the non-zero interaction constants $J_{ij}$ (d–f), that reproduce the respective data of the Ising model presented in Fig. 13 for the three different thresholds $\eta = \{0.10, 0.15, 0.20\}$, for a fully-excitatory system at criticality of size $N = 180$.</figcaption></figure>
<figure><img src="image_22.png" /><figcaption>**Figure 15.** Predictive capability of the partially-connected Ising model for the three-point correlation functions $T_{ijk}$. Comparison between the three-point correlation functions $T_{ijk}$ of the partially-connected Ising-like model (y-axes) and the IF network (x-axes), for the three thresholds $\eta = \{0.10, 0.15, 0.20\}$, for a fully-excitatory system at criticality of size $N = 180$. Results are averaged over $N_b = 10^7$ time bins for the IF model (IF), and over $M_c = 3 \cdot 10^7$ spin configurations for the Ising model (BM). Error bars are given by the standard error, and most are smaller or equal to the symbol size.</figcaption></figure>

<figure><img src="image_23.png" /><figcaption>**Figure 16.** Predictive capability of the partially-connected Ising model for the probability of simultaneous firing $P(K)$. Comparison between the simultaneous firing/up-state probability $P(K)$ of the Ising-like model (red) and the IF model network (blue), for the three thresholds $\eta = \{0.10, 0.15, 0.20\}$, for a fully-excitatory system at criticality of size $N = 180$. Results are averaged over $N_b = 10^7$ time bins for the IF model (IF), and over $M_c = 3 \cdot 10^7$ spin configurations for the Ising model (BM). Error bars are given by the standard error, and most are smaller or equal to the symbol size.</figcaption></figure>

with $N$ is also still clearly visible when comparing to the results for a smaller system size with $N = 120$ (black symbols), which in turn considers the full set $\{J_{ij}\}$.

Discussion

In this study we apply for the first time the Maximum Entropy Modelling (MEM) method to neuronal activity data generated by a numerical model containing the fundamental biological features of living networks and able to tune the system at criticality. Indeed, it has been widely proposed in the literature that the brain could be considered as a system acting close to a critical point, however in the majority of previous papers this feature was not taken in account clearly. The advantage of the numerical study is the possibility to have a clear knowledge of the state of the neuronal system, to control the percentage of inhibitory neurons and to implement artificially, by considering subnetworks, the subsampling limitation in experimental measurements. We mapped the local and pairwise information of IF complex neural networks into Ising models with frustrated spins. Independently of the system size $N$, the local fields $h_i$ are mostly negative (Fig. 4a–d) and the distribution of interaction constants $J_{ij}$ (Fig. 4e–h) is bimodal for small systems, but tends to a normal distribution centered around zero as $N$
<figure><img src="image_22.png" /><figcaption>**Figure 15.** Predictive capability of the partially-connected Ising model for the three-point correlation functions $T_{ijk}$. Comparison between the three-point correlation functions $T_{ijk}$ of the partially-connected Ising-like model (y-axes) and the IF network (x-axes), for the three thresholds $\eta = \{0.10, 0.15, 0.20\}$, for a fully-excitatory system at criticality of size $N = 180$. Results are averaged over $N_b = 10^7$ time bins for the IF model (IF), and over $M_c = 3 \cdot 10^7$ spin configurations for the Ising model (BM). Error bars are given by the standard error, and most are smaller or equal to the symbol size.</figcaption></figure>

<figure><img src="image_23.png" /><figcaption>**Figure 16.** Predictive capability of the partially-connected Ising model for the probability of simultaneous firing $P(K)$. Comparison between the simultaneous firing/up-state probability $P(K)$ of the Ising-like model (red) and the IF model network (blue), for the three thresholds $\eta = \{0.10, 0.15, 0.20\}$, for a fully-excitatory system at criticality of size $N = 180$. Results are averaged over $N_b = 10^7$ time bins for the IF model (IF), and over $M_c = 3 \cdot 10^7$ spin configurations for the Ising model (BM). Error bars are given by the standard error, and most are smaller or equal to the symbol size.</figcaption></figure>

with $N$ is also still clearly visible when comparing to the results for a smaller system size with $N = 120$ (black symbols), which in turn considers the full set $\{J_{ij}\}$.

Discussion

In this study we apply for the first time the Maximum Entropy Modelling (MEM) method to neuronal activity data generated by a numerical model containing the fundamental biological features of living networks and able to tune the system at criticality. Indeed, it has been widely proposed in the literature that the brain could be considered as a system acting close to a critical point, however in the majority of previous papers this feature was not taken in account clearly. The advantage of the numerical study is the possibility to have a clear knowledge of the state of the neuronal system, to control the percentage of inhibitory neurons and to implement artificially, by considering subnetworks, the subsampling limitation in experimental measurements. We mapped the local and pairwise information of IF complex neural networks into Ising models with frustrated spins. Independently of the system size $N$, the local fields $h_i$ are mostly negative (Fig. 4a–d) and the distribution of interaction constants $J_{ij}$ (Fig. 4e–h) is bimodal for small systems, but tends to a normal distribution centered around zero as $N$
<figure><img src="image_24.png" /></figure>

**Figure 17.** Thermodynamic functions of partially-connected Ising models associated with an IF model network with $N = 180$. Average magnetization per spin $m$ (a), susceptibility $\chi$ (b) and specific heat $C_v$ (c) as a function of the temperature $T \in [0.1, 3.0]$, simulated using the learned parameters shown in Fig. 14 for the three different thresholds $\eta = \{0.10, 0.15, 0.20\}$, for a fully-excitatory IF network at criticality of size $N = 180$. The black symbols show the results for a learned Ising model associated with an IF network of size $N = 120$, considering the full set $\{J_{ij}\}$. As in Figs. 8 and 12, the cloud of random values for $T < 1$ suggests the presence of a spin-glass phase. Results are averaged over $M_c = 3 \cdot 10^6$ spin configurations. Error bars are given by the standard error and are overall smaller or equal to the symbol size.

increases. These Ising systems display a spin glass phase at low temperatures, independently of $N$ (Fig. 8), and the susceptibility and specific heat tend to diverge with the system size, with a maximum near $T = T_0 = 1$, the effective temperature used in the BM learning to fit the IF network neuronal data. This is an indication that the Ising model analogs of the IF networks operate near a critical point. Furthermore, at least for $N \le 120$, the height of the maximum in the susceptibility is sensible to changes in the details of the connectivity of the network (Fig. 8b). We remark that the Ising models do not predict well unconstrained quantities such as the three-point correlation functions (Fig. 6) and the probability of simultaneous firing even for whole networks (Fig. 7). This discrepancy is enhanced when considering only a subset of the neurons (Fig. 10), likely due to inputs from neurons outside the subnetwork which are not considered in the MEM mapping. This indicates that additional caution should be taken when analysing real neural networks using the MEM approach, which often consider only a subset of the total neuronal population. The presence of inhibitory neurons in the IF networks leads to reduced thermal fluctuations in the associated Ising models, evidenced by a decrease of the maxima near the critical point (Fig. 12).

We have to stress that, as verified in the present case, the potential of the MEM approach is limited by its CPU time demand when considering large system sizes, as the computation time complexity increases as $\propto N^2$ since the Ising model is fully-connected and has high heterogeneity in the distributions of the fields $h_i$ and interaction constants $J_{ij}$. To circumvent this issue, we considered a partially-connected Ising model with only a subset of the total pairs of interaction constants $\{J_{ij}\}$ obtained considering only the largest correlation functions measured in the IF network, allowing to study a system of size $N = 180$ (Figs. 13, 14, 15, 16 and 17). Similar MEM procedures allow the study of much larger sizes, such as the so-called restricted Boltzmann Machine⁵⁹ or a random projection model⁶⁰, but the resulting maximum entropy distributions from these approaches are no longer analogous to that of an Ising model, and the thermodynamical interpretation is no longer applicable³⁵. A very recent study⁶¹ shows that it might be possible to consider large networks considering only connections which contain the maximum mutual information among all pairs of neurons. As future developments, we remark that the IF model allows the investigation of systems of off criticality by appropriately tuning the short-term plasticity parameter. The study of neuronal systems at and off criticality has recently revealed intriguing behavior typical of thermodynamic systems, as the Fluctuation-Dissipation relations⁶²,⁶³. This approach can then be used to study thermodynamic properties of networks in the sub- and supercritical state, brain states associated with pathological conditions⁶⁴.

# Methods

## Integrate-and-fire model

The model considers $N$ neurons placed randomly inside a cubic space of side $L = \sqrt{N/\rho}$, where $\rho = 0.016$⁵⁷ is the density of the neurons. Connections are directed and weighted, with dynamic synaptic strengths $w_{ij} \in [0, 1]$ between pre- and post-synaptic neurons $i$ and $j$. Each neuron has at least one incoming connection, and the distribution of out-going degrees $k_{\text{out}} \in [2, 20]$ is a power law, i.e. $P(k_{\text{out}}) \propto k_{\text{out}}^{-2}$, following experimental measurements in functional networks⁶⁵. The probability that two neurons are connected decays exponentially with their euclidean distance $r$, $P(r) \propto e^{-r/r_0}$, where $r_0 = 5$ is a characteristic length⁶⁶. Neurons can be excitatory or inhibitory, with fraction $p_{\text{in}}$. We implement synaptic plasticity, in which the strengths $w_{ij}$ can change dynamically with time. We consider both short-term and long-term plasticity (STP and LTP for short, respectively), so we separate the synaptic strengths $w_{ij}(t) = u_i(t)g_{ij}$ into the two components $u_i(t)$ and $g_{ij}$. STP refers to the modification of synaptic strengths on a short timescale, of the order of milliseconds, and we denote by $u_i(t) \in [0, 1]$ the normalized amount of neurotransmitters available to neuron $i$ at time $t$, representing the so-called readily releasable pool of neurotransmitter vesicles⁶⁷. On the other hand, LTP is a Hebbian-like process which strengthens or weakens connections strengths $g_{ij} \in (0, 1]$ depending on their usage over time, acting on much longer timescales, ranging from minutes to hours or even years⁶⁸, so we regard this term as constant compared to the timescale of the dynamics. Each neuron $i$ is characterized by a membrane potential $v_i$. A neuron $i$ will fire at
<figure><img src="image_24.png" /></figure>

**Figure 17.** Thermodynamic functions of partially-connected Ising models associated with an IF model network with $N = 180$. Average magnetization per spin $m$ (a), susceptibility $\chi$ (b) and specific heat $C_v$ (c) as a function of the temperature $T \in [0.1, 3.0]$, simulated using the learned parameters shown in Fig. 14 for the three different thresholds $\eta = \{0.10, 0.15, 0.20\}$, for a fully-excitatory IF network at criticality of size $N = 180$. The black symbols show the results for a learned Ising model associated with an IF network of size $N = 120$, considering the full set $\{J_{ij}\}$. As in Figs. 8 and 12, the cloud of random values for $T < 1$ suggests the presence of a spin-glass phase. Results are averaged over $M_c = 3 \cdot 10^6$ spin configurations. Error bars are given by the standard error and are overall smaller or equal to the symbol size.

increases. These Ising systems display a spin glass phase at low temperatures, independently of $N$ (Fig. 8), and the susceptibility and specific heat tend to diverge with the system size, with a maximum near $T = T_0 = 1$, the effective temperature used in the BM learning to fit the IF network neuronal data. This is an indication that the Ising model analogs of the IF networks operate near a critical point. Furthermore, at least for $N \le 120$, the height of the maximum in the susceptibility is sensible to changes in the details of the connectivity of the network (Fig. 8b). We remark that the Ising models do not predict well unconstrained quantities such as the three-point correlation functions (Fig. 6) and the probability of simultaneous firing even for whole networks (Fig. 7). This discrepancy is enhanced when considering only a subset of the neurons (Fig. 10), likely due to inputs from neurons outside the subnetwork which are not considered in the MEM mapping. This indicates that additional caution should be taken when analysing real neural networks using the MEM approach, which often consider only a subset of the total neuronal population. The presence of inhibitory neurons in the IF networks leads to reduced thermal fluctuations in the associated Ising models, evidenced by a decrease of the maxima near the critical point (Fig. 12).

We have to stress that, as verified in the present case, the potential of the MEM approach is limited by its CPU time demand when considering large system sizes, as the computation time complexity increases as $\propto N^2$ since the Ising model is fully-connected and has high heterogeneity in the distributions of the fields $h_i$ and interaction constants $J_{ij}$. To circumvent this issue, we considered a partially-connected Ising model with only a subset of the total pairs of interaction constants $\{J_{ij}\}$ obtained considering only the largest correlation functions measured in the IF network, allowing to study a system of size $N = 180$ (Figs. 13, 14, 15, 16 and 17). Similar MEM procedures allow the study of much larger sizes, such as the so-called restricted Boltzmann Machine⁵⁹ or a random projection model⁶⁰, but the resulting maximum entropy distributions from these approaches are no longer analogous to that of an Ising model, and the thermodynamical interpretation is no longer applicable³⁵. A very recent study⁶¹ shows that it might be possible to consider large networks considering only connections which contain the maximum mutual information among all pairs of neurons. As future developments, we remark that the IF model allows the investigation of systems of off criticality by appropriately tuning the short-term plasticity parameter. The study of neuronal systems at and off criticality has recently revealed intriguing behavior typical of thermodynamic systems, as the Fluctuation-Dissipation relations⁶²,⁶³. This approach can then be used to study thermodynamic properties of networks in the sub- and supercritical state, brain states associated with pathological conditions⁶⁴.

# Methods

## Integrate-and-fire model

The model considers $N$ neurons placed randomly inside a cubic space of side $L = \sqrt{N/\rho}$, where $\rho = 0.016$⁵⁷ is the density of the neurons. Connections are directed and weighted, with dynamic synaptic strengths $w_{ij} \in [0, 1]$ between pre- and post-synaptic neurons $i$ and $j$. Each neuron has at least one incoming connection, and the distribution of out-going degrees $k_{\text{out}} \in [2, 20]$ is a power law, i.e. $P(k_{\text{out}}) \propto k_{\text{out}}^{-2}$, following experimental measurements in functional networks⁶⁵. The probability that two neurons are connected decays exponentially with their euclidean distance $r$, $P(r) \propto e^{-r/r_0}$, where $r_0 = 5$ is a characteristic length⁶⁶. Neurons can be excitatory or inhibitory, with fraction $p_{\text{in}}$. We implement synaptic plasticity, in which the strengths $w_{ij}$ can change dynamically with time. We consider both short-term and long-term plasticity (STP and LTP for short, respectively), so we separate the synaptic strengths $w_{ij}(t) = u_i(t)g_{ij}$ into the two components $u_i(t)$ and $g_{ij}$. STP refers to the modification of synaptic strengths on a short timescale, of the order of milliseconds, and we denote by $u_i(t) \in [0, 1]$ the normalized amount of neurotransmitters available to neuron $i$ at time $t$, representing the so-called readily releasable pool of neurotransmitter vesicles⁶⁷. On the other hand, LTP is a Hebbian-like process which strengthens or weakens connections strengths $g_{ij} \in (0, 1]$ depending on their usage over time, acting on much longer timescales, ranging from minutes to hours or even years⁶⁸, so we regard this term as constant compared to the timescale of the dynamics. Each neuron $i$ is characterized by a membrane potential $v_i$. A neuron $i$ will fire at
some time $t$ when its potential $v_i$ surpasses a threshold $v_c = 1$. Activity will then propagate to all post-synaptic neurons $j$ with incoming connections from $i$ according to the following equations³⁹:

$$v_j(t+1) = v_j(t) \pm v_i(t) u_i(t) g_{ij}, \qquad (9)$$

$$u_i(t+1) = u_i(t) \cdot (1 - \delta u), \qquad (10)$$

$$v_i(t+1) = 0, \qquad (11)$$

where + and – stands for excitatory and inhibitory pre-synaptic neuron, respectively, and $\delta u = 0.05$ represents the fractional amount of neurotransmitter released at each neuronal firing⁶⁹. The timestep unit corresponds roughly to the joint interval of synaptic and axonal delay, i.e. the time interval between the generation of the action potential at the pre-synaptic neuron and the membrane potential change at the post-synaptic one, and is of the order of 10 milliseconds⁷⁰. We set a minimum value $v_{min} = -1$ for the membrane potential of each neuron, to prevent the possibility of a $v_i$ being systematically decreased to overly negative values when $p_{in} > 0$. After firing, a neuron enters into a refractory period of a single timestep $t_r = 1$ during which it is incapable of receiving or eliciting any activity. To keep the activity ongoing, we implement a small external stimulation. Namely, at every timestep, even during avalanches, a voltage input $\delta v = 0.1 v_c$ is added to a randomly chosen neuron. From Eq. (10) we see that the amount of synaptic resources gradually depletes as a neuron fires, eventually rendering it incapable of transmitting further signals. In real systems, this is counteracted by a recovery mechanism, where the readily releasable pool is slowly recharged over the span of seconds⁷¹. We assume a separation of timescales⁷², and recharge simultaneously the $u_i$ of all neurons by a certain amount $\delta u_{rec}$ only at the end of every avalanche, $u_i(t) \rightarrow u_i(t) + \delta u_{rec}$. By changing the value of $\delta u_{rec}$ one can adjust the dynamical state of the system³⁹,⁵⁷. Before performing any measurements, we let the dynamics evolve for a certain number of timesteps in order to shape the distribution of strengths $g_{ij}$ by LTP. We start by setting them uniformly distributed in the interval $g_{ij} \in [0.04, 0.06]$. According to the rules of Hebbian plasticity⁷³,⁷⁴, if some neuron $i$ frequently stimulates another neuron $j$, then the synapse from $i$ to $j$ will be strengthened. In this case, we increase the strength of the synapses $g_{ij}$ proportionally to the voltage variation induced in the post-synaptic neuron $j$ due to $i$ as $g_{ij}(t+1) = g_{ij}(t) + \delta g_j$, where $\delta g_j = \beta |v_j(t+1) - v_i(t)|$ and $\beta = 0.04$ sets the rate of this adaptation. On the other hand, synapses that are rarely active tend to weaken over time⁷⁴. Therefore, at the end of each avalanche, we decrease all terms $g_{ij}$ by the average increase in strength per synapse, $g_{ij}(t+1) = g_{ij}(t) - \frac{1}{N_s} \sum \delta g_j$ where $N_s$ is the number of synapses. In real networks, synapses that weaken consistently are eventually pruned. To avoid modifying the scale-free structure of the network, We modify the strengths $g_{ij}$ either for a fixed number $N_{aval} = 10^4$ of avalanches or until a strength $g_{ij}$ first reaches a minimum value $g_{min} = 10^{-5}$, where we then set that strength to $g_{ij} = g_{min}$.

## Maximum entropy modelling

The binarization of the IF model firing dynamics introduces the notion of the probability $P_{IF}(\sigma)$ to observe, during a time bin of duration $\Delta t_b$, any of the $2^N$ possible patterns of firing states $\sigma = \{\sigma_1, \sigma_2, ..., \sigma_N\}$ in a network of size $N$, with each $\sigma_i \in \{-1, 1\}$. We are interested in defining this $P_{IF}(\sigma)$ in a way that is consistent with the expectation values of the average local activities $\langle \sigma_i \rangle$ and two-point activities $\langle \sigma_i \sigma_j \rangle$ measured in the IF model for a given network. This surmounts to finding a probability distribution $P(\sigma)$ that maximizes the entropy¹⁴ $S = - \sum_{\sigma} P(\sigma) \ln[P(\sigma)]$, where $\sum_{\sigma}$ indicates a sum over all possible outcomes of $\sigma$, while subject to the constraints $\langle \sigma_i \rangle = \sum_{\sigma} \sigma_i P(\sigma)$ and $\langle \sigma_i \sigma_j \rangle = \sum_{\sigma} \sigma_i \sigma_j P(\sigma)$. Solving this problem using the method of Lagrangian multipliers¹⁵ (see Supplementary Information for the derivation) yields the distribution of spin states of a generalized Ising model at unit temperature⁷,

$$P(\sigma) = \frac{1}{Z} e^{-H(\sigma)}, \qquad (12)$$

$$Z = \sum_{\sigma} e^{-H(\sigma)}, \qquad (13)$$

$$H(\sigma) = - \sum_{i}^{N} h_i \sigma_i - \sum_{i}^{N} \sum_{j<i}^{N} J_{ij} \sigma_i \sigma_j. \qquad (14)$$

where $\{h_i\}$ and $\{J_{ij}\}$ are the Lagrangian multipliers, $H(\sigma)$ is the Hamiltonian or energy function and $Z$ is the partition function, whose sum runs over all $2^N$ possible configurations of spin states $\sigma$. This prompts us to interpret this maximum entropy description as a mapping from the temporal averages of the neuronal network to that of a fully-connected spin lattice. Under this conceptual view, $h_i$ is analogous to a local external field acting on spin $i$, whereas $J_{ij}$ is an interaction constant between spins $i$ and $j$. The task is now to find the values of $h_i$ and $J_{ij}$ that reproduce the measured expectation values from the IF model. This is a particular example of the inverse Ising problem¹⁴, also known as Boltzmann Machine (BM) learning¹⁴, which consists in inferring the Hamiltonian of a certain complex multi-component system from its observed statistics. In principle, each parameter $h_i$ and $J_{ij}$ can be determined from the derivative of the logarithm of the partition function (13), $\langle \sigma_i \rangle = \frac{\partial}{\partial h_i} \ln[Z]$ and $\langle \sigma_i \sigma_j \rangle = \frac{\partial}{\partial J_{ij}} \ln[Z]$. However, the number of terms in $Z$ grows exponentially with $N$, as $2^N$, and an analytical
some time $t$ when its potential $v_i$ surpasses a threshold $v_c = 1$. Activity will then propagate to all post-synaptic neurons $j$ with incoming connections from $i$ according to the following equations³⁹:

$$v_j(t+1) = v_j(t) \pm v_i(t) u_i(t) g_{ij}, \quad (9)$$

$$u_i(t+1) = u_i(t) \cdot (1 - \delta u), \quad (10)$$

$$v_i(t+1) = 0, \quad (11)$$

where + and – stands for excitatory and inhibitory pre-synaptic neuron, respectively, and $\delta u = 0.05$ represents the fractional amount of neurotransmitter released at each neuronal firing⁶⁹. The timestep unit corresponds roughly to the joint interval of synaptic and axonal delay, i.e. the time interval between the generation of the action potential at the pre-synaptic neuron and the membrane potential change at the post-synaptic one, and is of the order of 10 milliseconds⁷⁰. We set a minimum value $v_{min} = -1$ for the membrane potential of each neuron, to prevent the possibility of a $v_i$ being systematically decreased to overly negative values when $p_{in} > 0$. After firing, a neuron enters into a refractory period of a single timestep $t_r = 1$ during which it is incapable of receiving or eliciting any activity. To keep the activity ongoing, we implement a small external stimulation. Namely, at every timestep, even during avalanches, a voltage input $\delta v = 0.1 v_c$ is added to a randomly chosen neuron. From Eq. (10) we see that the amount of synaptic resources gradually depletes as a neuron fires, eventually rendering it incapable of transmitting further signals. In real systems, this is counteracted by a recovery mechanism, where the readily releasable pool is slowly recharged over the span of seconds⁷¹. We assume a separation of timescales⁷², and recharge simultaneously the $u_i$ of all neurons by a certain amount $\delta u_{rec}$ only at the end of every avalanche, $u_i(t) \rightarrow u_i(t) + \delta u_{rec}$. By changing the value of $\delta u_{rec}$ one can adjust the dynamical state of the system³⁹,⁵⁷. Before performing any measurements, we let the dynamics evolve for a certain number of timesteps in order to shape the distribution of strengths $g_{ij}$ by LTP. We start by setting them uniformly distributed in the interval $g_{ij} \in [0.04, 0.06]$. According to the rules of Hebbian plasticity⁷³,⁷⁴, if some neuron $i$ frequently stimulates another neuron $j$, then the synapse from $i$ to $j$ will be strengthened. In this case, we increase the strength of the synapses $g_{ij}$ proportionally to the voltage variation induced in the post-synaptic neuron $j$ due to $i$ as $g_{ij}(t+1) = g_{ij}(t) + \delta g_j$, where $\delta g_j = \beta |v_j(t+1) - v_i(t)|$ and $\beta = 0.04$ sets the rate of this adaptation. On the other hand, synapses that are rarely active tend to weaken over time⁷⁴. Therefore, at the end of each avalanche, we decrease all terms $g_{ij}$ by the average increase in strength per synapse, $g_{ij}(t+1) = g_{ij}(t) - \frac{1}{N_s} \sum \delta g_j$ where $N_s$ is the number of synapses. In real networks, synapses that weaken consistently are eventually pruned. To avoid modifying the scale-free structure of the network, We modify the strengths $g_{ij}$ either for a fixed number $N_{aval} = 10^4$ of avalanches or until a strength $g_{ij}$ first reaches a minimum value $g_{min} = 10^{-5}$, where we then set that strength to $g_{ij} = g_{min}$.

## Maximum entropy modelling

The binarization of the IF model firing dynamics introduces the notion of the probability $P_{IF}(\sigma)$ to observe, during a time bin of duration $\Delta t_b$, any of the $2^N$ possible patterns of firing states $\sigma = \{\sigma_1, \sigma_2, ..., \sigma_N\}$ in a network of size $N$, with each $\sigma_i \in \{-1, 1\}$. We are interested in defining this $P_{IF}(\sigma)$ in a way that is consistent with the expectation values of the average local activities $\langle \sigma_i \rangle$ and two-point activities $\langle \sigma_i \sigma_j \rangle$ measured in the IF model for a given network. This surmounts to finding a probability distribution $P(\sigma)$ that maximizes the entropy¹⁴ $S = - \sum_{\sigma} P(\sigma) \ln[P(\sigma)]$, where $\sum_{\sigma}$ indicates a sum over all possible outcomes of $\sigma$, while subject to the constraints $\langle \sigma_i \rangle = \sum_{\sigma} \sigma_i P(\sigma)$ and $\langle \sigma_i \sigma_j \rangle = \sum_{\sigma} \sigma_i \sigma_j P(\sigma)$. Solving this problem using the method of Lagrangian multipliers¹⁵ (see Supplementary Information for the derivation) yields the distribution of spin states of a generalized Ising model at unit temperature⁷,

$$P(\sigma) = \frac{1}{Z} e^{-H(\sigma)}, \quad (12)$$

$$Z = \sum_{\sigma} e^{-H(\sigma)}, \quad (13)$$

$$H(\sigma) = - \sum_{i}^{N} h_i \sigma_i - \sum_{i}^{N} \sum_{j<i}^{N} J_{ij} \sigma_i \sigma_j. \quad (14)$$

where $\{h_i\}$ and $\{J_{ij}\}$ are the Lagrangian multipliers, $H(\sigma)$ is the Hamiltonian or energy function and $Z$ is the partition function, whose sum runs over all $2^N$ possible configurations of spin states $\sigma$. This prompts us to interpret this maximum entropy description as a mapping from the temporal averages of the neuronal network to that of a fully-connected spin lattice. Under this conceptual view, $h_i$ is analogous to a local external field acting on spin $i$, whereas $J_{ij}$ is an interaction constant between spins $i$ and $j$. The task is now to find the values of $h_i$ and $J_{ij}$ that reproduce the measured expectation values from the IF model. This is a particular example of the inverse Ising problem¹⁴, also known as Boltzmann Machine (BM) learning¹⁴, which consists in inferring the Hamiltonian of a certain complex multi-component system from its observed statistics. In principle, each parameter $h_i$ and $J_{ij}$ can be determined from the derivative of the logarithm of the partition function (13), $\langle \sigma_i \rangle = \frac{\partial}{\partial h_i} \ln[Z]$ and $\langle \sigma_i \sigma_j \rangle = \frac{\partial}{\partial J_{ij}} \ln[Z]$. However, the number of terms in $Z$ grows exponentially with $N$, as $2^N$, and an analytical
approach becomes intractable when $N \gtrsim 20$. To proceed, notice that we are trying to describe an empirical distribution $P_{IF}(\sigma)$, as observed from the measured temporal averages of the IF model, using the analytical description $P(\sigma)$, as given by Eq. (12). We want to choose the {$h_i$} and {$J_{ij}$} that minimize the loss in information when using $P(\sigma)$ as a proxy for $P_{IF}(\sigma)$. Specifically, this means we want the set of {$h_i$} and {$J_{ij}$} that minimize the so-called Kullback-Leibler divergence¹⁴ between these distributions,

$$D_{KL}(P_{IF}(\sigma) || P(\sigma)) = \sum_{\sigma} P_{IF}(\sigma) \ln \left[ \frac{P_{IF}(\sigma)}{P(\sigma)} \right]. \quad (15)$$

From the minimum condition equations one can show that $\frac{\partial D_{KL}}{\partial h_i} = 0 \Rightarrow \langle \sigma_i \rangle^{(IF)} = \langle \sigma_i \rangle^{(BM)}$ and $\frac{\partial D_{KL}}{\partial J_{ij}} = 0 \Rightarrow \langle \sigma_i \sigma_j \rangle^{(IF)} = \langle \sigma_i \sigma_j \rangle^{(BM)}$, where $\langle \sigma_i \rangle^{(IF)} = \sum_{\sigma} \sigma_i P_{IF}(\sigma)$ are the empirical averages measured in the IF model and $\langle \sigma_i \rangle^{(BM)} = \sum_{\sigma} \sigma_i P(\sigma)$ are the ones estimated from the Ising distribution (12), and analogously for the $\langle \sigma_i \sigma_j \rangle$. Therefore, a suitable method to search for the fields $h_i$ and interaction constants $J_{ij}$ of the Hamiltonian (14) is the following iterative scheme⁷,

$$h_i(x+1) = h_i(x) - \eta_h(x) \cdot (\langle \sigma_i \rangle^{(BM)} - \langle \sigma_i \rangle^{(IF)}), \quad (16)$$

$$J_{ij}(x+1) = J_{ij}(x) - \eta_J(x) \cdot (\langle \sigma_i \sigma_j \rangle^{(BM)} - \langle \sigma_i \sigma_j \rangle^{(IF)}), \quad (17)$$

where $\eta_h(x) = 2\eta_J(x) \propto x^{-\alpha}$ are decreasing learning rates, with $\alpha = 0.4$ for $N \le 40$, $\alpha = 0.6$ if $40 < N \le 120$, and $\alpha = 1.0$ otherwise. We set a smaller learning rate for the {$J_{ij}$} since their number ($\sim N^2$) is much larger when compared to the {$h_i$} ($N$), so we update their values at a slower rate to avoid divergences during the learning procedure. At each iteration $x$, the data sets $\{\langle \sigma_i \rangle^{(IF)}\}$ and $\{\langle \sigma_i \sigma_j \rangle^{(IF)}\}$ generated by the IF model are compared to those estimated by sampling the distribution (12), $\{\langle \sigma_i \rangle^{(BM)}\}$ and $\{\langle \sigma_i \sigma_j \rangle^{(BM)}\}$, using the sets of fields {$h_i(x)$} and coupling constants {$J_{ij}(x)$}. We start with $h_i(x = 1) = \langle \sigma_i \rangle^{(IF)}$ and $J_{ij}(x = 1) = 0$ and then iterate Eqs. (16) and (17) typically until $x = N_{BM} \sim 60000$. At each iteration, $\{\langle \sigma_i \rangle^{(BM)}\}$ and $\{\langle \sigma_i \sigma_j \rangle^{(BM)}\}$ are estimated from Monte Carlo simulations, using the Metropolis algorithm⁵⁵, by averaging over $M_c = 3 \cdot 10^5$ spin configurations. We disregard the first $150N$ configurations for systems with $N \le 120$, or up to $10^5 N^2$ configurations for the $N = 180$ case, in order to reduce correlations with the initial state, and use only every $2N$-th configuration for averaging, to reduce autocorrelations. At the end of the learning routine, we study the generalized Ising models with the set of fitted parameters {$h_i$} and {$J_{ij}$} by sampling the distribution (12), averaging over an increased amount of spin configurations $M_c = 3 \cdot 10^6$ for systems with $N \le 120$, or up to $M_c = 3 \cdot 10^7$ configurations for the system with $N = 180$, to reduce error bars.

## Data availability

The datasets generated and analyzed during the current study are available from the corresponding author on reasonable request.

Received: 29 January 2024; Accepted: 18 April 2024

Published online: 25 April 2024

## References

1. Hodgkin, A. L. & Huxley, A. F. A quantitative description of membrane current and its application to conduction and excitation in nerve. *J. Physiol.* **117**, 500. https://doi.org/10.1113/jphysiol.1952.sp004764 (1952).

2. Dayan, P. & Abbott, L. F. *Theoretical Neuroscience: Computational and Mathematical Modeling of Neural Systems* (The MIT Press, 2005).

3. Little, W. A. The existence of persistent states in the brain. *Math. Biosci.* **19**, 101–120. https://doi.org/10.1016/0025-5564(74)90031-5 (1974).

4. Hopfield, J. J. Neural networks and physical systems with emergent collective computational abilities. *Proc. Natl. Acad. Sci.* **79**, 2554–2558. https://doi.org/10.1073/pnas.79.8.2554 (1982).

5. Amit, D. J., Gutfreund, H. & Sompolinsky, H. Spin-glass models of neural networks. *Phys. Rev. A* **32**, 1007–1018. https://doi.org/10.1103/PhysRevA.32.1007 (1985).

6. Schneidman, E., Berry, M. J., Segev, R. & Bialek, W. Weak pairwise correlations imply strongly correlated network states in a neural population. *Nature* **440**, 1007–1012. https://doi.org/10.1038/nature04701 (2006).

7. Tkacik, G., Schneidman, E., Berry II, M. J. & Bialek, W. Spin glass models for a network of real neurons (2009). *arXiv:0912.5409 [q-bio]*.

8. Tkacik, G. et al. The simplest maximum entropy model for collective behavior in a neural network. *J. Stat. Mech. Theory Exp.* **2013**, P03011. https://doi.org/10.1088/1742-5468/2013/03/P03011 (2013).

9. Tkacik, G. et al. Searching for Collective Behavior in a Large Network of Sensory Neurons. *PLOS Comput. Biol.* **10**, e1003408. https://doi.org/10.1371/journal.pcbi.1003408 (2014).

10. Tkacik, G. et al. Thermodynamics and signatures of criticality in a network of neurons. *Proc. Natl. Acad. Sci.* **112**, 11508–11513. https://doi.org/10.1073/pnas.1514188112 (2015).

11. Meshulam, L., Gauthier, J. L., Brody, C. D., Tank, D. W. & Bialek, W. Collective Behavior of Place and Non-place Neurons in the Hippocampal Network. *Neuron* **96**, 1178–1191.e4. https://doi.org/10.1016/j.neuron.2017.10.027 (2017).

12. Chen, X., Randi, F., Leifer, A. M. & Bialek, W. Searching for collective behavior in a small brain. *Phys. Rev. E* **99**, 052418. https://doi.org/10.1103/PhysRevE.99.052418 (2019).

13. Meshulam, L., Gauthier, J. L., Brody, C. D., Tank, D. W. & Bialek, W. Successes and failures of simple statistical physics models for a network of real neurons, https://doi.org/10.48550/arXiv.2112.14735 (2023). https://doi.org/10.48550/arXiv.2112.14735 [physics, q-bio].
approach becomes intractable when $N \gtrsim 20$. To proceed, notice that we are trying to describe an empirical distribution $P_{IF}(\sigma)$, as observed from the measured temporal averages of the IF model, using the analytical description $P(\sigma)$, as given by Eq. (12). We want to choose the {$h_i$} and {$J_{ij}$} that minimize the loss in information when using $P(\sigma)$ as a proxy for $P_{IF}(\sigma)$. Specifically, this means we want the set of {$h_i$} and {$J_{ij}$} that minimize the so-called Kullback-Leibler divergence¹⁴ between these distributions,

$$D_{KL}(P_{IF}(\sigma) || P(\sigma)) = \sum_{\sigma} P_{IF}(\sigma) \ln \left[ \frac{P_{IF}(\sigma)}{P(\sigma)} \right]. \qquad (15)$$

From the minimum condition equations one can show that $\frac{\partial D_{KL}}{\partial h_i} = 0 \Rightarrow \langle \sigma_i \rangle^{(IF)} = \langle \sigma_i \rangle^{(BM)}$ and $\frac{\partial D_{KL}}{\partial J_{ij}} = 0 \Rightarrow \langle \sigma_i \sigma_j \rangle^{(IF)} = \langle \sigma_i \sigma_j \rangle^{(BM)}$, where $\langle \sigma_i \rangle^{(IF)} = \sum_{\sigma} \sigma_i P_{IF}(\sigma)$ are the empirical averages measured in the IF model and $\langle \sigma_i \rangle^{(BM)} = \sum_{\sigma} \sigma_i P(\sigma)$ are the ones estimated from the Ising distribution (12), and analogously for the $\langle \sigma_i \sigma_j \rangle$. Therefore, a suitable method to search for the fields $h_i$ and interaction constants $J_{ij}$ of the Hamiltonian (14) is the following iterative scheme⁷,

$$h_i(x+1) = h_i(x) - \eta_h(x) \cdot (\langle \sigma_i \rangle^{(BM)} - \langle \sigma_i \rangle^{(IF)}), \qquad (16)$$

$$J_{ij}(x+1) = J_{ij}(x) - \eta_J(x) \cdot (\langle \sigma_i \sigma_j \rangle^{(BM)} - \langle \sigma_i \sigma_j \rangle^{(IF)}), \qquad (17)$$

where $\eta_h(x) = 2\eta_J(x) \propto x^{-\alpha}$ are decreasing learning rates, with $\alpha = 0.4$ for $N \le 40$, $\alpha = 0.6$ if $40 < N \le 120$, and $\alpha = 1.0$ otherwise. We set a smaller learning rate for the {$J_{ij}$} since their number ($\sim N^2$) is much larger when compared to the {$h_i$} ($N$), so we update their values at a slower rate to avoid divergences during the learning procedure. At each iteration $x$, the data sets {$\langle \sigma_i \rangle^{(IF)}$} and {$\langle \sigma_i \sigma_j \rangle^{(IF)}$} generated by the IF model are compared to those estimated by sampling the distribution (12), {$\langle \sigma_i \rangle^{(BM)}$} and {$\langle \sigma_i \sigma_j \rangle^{(BM)}$}, using the sets of fields {$h_i(x)$} and coupling constants {$J_{ij}(x)$}. We start with $h_i(x = 1) = \langle \sigma_i \rangle^{(IF)}$ and $J_{ij}(x = 1) = 0$ and then iterate Eqs. (16) and (17) typically until $x = N_{BM} \sim 60000$. At each iteration, {$\langle \sigma_i \rangle^{(BM)}$} and {$\langle \sigma_i \sigma_j \rangle^{(BM)}$} are estimated from Monte Carlo simulations, using the Metropolis algorithm⁵⁵, by averaging over $M_c = 3 \cdot 10^5$ spin configurations. We disregard the first $150N$ configurations for systems with $N \le 120$, or up to $10^5 N^2$ configurations for the $N = 180$ case, in order to reduce correlations with the initial state, and use only every $2N$-th configuration for averaging, to reduce autocorrelations. At the end of the learning routine, we study the generalized Ising models with the set of fitted parameters {$h_i$} and {$J_{ij}$} by sampling the distribution (12), averaging over an increased amount of spin configurations $M_c = 3 \cdot 10^6$ for systems with $N \le 120$, or up to $M_c = 3 \cdot 10^7$ configurations for the system with $N = 180$, to reduce error bars.

## Data availability

The datasets generated and analyzed during the current study are available from the corresponding author on reasonable request.

Received: 29 January 2024; Accepted: 18 April 2024

Published online: 25 April 2024

## References

1. Hodgkin, A. L. & Huxley, A. F. A quantitative description of membrane current and its application to conduction and excitation in nerve. *J. Physiol.* **117**, 500. https://doi.org/10.1113/jphysiol.1952.sp004764 (1952).

2. Dayan, P. & Abbott, L. F. *Theoretical Neuroscience: Computational and Mathematical Modeling of Neural Systems* (The MIT Press, 2005).

3. Little, W. A. The existence of persistent states in the brain. *Math. Biosci.* **19**, 101–120. https://doi.org/10.1016/0025-5564(74)90031-5 (1974).

4. Hopfield, J. J. Neural networks and physical systems with emergent collective computational abilities. *Proc. Natl. Acad. Sci.* **79**, 2554–2558. https://doi.org/10.1073/pnas.79.8.2554 (1982).

5. Amit, D. J., Gutfreund, H. & Sompolinsky, H. Spin-glass models of neural networks. *Phys. Rev. A* **32**, 1007–1018. https://doi.org/10.1103/PhysRevA.32.1007 (1985).

6. Schneidman, E., Berry, M. J., Segev, R. & Bialek, W. Weak pairwise correlations imply strongly correlated network states in a neural population. *Nature* **440**, 1007–1012. https://doi.org/10.1038/nature04701 (2006).

7. Tkacik, G., Schneidman, E., Berry II, M. J. & Bialek, W. Spin glass models for a network of real neurons (2009). arXiv:0912.5409 [q-bio].

8. Tkacik, G. et al. The simplest maximum entropy model for collective behavior in a neural network. *J. Stat. Mech. Theory Exp.* **2013**, P03011. https://doi.org/10.1088/1742-5468/2013/03/P03011 (2013).

9. Tkacik, G. et al. Searching for Collective Behavior in a Large Network of Sensory Neurons. *PLOS Comput. Biol.* **10**, e1003408. https://doi.org/10.1371/journal.pcbi.1003408 (2014).

10. Tkacik, G. et al. Thermodynamics and signatures of criticality in a network of neurons. *Proc. Natl. Acad. Sci.* **112**, 11508–11513. https://doi.org/10.1073/pnas.1514188112 (2015).

11. Meshulam, L., Gauthier, J. L., Brody, C. D., Tank, D. W. & Bialek, W. Collective Behavior of Place and Non-place Neurons in the Hippocampal Network. *Neuron* **96**, 1178–1191.e4. https://doi.org/10.1016/j.neuron.2017.10.027 (2017).

12. Chen, X., Randi, F., Leifer, A. M. & Bialek, W. Searching for collective behavior in a small brain. *Phys. Rev. E* **99**, 052418. https://doi.org/10.1103/PhysRevE.99.052418 (2019).

13. Meshulam, L., Gauthier, J. L., Brody, C. D., Tank, D. W. & Bialek, W. Successes and failures of simple statistical physics models for a network of real neurons, https://doi.org/10.48550/arXiv.2112.14735 (2023). arXiv:2112.14735 [physics, q-bio].
14. Nguyen, H. C., Zecchina, R. & Berg, J. Inverse statistical problems: from the inverse Ising problem to data science. *Adv. Phys.* **66**, 197–261. https://doi.org/10.1080/00018732.2017.1341604 (2017).

15. Jaynes, E. T. Information Theory and Statistical Mechanics. *Phys. Rev.* **106**, 62–79. https://doi.org/10.1103/PhysRev.106.620 (1957).

16. Wainwright, M. J. & Jordan, M. I. Graphical Models, Exponential Families, and Variational Inference. *Found. Trends Mach. Learn.* **1**, 1–305. https://doi.org/10.1561/2200000001 (2008).

17. Phillips, S. J., Anderson, R. P. & Schapire, R. E. Maximum entropy modeling of species geographic distributions. *Ecol. Model.* **190**, 231–259. https://doi.org/10.1016/j.ecolmodel.2005.03.026 (2006).

18. Morcos, F. et al. Direct-coupling analysis of residue coevolution captures native contacts across many protein families. *Proc. Natl. Acad. Sci.* **108**, E1293–E1301. https://doi.org/10.1073/pnas.1111471108 (2011).

19. Figliuzzi, M., Barrat-Charlaix, P. & Weigt, M. How Pairwise Coevolutionary Models Capture the Collective Residue Variability in Proteins? *Mol. Biol. Evol.* **35**, 1018–1027. https://doi.org/10.1093/molbev/msy007 (2018).

20. Hopf, T. A. et al. Mutation effects predicted from sequence co-variation. *Nat. Biotechnol.* **35**, 128–135. https://doi.org/10.1038/nbt.3769 (2017).

21. Bialek, W. et al. Social interactions dominate speed control in poising natural flocks near criticality. *Proc. Natl. Acad. Sci.* **111**, 7212–7217. https://doi.org/10.1073/pnas.1324045111 (2014).

22. Lai, J. et al. Maximum Entropy Analysis of Bird Diversity and Environmental Variables in Nanjing Megapolis, China. *Sustainability* **16**, 2139. https://doi.org/10.3390/su16052139 (2024).

23. Burleson-Lesser, K., Morone, F., DeGuzman, P., Parra, L. C. & Makse, H. A. Collective Behaviour in Video Viewing: A Thermo-dynamic Analysis of Gaze Position. *PLOS ONE* **12**, e0168995. https://doi.org/10.1371/journal.pone.0168995 (2017).

24. Torres, D. et al. Eye-tracking as a proxy for coherence and complexity of texts. *PLOS ONE* **16**, e0260236. https://doi.org/10.1371/journal.pone.0260236 (2021).

25. Cabrera, J. S. & Lee, H. S. Flood risk assessment for Davao Oriental in the Philippines using geographic information system-based multi-criteria analysis and the maximum entropy model. *J. Flood Risk Manag.* **13**, e12607. https://doi.org/10.1111/jfr3.12607 (2020).

26. Daniotti, S., Monechi, B. & Ubaldi, E. A maximum entropy approach for the modelling of car-sharing parking dynamics. *Sci. Rep.* **13**, 2993. https://doi.org/10.1038/s41598-023-30134-9 (2023).

27. De Paola, P. Real Estate Valuations with Small Dataset: A Novel Method Based on the Maximum Entropy Principle and Lagrange Multipliers. *Real Estate* **1**, 26–40. https://doi.org/10.3390/realestate1010003 (2024).

28. Ganmor, E., Segev, R. & Schneidman, E. The architecture of functional interaction networks in the retina. *J. Neurosci.* **31**, 3044–3054. https://doi.org/10.1523/JNEUROSCI.3682-10.2011 (2011).

29. Humplik, J. & Tkačik, G. Probabilistic models for neural populations that naturally capture global coupling and criticality. *PLOS Comput. Biol.* **13**, e1005763. https://doi.org/10.1371/journal.pcbi.1005763 (2017).

30. Delamare, G. & Ferrari, U. Time-dependent maximum entropy model for populations of retinal ganglion cells. https://doi.org/10.1101/2022.07.13.498395 (2022).

31. Nonnenmacher, M., Behrens, C., Berens, P., Bethge, M. & Macke, J. H. Signatures of criticality arise from random subsampling in simple population models. *PLOS Comput. Biol.* **13**, e1005718. https://doi.org/10.1371/journal.pcbi.1005718 (2017).

32. Hahn, G. et al. Spontaneous cortical activity is transiently poised close to criticality. *PLoS Comput. Biol.* **13**, e1005543. https://doi.org/10.1371/journal.pcbi.1005543 (2017).

33. Sampaio Filho, C. I. et al. Ising-like model replicating time-averaged spiking behaviour of in vivo and in vitro neuronal networks (2023). Preprint.

34. White, J. G., Southgate, E., Thomson, J. N. & Brenner, S. The structure of the nervous system of the nematode Caenorhabditis elegans. *Philos. Trans. R. Soc. Lond. B Biol. Sci.* **314**, 1–340. https://doi.org/10.1098/rstb.1986.0056 (1986).

35. Gardella, C., Marre, O. & Mora, T. Modeling the correlated activity of neural populations: A review. *Neural Comput.* **31**, 233–269. https://doi.org/10.1162/neco_a_01154 (2019).

36. Buccino, A. P., Garcia, S. & Yger, P. Spike sorting: new trends and challenges of the era of high-density probes. *Prog. Biomed. Eng.* **4**, 022005. https://doi.org/10.1088/2516-1091/ac6b96 (2022).

37. Ventura, V. & Gerkin, R. C. Accurately estimating neuronal correlation requires a new spike-sorting paradigm. *Proc. Natl. Acade. Sci.* **109**, 7230–7235. https://doi.org/10.1073/pnas.1115236109 (2012).

38. Pastore, V. P., Massobrio, P., Godjoski, A. & Martinioia, S. Identification of excitatory-inhibitory links and network topology in large-scale neuronal assemblies from multi-electrode recordings. *PLOS Comput. Biol.* **14**, e1006381. https://doi.org/10.1371/journal.pcbi.1006381 (2018).

39. Michiels van Kessenich, L., Luković, M., de Arcangelis, L. & Herrmann, H. J. Critical neural networks with short- and long-term plasticity. *Phys. Rev. E* **97**, 032312. https://doi.org/10.1103/PhysRevE.97.032312 (2018).

40. Beggs, J. M. & Plenz, D. Neuronal Avalanches in Neocortical Circuits. *J. Neurosci.* **23**, 11167–11177. https://doi.org/10.1523/JNEUR.OSCI.23-35-11167.2003 (2003).

41. de Arcangelis, L., Perrone-Capano, C. & Herrmann, H. J. Self-organized criticality model for brain plasticity. *Phys. Rev. Lett.* **96**, 028107. https://doi.org/10.1103/PhysRevLett.96.028107 (2006).

42. Levina, A., Herrmann, J. M. & Geisel, T. Dynamical synapses causing self-organized criticality in neural networks. *Nat. Phys.* **3**, 857–860. https://doi.org/10.1038/nphys758 (2007).

43. Shew, W. L., Yang, H., Yu, S., Roy, R. & Plenz, D. Information capacity and transmission are maximized in balanced cortical networks with neuronal avalanches. *J. Neurosci.* **31**, 55–63. https://doi.org/10.1523/JNEUROSCI.4637-10.2011 (2011).

44. Haimovici, A., Tagliazucchi, E., Balenzuela, P. & Chialvo, D. R. Brain organization into resting state networks emerges at criticality on a model of the human connectome. *Phys. Rev. Lett.* **110**, 178101. https://doi.org/10.1103/PhysRevLett.110.178101 (2013).

45. Priesemann, V., Valderrama, M., Wibral, M. & Quyen, M. L. V. Neuronal avalanches differ from wakefulness to deep sleep–evidence from intracranial depth recordings in humans. *PLOS Comput. Biol.* **9**, e1002985. https://doi.org/10.1371/journal.pcbi.1002985 (2013).

46. Shriki, O. et al. Neuronal avalanches in the resting MEG of the human brain. *J. Neurosci.* **33**, 7079–7090. https://doi.org/10.1523/JNEUROSCI.4286-12.2013 (2013).

47. Petermann, T. et al. Spontaneous cortical activity in awake monkeys composed of neuronal avalanches. *Proc. Natl. Acad. Sci.* **106**, 15921–15926. https://doi.org/10.1073/pnas.0904089106 (2009).

48. Mora, T., Deny, S. & Marre, O. Dynamical criticality in the collective activity of a population of retinal neurons. *Phys. Rev. Lett.* **114**, 078105. https://doi.org/10.1103/PhysRevLett.114.078105 (2015).

49. Shew, W. L. et al. Adaptation to sensory input tunes visual cortex to criticality. *Nat. Phys.* **11**, 659–663. https://doi.org/10.1038/nphys3370 (2015).

50. Ponce-Alvarez, A., Jouary, A., Privat, M., Deco, G. & Sumbre, G. Whole-brain neuronal activity displays crackling noise dynamics. *Neuron* **100**, 1446–1459. https://doi.org/10.1016/j.neuron.2018.10.045 (2018).

51. Michiels van Kessenich, L., de Arcangelis, L. & Herrmann, H. J. Synaptic plasticity and neuronal refractory time cause scaling behaviour of neuronal avalanches. *Sci. Rep.* **6**, 32071. https://doi.org/10.1038/srep32071 (2016).

52. Zapperi, S., Lauritsen, K. B. & Stanley, H. E. Self-organized branching processes: Mean-field theory for avalanches. *Phys. Rev. Lett.* **75**, 4071–4074. https://doi.org/10.1103/PhysRevLett.75.4071 (1995).

53. Sherrington, D. & Kirkpatrick, S. Solvable model of a spin-glass. *Phys. Rev. Lett.* **35**, 1792–1796. https://doi.org/10.1103/PhysRevLett.35.1792 (1975).
14. Nguyen, H. C., Zecchina, R. & Berg, J. Inverse statistical problems: from the inverse Ising problem to data science. *Adv. Phys.* **66**, 197–261. https://doi.org/10.1080/00018732.2017.1341604 (2017).

15. Jaynes, E. T. Information Theory and Statistical Mechanics. *Phys. Rev.* **106**, 62–79. https://doi.org/10.1103/PhysRev.106.620 (1957).

16. Wainwright, M. J. & Jordan, M. I. Graphical Models, Exponential Families, and Variational Inference. *Found. Trends Mach. Learn.* **1**, 1–305. https://doi.org/10.1561/2200000001 (2008).

17. Phillips, S. J., Anderson, R. P. & Schapire, R. E. Maximum entropy modeling of species geographic distributions. *Ecol. Model.* **190**, 231–259. https://doi.org/10.1016/j.ecolmodel.2005.03.026 (2006).

18. Morcos, F. et al. Direct-coupling analysis of residue coevolution captures native contacts across many protein families. *Proc. Natl. Acad. Sci.* **108**, E1293–E1301. https://doi.org/10.1073/pnas.1111471108 (2011).

19. Figliuzzi, M., Barrat-Charlaix, P. & Weigt, M. How Pairwise Coevolutionary Models Capture the Collective Residue Variability in Proteins? *Mol. Biol. Evol.* **35**, 1018–1027. https://doi.org/10.1093/molbev/msy007 (2018).

20. Hopf, T. A. et al. Mutation effects predicted from sequence co-variation. *Nat. Biotechnol.* **35**, 128–135. https://doi.org/10.1038/nbt.3769 (2017).

21. Bialek, W. et al. Social interactions dominate speed control in poising natural flocks near criticality. *Proc. Natl. Acad. Sci.* **111**, 7212–7217. https://doi.org/10.1073/pnas.1324045111 (2014).

22. Lai, J. et al. Maximum Entropy Analysis of Bird Diversity and Environmental Variables in Nanjing Megapolis, China. *Sustainability* **16**, 2139. https://doi.org/10.3390/su16052139 (2024).

23. Burleson-Lesser, K., Morone, F., DeGuzman, P., Parra, L. C. & Makse, H. A. Collective Behaviour in Video Viewing: A Thermo-dynamic Analysis of Gaze Position. *PLOS ONE* **12**, e0168995. https://doi.org/10.1371/journal.pone.0168995 (2017).

24. Torres, D. et al. Eye-tracking as a proxy for coherence and complexity of texts. *PLOS ONE* **16**, e0260236. https://doi.org/10.1371/journal.pone.0260236 (2021).

25. Cabrera, J. S. & Lee, H. S. Flood risk assessment for Davao Oriental in the Philippines using geographic information system-based multi-criteria analysis and the maximum entropy model. *J. Flood Risk Manag.* **13**, e12607. https://doi.org/10.1111/jfr3.12607 (2020).

26. Daniotti, S., Monechi, B. & Ubaldi, E. A maximum entropy approach for the modelling of car-sharing parking dynamics. *Sci. Rep.* **13**, 2993. https://doi.org/10.1038/s41598-023-30134-9 (2023).

27. De Paola, P. Real Estate Valuations with Small Dataset: A Novel Method Based on the Maximum Entropy Principle and Lagrange Multipliers. *Real Estate* **1**, 26–40. https://doi.org/10.3390/realestate1010003 (2024).

28. Ganmor, E., Segev, R. & Schneidman, E. The architecture of functional interaction networks in the retina. *J. Neurosci.* **31**, 3044–3054. https://doi.org/10.1523/JNEUROSCI.3682-10.2011 (2011).

29. Humplik, J. & Tkačik, G. Probabilistic models for neural populations that naturally capture global coupling and criticality. *PLOS Comput. Biol.* **13**, e1005763. https://doi.org/10.1371/journal.pcbi.1005763 (2017).

30. Delamare, G. & Ferrari, U. Time-dependent maximum entropy model for populations of retinal ganglion cells. https://doi.org/10.1101/2022.07.13.498395 (2022).

31. Nonnenmacher, M., Behrens, C., Berens, P., Bethge, M. & Macke, J. H. Signatures of criticality arise from random subsampling in simple population models. *PLOS Comput. Biol.* **13**, e1005718. https://doi.org/10.1371/journal.pcbi.1005718 (2017).

32. Hahn, G. et al. Spontaneous cortical activity is transiently poised close to criticality. *PLoS Comput. Biol.* **13**, e1005543. https://doi.org/10.1371/journal.pcbi.1005543 (2017).

33. Sampaio Filho, C. I. et al. Ising-like model replicating time-averaged spiking behaviour of in vivo and in vitro neuronal networks (2023). Preprint.

34. White, J. G., Southgate, E., Thomson, J. N. & Brenner, S. The structure of the nervous system of the nematode Caenorhabditis elegans. *Philos. Trans. R. Soc. Lond. B Biol. Sci.* **314**, 1–340. https://doi.org/10.1098/rstb.1986.0056 (1986).

35. Gardella, C., Marre, O. & Mora, T. Modeling the correlated activity of neural populations: A review. *Neural Comput.* **31**, 233–269. https://doi.org/10.1162/neco_a_01154 (2019).

36. Buccino, A. P., Garcia, S. & Yger, P. Spike sorting: new trends and challenges of the era of high-density probes. *Prog. Biomed. Eng.* **4**, 022005. https://doi.org/10.1088/2516-1091/ac6b96 (2022).

37. Ventura, V. & Gerkin, R. C. Accurately estimating neuronal correlation requires a new spike-sorting paradigm. *Proc. Natl. Acade. Sci.* **109**, 7230–7235. https://doi.org/10.1073/pnas.1115236109 (2012).

38. Pastore, V. P., Massobrio, P., Godjoski, A. & Martinioia, S. Identification of excitatory-inhibitory links and network topology in large-scale neuronal assemblies from multi-electrode recordings. *PLOS Comput. Biol.* **14**, e1006381. https://doi.org/10.1371/journal.pcbi.1006381 (2018).

39. Michiels van Kessenich, L., Luković, M., de Arcangelis, L. & Herrmann, H. J. Critical neural networks with short- and long-term plasticity. *Phys. Rev. E* **97**, 032312. https://doi.org/10.1103/PhysRevE.97.032312 (2018).

40. Beggs, J. M. & Plenz, D. Neuronal Avalanches in Neocortical Circuits. *J. Neurosci.* **23**, 11167–11177. https://doi.org/10.1523/JNEUR.OSCI.23-35-11167.2003 (2003).

41. de Arcangelis, L., Perrone-Capano, C. & Herrmann, H. J. Self-organized criticality model for brain plasticity. *Phys. Rev. Lett.* **96**, 028107. https://doi.org/10.1103/PhysRevLett.96.028107 (2006).

42. Levina, A., Herrmann, J. M. & Geisel, T. Dynamical synapses causing self-organized criticality in neural networks. *Nat. Phys.* **3**, 857–860. https://doi.org/10.1038/nphys758 (2007).

43. Shew, W. L., Yang, H., Yu, S., Roy, R. & Plenz, D. Information capacity and transmission are maximized in balanced cortical networks with neuronal avalanches. *J. Neurosci.* **31**, 55–63. https://doi.org/10.1523/JNEUROSCI.4637-10.2011 (2011).

44. Haimovici, A., Tagliazucchi, E., Balenzuela, P. & Chialvo, D. R. Brain organization into resting state networks emerges at criticality on a model of the human connectome. *Phys. Rev. Lett.* **110**, 178101. https://doi.org/10.1103/PhysRevLett.110.178101 (2013).

45. Priesemann, V., Valderrama, M., Wibral, M. & Quyen, M. L. V. Neuronal avalanches differ from wakefulness to deep sleep–evidence from intracranial depth recordings in humans. *PLOS Comput. Biol.* **9**, e1002985. https://doi.org/10.1371/journal.pcbi.1002985 (2013).

46. Shriki, O. et al. Neuronal avalanches in the resting MEG of the human brain. *J. Neurosci.* **33**, 7079–7090. https://doi.org/10.1523/JNEUROSCI.4286-12.2013 (2013).

47. Petermann, T. et al. Spontaneous cortical activity in awake monkeys composed of neuronal avalanches. *Proc. Natl. Acad. Sci.* **106**, 15921–15926. https://doi.org/10.1073/pnas.0904089106 (2009).

48. Mora, T., Deny, S. & Marre, O. Dynamical criticality in the collective activity of a population of retinal neurons. *Phys. Rev. Lett.* **114**, 078105. https://doi.org/10.1103/PhysRevLett.114.078105 (2015).

49. Shew, W. L. et al. Adaptation to sensory input tunes visual cortex to criticality. *Nat. Phys.* **11**, 659–663. https://doi.org/10.1038/nphys3370 (2015).

50. Ponce-Alvarez, A., Jouary, A., Privat, M., Deco, G. & Sumbre, G. Whole-brain neuronal activity displays crackling noise dynamics. *Neuron* **100**, 1446–1459. https://doi.org/10.1016/j.neuron.2018.10.045 (2018).

51. Michiels van Kessenich, L., de Arcangelis, L. & Herrmann, H. J. Synaptic plasticity and neuronal refractory time cause scaling behaviour of neuronal avalanches. *Sci. Rep.* **6**, 32071. https://doi.org/10.1038/srep32071 (2016).

52. Zapperi, S., Lauritsen, K. B. & Stanley, H. E. Self-organized branching processes: Mean-field theory for avalanches. *Phys. Rev. Lett.* **75**, 4071–4074. https://doi.org/10.1103/PhysRevLett.75.4071 (1995).

53. Sherrington, D. & Kirkpatrick, S. Solvable model of a spin-glass. *Phys. Rev. Lett.* **35**, 1792–1796. https://doi.org/10.1103/PhysRevLett.35.1792 (1975).
54. Mezard, M., Parisi, G. & Virasoro, M. *Spin Glass Theory and Beyond: An Introduction to the Replica Method and Its Applications*, vol. 9 of *World Scientific Lecture Notes in Physics* (World Scientific, 1986).

55. Böttcher, L. & Herrmann, H. J. *Computational Statistical Physics* (Cambridge University Press, Cambridge, 2021).

56. Raimo, D., Sarracino, A. & de Arcangelis, L. Role of inhibitory neurons in temporal correlations of critical and supercritical spontaneous activity. *Phys. A* **565**, 125555. https://doi.org/10.1016/j.physa.2020.125555 (2021).

57. 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. https://doi.org/10.1103/PhysRevE.106.024304 (2022).

58. Yeh, F.-C. et al. Maximum entropy approaches to living neural networks. *Entropy* **12**, 89–106. https://doi.org/10.3390/e12010089 (2010).

59. van der Plas, T. L. et al. Neural assemblies uncovered by generative modeling explain whole-brain activity statistics and reflect structural connectivity. *eLife* **12**, e83139. https://doi.org/10.7554/eLife.83139 (2023).

60. Maoz, O., Tkačik, G., Esteki, M. S., Kiani, R. & Schneidman, E. Learning probabilistic neural representations with randomly connected circuits. *Proc. Natl. Acad. Sci.* **117**, 25066–25073. https://doi.org/10.1073/pnas.1912804117 (2020).

61. Lynn, C. W., Yu, Q., Pang, R., Bialek, W. & Palmer, S. E. Exactly solvable statistical physics models for large neuronal populations, https://doi.org/10.48550/arXiv.2310.10860 (2023). arXiv:2310.10860 [cond-mat, physics:physics, q-bio].

62. Sarracino, A., Arviv, O., Shriki, O. & de Arcangelis, L. Predicting brain evoked response to external stimuli from temporal correlations of spontaneous activity. *Phys. Rev. Res.* **2**, 033355. https://doi.org/10.1103/PhysRevResearch.2.033355 (2020).

63. Nandi, M. K., de Candia, A., Sarracino, A., Herrmann, H. J. & de Arcangelis, L. Fluctuation-dissipation relations in the imbalanced Wilson-Cowan model. *Phys. Rev. E107*, 064307. https://doi.org/10.1103/PhysRevE.107.064307 (2023).

64. Berger, D., Varriale, E., van Kessenich, L. M., Herrmann, H. J. & de Arcangelis, L. Three cooperative mechanisms required for recovery after brain damage. *Sci. Rep.* **9**, 15858. https://doi.org/10.1038/s41598-019-50946-y (2019).

65. Eguiluz, V. M., Chialvo, D. R., Cecchi, G. A., Baliki, M. & Apkarian, A. V. Scale-free brain functional networks. *Phys. Rev. Lett.* **94**, 018102. https://doi.org/10.1103/PhysRevLett.94.018102 (2005).

66. Roerig, B. & Chen, B. Relationships of local inhibitory and excitatory circuits to orientation preference maps in ferret visual cortex. *Cereb. Cortex* **12**, 187–198. https://doi.org/10.1093/cercor/12.2.187 (2002).

67. Kaeser, P. S. & Regehr, W. G. The readily releasable pool of synaptic vesicles. *Curr. Opin. Neurobiol.* **43**, 63–70. https://doi.org/10.1016/j.conb.2016.12.012 (2017).

68. Abraham, W. C. How long will long-term potentiation last? *Philos. Trans. R. Soc. Lond. Ser. B Biol. Sci.* **358**, 735–744. https://doi.org/10.1098/rstb.2002.1222 (2003).

69. Ikeda, K. & Bekkers, J. M. Counting the number of releasable synaptic vesicles in a presynaptic terminal. *Proc. Natl. Acad. Sci.* **106**, 2945–2950. https://doi.org/10.1073/pnas.0811017106 (2009).

70. Lemaréchal, J.-D. et al. A brain atlas of axonal and synaptic delays based on modelling of cortico-cortical evoked potentials. *Brain* **145**, 1653–1667. https://doi.org/10.1093/brain/awab362 (2021).

71. Markram, H. & Tsodyks, M. Redistribution of synaptic efficacy between neocortical pyramidal neurons. *Nature* **382**, 807–810. https://doi.org/10.1038/382807a0 (1996).

72. Zeraati, R., Priesemann, V. & Levina, A. Self-Organization Toward Criticality by Synaptic Plasticity. *Front. Phys.* **9** (2021).

73. Hebb, D. O. *The Organization of Behavior* (John Wiley & Sons Inc, New York, 1949).

74. Bi, G.-q. & Poo, M.-m. Synaptic Modifications in Cultured Hippocampal Neurons: Dependence on Spike Timing, Synaptic Strength, and Postsynaptic Cell Type. *J. Neurosci.* **18**, 10464–10472. https://doi.org/10.1523/JNEUROSCI.18-24-10464.1998 (1998).

## Acknowledgements

L.d.A. acknowledges support from the Italian MUR project PRIN2017WZFTZP and from NEXTGENERATIONEU (NGEU) funded by the Ministry of University and Research (MUR), National Recovery and Resilience Plan (NRRP), and project MNESYS (PE0000006)-A multiscale integrated approach to the study of the nervous system in health and disease (DN. 1553 11.10.2022). H.J.H. thanks FUNCAP and INCT-SC for financial support. J.S.A.J. thanks the Brazilian agencies CNPq, CAPES and FUNCAP for financial support.

## Author contributions

T.S.A.N.S. conducted the numerical simulations. T.S.A.N.S. and C.I.N.S.F. developed the software. H.J.H., J.S.A.J. and L.d.A. conceived the research. T.S.A.N.S., H.J.H. and L.d.A. analysed the results and wrote the manuscript. All authors contributed to the review of the manuscript.

## 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-60117-3.

**Correspondence** and requests for materials should be addressed to T.S.A.N.S.

**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_25.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/.
54. Mezard, M., Parisi, G. & Virasoro, M. *Spin Glass Theory and Beyond: An Introduction to the Replica Method and Its Applications*, vol. 9 of *World Scientific Lecture Notes in Physics* (World Scientific, 1986).

55. Böttcher, L. & Herrmann, H. J. *Computational Statistical Physics* (Cambridge University Press, Cambridge, 2021).

56. Raimo, D., Sarracino, A. & de Arcangelis, L. Role of inhibitory neurons in temporal correlations of critical and supercritical spontaneous activity. *Phys. A* **565**, 125555. https://doi.org/10.1016/j.physa.2020.125555 (2021).

57. 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. https://doi.org/10.1103/PhysRevE.106.024304 (2022).

58. Yeh, F.-C. et al. Maximum entropy approaches to living neural networks. *Entropy* **12**, 89–106. https://doi.org/10.3390/e12010089 (2010).

59. van der Plas, T. L. et al. Neural assemblies uncovered by generative modeling explain whole-brain activity statistics and reflect structural connectivity. *eLife* **12**, e83139. https://doi.org/10.7554/eLife.83139 (2023).

60. Maoz, O., Tkačik, G., Esteki, M. S., Kiani, R. & Schneidman, E. Learning probabilistic neural representations with randomly connected circuits. *Proc. Natl. Acad. Sci.* **117**, 25066–25073. https://doi.org/10.1073/pnas.1912804117 (2020).

61. Lynn, C. W., Yu, Q., Pang, R., Bialek, W. & Palmer, S. E. Exactly solvable statistical physics models for large neuronal populations, https://doi.org/10.48550/arXiv.2310.10860 (2023). arXiv:2310.10860 [cond-mat, physics:physics, q-bio].

62. Sarracino, A., Arxiv, O., Shriki, O. & de Arcangelis, L. Predicting brain evoked response to external stimuli from temporal correlations of spontaneous activity. *Phys. Rev. Res.* **2**, 033355. https://doi.org/10.1103/PhysRevResearch.2.033355 (2020).

63. Nandi, M. K., de Candia, A., Sarracino, A., Herrmann, H. J. & de Arcangelis, L. Fluctuation-dissipation relations in the imbalanced Wilson-Cowan model. *Phys. Rev. E107*, 064307. https://doi.org/10.1103/PhysRevE.107.064307 (2023).

64. Berger, D., Varriale, E., van Kessenich, L. M., Herrmann, H. J. & de Arcangelis, L. Three cooperative mechanisms required for recovery after brain damage. *Sci. Rep.* **9**, 15858. https://doi.org/10.1038/s41598-019-50946-y (2019).

65. Eguiluz, V. M., Chialvo, D. R., Cecchi, G. A., Baliki, M. & Apkarian, A. V. Scale-free brain functional networks. *Phys. Rev. Lett.* **94**, 018102. https://doi.org/10.1103/PhysRevLett.94.018102 (2005).

66. Roerig, B. & Chen, B. Relationships of local inhibitory and excitatory circuits to orientation preference maps in ferret visual cortex. *Cereb. Cortex* **12**, 187–198. https://doi.org/10.1093/cercor/12.2.187 (2002).

67. Kaeser, P. S. & Regehr, W. G. The readily releasable pool of synaptic vesicles. *Curr. Opin. Neurobiol.* **43**, 63–70. https://doi.org/10.1016/j.conb.2016.12.012 (2017).

68. Abraham, W. C. How long will long-term potentiation last? *Philos. Trans. R. Soc. Lond. Ser. B Biol. Sci.* **358**, 735–744. https://doi.org/10.1098/rstb.2002.1222 (2003).

69. Ikeda, K. & Bekkers, J. M. Counting the number of releasable synaptic vesicles in a presynaptic terminal. *Proc. Natl. Acad. Sci.* **106**, 2945–2950. https://doi.org/10.1073/pnas.0811017106 (2009).

70. Lemaréchal, J.-D. et al. A brain atlas of axonal and synaptic delays based on modelling of cortico-cortical evoked potentials. *Brain* **145**, 1653–1667. https://doi.org/10.1093/brain/awab362 (2021).

71. Markram, H. & Tsodyks, M. Redistribution of synaptic efficacy between neocortical pyramidal neurons. *Nature* **382**, 807–810. https://doi.org/10.1038/382807a0 (1996).

72. Zeraati, R., Priesemann, V. & Levina, A. Self-Organization Toward Criticality by Synaptic Plasticity. *Front. Phys.* **9** (2021).

73. Hebb, D. O. *The Organization of Behavior* (John Wiley & Sons Inc, New York, 1949).

74. Bi, G.-q. & Poo, M.-m. Synaptic Modifications in Cultured Hippocampal Neurons: Dependence on Spike Timing, Synaptic Strength, and Postsynaptic Cell Type. *J. Neurosci.* **18**, 10464–10472. https://doi.org/10.1523/JNEUROSCI.18-24-10464.1998 (1998).

## Acknowledgements

L.d.A. acknowledges support from the Italian MUR project PRIN2017WZFTZP and from NEXTGENERATIONEU (NGEU) funded by the Ministry of University and Research (MUR), National Recovery and Resilience Plan (NRRP), and project MNESYS (PE0000006)-A multiscale integrated approach to the study of the nervous system in health and disease (DN. 1553 11.10.2022). H.J.H. thanks FUNCAP and INCT-SC for financial support. J.S.A.J. thanks the Brazilian agencies CNPq, CAPES and FUNCAP for financial support.

## Author contributions

T.S.A.N.S. conducted the numerical simulations. T.S.A.N.S. and C.I.N.S.F. developed the software. H.J.H., J.S.A.J. and L.d.A. conceived the research. T.S.A.N.S., H.J.H. and L.d.A. analysed the results and wrote the manuscript. All authors contributed to the review of the manuscript.

## 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-60117-3.

**Correspondence** and requests for materials should be addressed to T.S.A.N.S.

**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_25.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/.