Vito Dichio

Fabrizio De

Vico Fallani

Statistical models of complex brain networks: a maximum entropy approach

Keywords: Statistical modeling, Complex systems, Exponential random graph model, Brain networks, Inference, Maximum entropy principle . . . . . . . . . . 19 4.3 Predicting states in temporal brain networks . . . . . . . . . . . . . statistics . . . . . . . . 26 5.3 Explanation: ERGMs and Stat.Mech. . . . .

The brain is a highly complex system. Most of such complexity stems from the intermingled connections between its parts, which give rise to rich dynamics and to the emergence of high-level cognitive functions. Disentangling the underlying network structure is crucial to understand the brain functioning under both healthy and pathological conditions.

Yet, analyzing brain networks is challenging, in part because their structure represents only one possible realization of a generative stochastic process which is in general unknown. Having a formal way to cope with such intrinsic variability is therefore central for the characterization of brain network properties.

Addressing this issue entails the development of appropriate tools mostly adapted from network science and statistics. Here, we focus on a particular class of maximum entropy models for networks, i.e. exponential random graph models (ERGMs), as a parsimonious approach to identify the local connection mechanisms behind observed global network structure. Efforts are reviewed on the quest for basic organizational properties of human brain networks, as well as on the identification of predictive biomarkers of neurological diseases such as stroke.

We conclude with a discussion on how emerging results and tools from statistical graph modeling, associated with forthcoming improvements in experimental data acquisition, could lead to a finer probabilistic description of complex systems in network neuroscience.

Introduction

The human brain is a biological system of tremendous complexity. At different scales of neuronal organization, the paradigm of a system "made up of a large number of parts that interact in a nonsimple way" [START_REF] Herbert | The architecture of complexity[END_REF] turns out to be an apt abstraction. Notably, it suits neurons interacting through synapses at the microscale as well as brain regions' activity coordinating at the macroscale and resulting in the rich spectrum of mind's functional states. The study of the brain as a complex system has flourished in the last 20 years, drawing analogies from the physics of disordered systems, graph theory, dynamical systems and fueled by the modern deluge of data upon nearly all scientific fields [START_REF] Christopher | The physics of brain network structure, function and control[END_REF].

By explicitly representing the interactions between the system's components, networks, or graphs, constitute a natural and powerful way to inquire its organizing principles. The description of the brain under the lens of network science has lead to a number of fundamental results. Topological network properties, such as node centrality, modular organization and global efficiency, play a fundamental role in the emergence of basic physiological functions, as well as on the apperance of many brain diseases [START_REF] Bullmore | Complex brain networks: graph theoretical analysis of structural and functional systems[END_REF][START_REF] Bassett | Network neuroscience[END_REF][START_REF] Cornelis | Modern network science of neurological disorders[END_REF].

Beyond purely descriptive analyses, important questions on the network structure include what is the underlying generative process and how local wiring rules result in the observed large-scale properties [START_REF] Bassett | On the nature and use of models in network neuroscience[END_REF]. To this end, network models based on random edge rewiring or nodal preferential attachment rules have been initially investigated and allowed to assess nontrivial properties of brain networks, such as integration and segregation of information or the presence of few nodes with a high number of connections, i.e. the so-called hubs [START_REF] Sporns | Complex networks: small-world and scale-free architectures[END_REF]. In addition to purely topological models, brain networks can be more finely modeled by explicitly taking into account the spatial position of the nodes so as to penalize the cost of having long distance connections and minimize the associated metabolic consumption [START_REF] Richard F Betzel | Generative models of the human connectome[END_REF][START_REF] Bullmore | The economy of brain network organization[END_REF].

While these models allow to identify putative mechanisms involved in generating the observed brain networks, they implicitly make the assumption that the connections have a physical meaning, i.e. they are tangible quantities. However, it is important to remind that almost all brain networks are currently inferred from experimental data through advanced tools from image and signal processing (Fig. 1a) [START_REF] De | Graph analysis of functional brain networks: practical issues in translational neuroscience[END_REF]. In particular,

𝑃(𝐺)

𝐺 !"#$"$# 𝐺 "%&#!' 𝐺 * links are statistical estimates of possibly existing anatomical pathways or functional interactions between different brain areas and they are affected by some uncertainty (Fig. 1b).

Statistical modeling has been introduced in an effort to identify the probability distribution of all the possible network realizations associated with an observed one (Fig. 1c). Approaches of this kind include generating a-posteriori network surrogates preserving some properties (e.g., node degree distribution [START_REF] Maslov | Specificity and stability in topology of protein networks[END_REF]), or apriori reproducing, for example, the underlying modular structure of networks via stochastic block modeling [START_REF] Richard F Betzel | Diversity of meso-scale architecture in human and non-human connectomes[END_REF][START_REF] Faskowitz | Weighted stochastic block models of the human connectome across the life span[END_REF]. These approaches have effectively improved the characterization of brain networks with respect to random null-models and provided a statistical framework to estimate confidence intervals. Nevertheless, they remain quite limited in terms of the number of local network properties that can be simultaneously tested and modeled. This is a crucial aspect as in general we don't know how many and which are the local connection properties that have generated the global observed structure. Exponential random graph models (ERGMs) represent an intriguing solution to this limitation, as demonstrated by the recent development in the field. Broadly speaking, ERGMs belong to the family of maximum entropy (Max.Ent.) models [START_REF] Edwin | Information theory and statistical mechanics[END_REF][START_REF] Cimini | The statistical physics of real-world networks[END_REF]. Given an observed network G * , the ERGM probability distribution is the one with the largest information entropy that statistically reproduces some key structural properties, or features, of G * . The latter are associated with the model parameters, which can be inferred from the data.

A fortunate alignment among availability of experimental brain data, development of network theory and of powerful computational tools -mainly borrowed from social sciences [START_REF] David | Inference in curved exponential family models for networks[END_REF] -has allowed researchers to study the brain network structure and dynamics under the lens of the ERGMs. The goal of the present report is to track down these efforts and present them in a self-contained review with an unified accessible language. The ultimate ambition is to show their usefulness to a broad audience and encourage their exploitation to tackle current open questions in basic and clinical neuroscience.

The rest of the review is organized as follows. In (Sec.2) we introduce brain networks, briefly describe how to build them from experimental data and which are their main topological properties. In (Sec.3) we focus on a minimal version of the ERGM, we present both its theoretical and computational aspects, we linger on the interpretation of its parameters and introduce some widely used ERGM graph statistics. In (Sec.4) we review the existing literature of ERGM-based methods for neuroscience. We show how they have been used to characterize brain states, discriminate between different experimental conditions and model dynamic brain networks. In (Sec.5) the reader will find a discussion of the theoretical subtleties of the ERGMs with a focus on their applications to the brain. Finally, in (Sec.6) we conclude with a perspective on the future of the approach, and point out potentially useful alternatives.

Brain networks in a nutshell

The study of the brain from a network perspective is nowadays regarded as a vibrant interdisciplinary field of research, referred as network neuroscience [START_REF] Bassett | Network neuroscience[END_REF]. Such efforts have initially focused on quantifying the organizational properties of the nervous system by means of different network metrics [START_REF] Boccaletti | Complex networks: Structure and dynamics[END_REF][START_REF] Christopher | The physics of brain network structure, function and control[END_REF]. For example, it has been shown that brain networks tend to exhibit a pronounced modular structure i.e. they are partitioned in densely inter-connected communities linked by sparse inter-community connections [START_REF] Sporns | Modular brain networks[END_REF]. The high clustering together with the low average path length between pairs of nodes has suggested that neuronal interactions are organized in a small-world topology, optimizing the balance between segregation and integration of information [START_REF] Bassett | Small-world brain networks[END_REF][START_REF] Bassett | Small-world brain networks revisited[END_REF]. Moreover, there is mounting evidence on the presence of hubs [START_REF] Gong | Mapping anatomical connectivity patterns of human cerebral cortex using in vivo diffusion tensor imaging tractography[END_REF][START_REF] Achard | A resilient, low-frequency, small-world human brain functional network with highly connected association cortical hubs[END_REF] and metabolic constraints of the brain's wiring, with spatially closer regions supporting stronger patterns of connections [START_REF] Richard F Betzel | Generative models of the human connectome[END_REF][START_REF] Rosenbaum | The spatial structure of correlated neuronal variability[END_REF]. These brain network properties tend to span different spatial and temporal scales, are in common to several species, correlate with individual behavior and can be significantly affected by neurological disorders [START_REF] Richard | Multi-scale brain networks[END_REF][START_REF] Fornito | The connectomics of brain disorders[END_REF]. In the following, we will focus on brain networks derived from neuroimaging data resulting from the noninvasive recording of large-scale neural information in humans.

Building brain networks from experimental data

Neuroimaging broadly refers to a large set of techniques that allow to measure noninvasively the brain structure and function. The translation of these measurements into a network hinges on the definition of what nodes and edges are. In general, these definitions depend on experimental constraints as well as on the specific scientific question. Here, we will provide a gentle overview of the main approaches adopted so far in an effort to offer a minimal understanding of the nature of the networks that we will later explore. Should the readers be interested in more detailed technical descriptions, we refer them to some recent books and reviews [START_REF] De | Graph analysis of functional brain networks: practical issues in translational neuroscience[END_REF][START_REF] Jones | White matter integrity, fiber count, and other fallacies: the do's and don'ts of diffusion mri[END_REF][START_REF] Onerva Korhonen | Principles and open questions in functional brain network reconstruction[END_REF][START_REF] Fornito | Fundamentals of brain network analysis[END_REF].

A first type of brain networks aims at representing the structural wiring of the brain. Nodes are defined as different brain sites and edges represent their anatomical interconnections. To this purpose, diffusion tensor imaging (DTI) is a noninvasive technology that allows to infer the diffusion of water molecules through white matter tracts of the brain in the three-dimensional space [START_REF] Peter | In vivo fiber tractography using dt-mri data[END_REF]. The collected data are discretized into non-overlapping gray matter volumes (nodes) and the connection weights (edges) between two nodes are proportional to the anatomical properties of the fiber tracks (e.g. number, length).

The second type of brain networks aims instead at capturing the functional interactions between brain regions (nodes). These regions of interest (ROIs) typically correspond to cytoarchitectural landmarks known a-priori or more directly to sensors used for the measurements [START_REF] Gonzalez-Astudillo | Network-based brain-computer interfaces: principles and applications[END_REF]. For functional networks, interactions (edges) are defined as statistical similarities between the dynamics taking place in different nodes. Many technologies exist to record nodal dynamics as signals with different temporal and spatial resolutions. Among others, electroencephalography (EEG) and functional magnetic resonance imaging (fMRI) allow for noninvasive recordings of electrical activity through sensors placed over the scalp [START_REF] Craig | Electroencephalography: basic principles, clinical applications, and related fields[END_REF] and blood oxygen levels from 3D images of the brain, respectively [START_REF] Marcus E Raichle | Behind the scenes of functional brain imaging: a historical and physiological perspective[END_REF][START_REF] Salvador | Neurophysiological architecture of functional magnetic resonance images of human brain[END_REF]. Functional interactions are then inferred from the recorded signals using related measures such as Pearson or Spearman correlations, mutual information or Granger causality [START_REF] De | Graph analysis of functional brain networks: practical issues in translational neuroscience[END_REF][START_REF] Gonzalez-Astudillo | Network-based brain-computer interfaces: principles and applications[END_REF]. Thus, the study of functional networks is primarily intended to reveal the brain organizational properties from an information processing point of view.

At this stage, both structural and functional brain networks are characterized by an adjacency matrix A whose generic entry a ij is a scalar number, typically normalized between 0 and 1, representing the magnitude of the connection between nodes i and j. Thresholding procedures can be eventually adopted to filter out irrelevant links and obtain sparse unweighted networks. More details on why and how to filter information in complex brain networks can be found in [START_REF] Van Wijk | Comparing brain networks of different size and connectivity density using graph theory[END_REF][START_REF] De | A topological criterion for filtering information in complex brain networks[END_REF][START_REF] Garrison | The (in) stability of functional brain network measures across thresholds[END_REF].

The rationale of a statistical approach

The rationale for statistical modeling is based on the assumption that the observed system is a realization of a stochastic process, and that by characterizing the properties of this process we can gain insight into the underlying mechanisms governing the system's behavior. More specifically, the probability of observing each possible realization depends on endogenous system constraints, which can be either postulated or inferred from the data. This type of statistical uncertainty is different from that associated with exogenous experimental errors -here, the network construction process. Statistical uncertainty is due to the inherent variability of a system and is quantified using statistical measures. Uncertainty due to experimental errors arises from factors such as measurement noise and reflects the limitations of our measurement apparatus. In this review we deal with the former, discussion on the latter can be found elsewhere [START_REF] De | Graph analysis of functional brain networks: practical issues in translational neuroscience[END_REF].

In network neuroscience, a common first approach to assess the statistical significance of a network's properties has been to compare it a-posteriori to a null model. The latter is a simplified version of the network that preserves some basic structural features, but removes other non-random ones [START_REF] Váša | Null models in network neuroscience[END_REF]. By comparing a network to its null model, we can determine whether its observed properties are due to chance, or whether they reflect underlying mechanisms or constraints on the system. This comparison can be done using a variety of statistical tests, such as comparing the network's properties to those of a large ensemble of randomized networks, or comparing it to a null model that is specifically designed to test a particular hypothesis.

Such an approach however do not allow to test more specific hypotheses about the mechanisms that govern the network organization, e.g. the likelihood that a-priori defined structures are at the origin of the observed network. This kind of questions can be targeted by using instead an inferential approach. This involves (i) constructing a statistical model that describes the relationships between different variables in the system, and (ii) using this model to estimate the parameters that govern the behavior of the system. By comparing the predictions of the model to the observed data, we can determine whether the model accurately captures the underlying mechanisms that drive the system.

Furthermore complex systems typically arise from multiple organizing principles across different scales. Statistical (inferential) frameworks that can simultaneously account for these effects are crucial to understanding the global properties of the system. In network neuroscience, where the focus is on understanding the brain's structure and function, this approach is particularly useful [START_REF] Richard | Multi-scale brain networks[END_REF][START_REF] De | Multiscale topological properties of functional brain networks during motor imagery after stroke[END_REF][START_REF] Serrano | Extracting the multiscale backbone of complex weighted networks[END_REF]. By using inferential methods, we can identify the combinations of relevant network properties that give rise to the observed large-scale structures and predict how these properties will change under different conditions or perturbations.

Recent endeavors in exponential random graph models (ERGMs) have increasingly attracted the interest of the network science and neuroscience community due to their ability to match the aforementioned desiderata within a coherent and unifying framework. Hence, it is timely to discuss these emerging developments, and to seek to put them together into a common theoretical ground that can be used to tackle current open questions in modern neuroscience.

Exponential random graph model (ERGM)

The interest in the exponential class of probability distributions dates back to the dawn of modern statistics [START_REF] Bernard O Koopman | On distributions admitting a sufficient statistic[END_REF][START_REF] James | Sufficient statistics and intrinsic accuracy[END_REF] and was originally motivated by its properties with respect to sufficient statistics, as introduced by R.A. Fisher in the early '20s [START_REF] Fisher | On the mathematical foundations of theoretical statistics[END_REF][START_REF] Fisher | Theory of statistical estimation[END_REF]. A number of theoretical and computational properties makes them very appealing for the problem of statistical inference [START_REF] Barndorff-Nielsen | Information and exponential families: in statistical theory[END_REF][START_REF] Lawrence D Brown | Fundamentals of statistical exponential families: with applications in statistical decision theory[END_REF] and justifies their ubiquity in all branches of scientific research. In graph theory, the exponential family, here referred to as exponential random graph models (ERGMs), came on stage in the early '80s [START_REF] Paul | An exponential family of probability distributions for directed graphs[END_REF][START_REF] Frank | Markov Graphs[END_REF][START_REF] Strauss | On a general class of models for interaction[END_REF] building on the seminal work of J. Besag [START_REF] Besag | Spatial Interaction and the Statistical Analysis of Lattice Systems[END_REF] and further developing since then [START_REF] Wasserman | Logit models and logistic regressions for social networks: I. an introduction to markov graphs andp[END_REF][START_REF] Anderson | A p* primer: Logit models for social networks[END_REF]. Analogies with well-known methods of statistical mechanics have been recently explored in simple cases in conjunction with the outbreak of network science [START_REF] Bianconi | Entropy of network ensembles[END_REF][START_REF] Radicchi | Classical information theory of networks[END_REF][START_REF] Cimini | The statistical physics of real-world networks[END_REF]. ERGMs have recently become popular in several scientific fields e.g. epidemiology [START_REF] Matthew | The application of statistical network models in disease research[END_REF], sociology [START_REF] Mark | Statistical models for social networks: Inference and degeneracy[END_REF][START_REF] Robins | An introduction to exponential random graph (p*) models for social networks[END_REF][START_REF] Lusher | Exponential random graph models for social networks: Theory, methods, and applications[END_REF] as well as neuroscience. As a consequence an increasing number of publicly available related softwares has appeared, such as the R language packages collected in the statnet suite [START_REF] Mark S Handcock | statnet: Software tools for the representation, visualization, analysis and simulation of network data[END_REF][START_REF] Steven M Goodreau | A statnet tutorial[END_REF].

Our reference implementation throughout this section is the ergm package [START_REF] David R Hunter | ergm: A package to fit, simulate and diagnose exponential-family models for networks[END_REF].

Here, we give the theoretical minimum to familiarize the reader with ERGMs and their estimation methods. Notably, we will consider unweighted and undirected graphs as they represent the most common scenario in the current literature. Later in Sec.4 we will extend the discussion to more sophisticated cases, when appropriate.

Let G be a set of all finite graphs of N nodes with no self-loops, containing at most one single edge between two nodes. Each graph G ∈ G can be equivalently represented by a N × N adjacency matrix A containing Boolean values, i.e. a ij = {1, 0}.

At the hearth of the ERGM approach there is the definition of an appropriate probability mass function P (G) over the ensemble G . The crucial idea of ERGMs is to encode the relevant information about the graph G in a vector x(G) ∈ R r of r statistics or metrics, where each element x α (G) measures a different network property of interest.

Specifically, the probability of observing a generic graph G reads as

P (G|θ) = e θ•x(G) G∈G e θ•x( G) , (1)
where θ ∈ R r are the model parameters weighting the graph statistics and • is the dot product. A formal analogy with the canonical Boltzmann distribution is readily recognized by defining the Hamiltonian

H(G) = -θ • x(G) = - α θ α x α (G) (2)
Hence, the denominator of Eq.1 is a normalizing constant, which turns out to play the role of partition function Z and we can rewrite P (G|θ) = 1 Z e -H(G) . Given a specific observed network G * , one can infer the model parameters θ * so that the expected value of each graph statistics over the ensemble

G ⟨x⟩ = G∈G x(G)P (G|θ * ) (3)
statistically matches the observed value ≃ x(G * ). Note that (Eq.3) sets the number of model parameters from 2 ( N 2 ) -the number of possible edges, in the worst case -to a significantly smaller amount determined by the chosen graph statistics. By consequence, in general a tie-level matching between the networks generated with P (G|θ * ) and the observed graph G * is neither expected nor desired.

Model parameter estimation

Apart from very simple cases, ERGM parameters are hard to obtain in a closed form. In general, the nature and the number of graph statistics make the denominator of Eq.1 impossible to compute analytically. The evaluation of the partition function Z is a very well-known problem in many situations, like for example when trying to infer parameters from Boltzmann-like distributions [START_REF] Chau Nguyen | Inverse statistical problems: from the inverse ising problem to data science[END_REF][START_REF] Cocco | Inverse statistical physics of protein sequences: a key issues review[END_REF].

In the last 30 years there has been a tremendous theoretical and computational effort in developing methods to address this issue and achieve efficient estimation algorithms. In the following, we present a brief introduction to the key-idea behind many of these methods [START_REF] Charles | Markov chain monte carlo maximum likelihood[END_REF][START_REF] Charles | Constrained monte carlo maximum likelihood for dependent data[END_REF].

Given a graph G * and a set of statistics x(G), we can take a Bayesian perspective and argue that our knowledge on the parameters θ is better described by the posterior distribution

P (θ|G * ) = P (G * , θ) P (G * ) = P (G * |θ)P (θ) P (G * ) . (4)
In the case where no prior information is available for θ, we can take P (θ) to be uniformly distributed. From Eq.4 we see that the posterior distribution P (θ|G * ) is directly proportional to the distribution of the data G * given the parameters θ. Accordingly, our best guess on the parameters value is given by the Maximum Likelihood Estimator (MLE)

θ * = arg max θ L (G * |θ) , (5)
where

L (G * |θ) = log P (G * |θ) is the log-likelihood function.
The evaluation of (Eq.5) is hampered by the computation of the partition function Z inside P (G * |θ). As discussed above, this entails the need for non-exact numerical methods.

Rather than maximizing L (G * |θ) directly, the idea is to maximize a shifted loglikelihood L (

G * |θ) = L (G * |θ) -L (G * |θ 0 )
, where θ 0 is an arbitrarily parameter vector. ‡ A bit of algebra reveals that the new log-likelihood function can be written as

L (G * |θ) = (θ -θ 0 ) • x(G * ) -log ⟨e (θ-θ0)•x(G) ⟩ θ0 (6)
where the subscript ⟨•⟩ θ0 indicates the expectation value over the distribution P (G|θ 0 ). Note that the new log-likelihood L (G * |θ) is maximized by the same θ * that maximizes L (G * |θ).

In practice, the evaluation of the partition functions Z is bypassed and we just need to calculate the expectation ⟨•⟩ θ0 , which can be efficiently performed with standard Markov-chain Monte-Carlo (MCMC) approximations. Eventually, the whole procedure is repeated until the convergence to some stable solution θ * is reached (Tab. 1). In the following, we will refer to this method as MCMC-MLE.

ERGMs in practice

The increasing diffusion of ERGMs has been driven by a concurrent technical effort to develop efficient algorithms for carrying out the inference task. By far the most popular implementation available to date is the statnet suite of R-based software packages [START_REF] Mark S Handcock | statnet: Software tools for the representation, visualization, analysis and simulation of network data[END_REF][START_REF] Mark | Statnet: Tools for the statistical modeling of network data[END_REF]. To our knowledge, no other implementation has reached the same level of maturity and extending these techniques outside the R language is an open challenge for the next future. The reference implementation for the vast majority of existing packages is the ergm package [START_REF] David R Hunter | ergm: A package to fit, simulate and diagnose exponential-family models for networks[END_REF]. A number of other packages have been built on it, both for static [START_REF] Caimo | Bayesian inference for exponential random graph models[END_REF][START_REF] James D Wilson | Stochastic weighted graphs: Flexible model specification and simulation[END_REF][START_REF] Schweinberger | Hergm: Hierarchical exponential-family random graph models[END_REF][START_REF] Pavel | Inference for social network models from egocentrically sampled data, with application to understanding persistent racial disparities in hiv prevalence in the us[END_REF] and temporal networks [START_REF] Pavel | A separable model for dynamic networks[END_REF][START_REF] Leifeld | Temporal exponential random graph models with btergm: Estimation and bootstrap confidence intervals[END_REF].

An illustrative workflow based on the ergm package is implemented in [START_REF] Dichio | ergm minimal[END_REF], consisting of six fundamental steps. They include: i) creating a network object using the network package from the adjacency matrix and nodal/edge covariates (if any), ii) computing ERGM statistics for the experimental network (optional), iii) specifying ‡ In practice θ 0 should be close enough to θ to ensure convergence in a reasonable computational time. One popular choice is to take θ 0 as the maximum pseudo-likelihood estimation (MPLE) of θ under the additional hypothesis that the edges are mutually independent [START_REF] David | Inference in curved exponential family models for networks[END_REF][START_REF] Bruce | Statistical inference for valued-edge networks: The generalized exponential random graph model[END_REF].

θ (0) = θ P L τ = 0 while (conv = FALSE) do τ += 1 generate n graphs G k ∼ P (G|θ (τ -1) ) by MCMC θ (τ ) = arg max θ (θ -θ (τ -1) ) • x(G * ) -log( 1 n n k=1 e (θ-θ (τ -1) )•x(G k ) ) if (convergence criterion) then conv = TRUE end if end while
Table 1: MCMC-MLE pseudo-code. The initial choice of parameters θ (0) is arbitrary, here the pseudo-likelihood value θ P L . Several convergence criteria are possible, from the simplest ∥θ (τ ) -θ (τ -1) ∥ < ϵ for some ϵ > 0, up to more refined choices [START_REF] David | Inference in curved exponential family models for networks[END_REF][START_REF] Pavel N Krivitsky | Ergm 4: Computational improvements[END_REF]. and estimating the model, iv) assessing MCMC convergence (optional), v) evaluating the model fit using goodness of fit methods (GoF), and vi) simulating the model using estimated parameters (optional).

The computational time required for the model estimation might vary from a few seconds to several hours, depending on a number of factors including the type and number of graph statistics, the basic properties of the networks (e.g. size, sparsity) and several estimation settings that can be possibly adjusted. Because of the intrinsic characteristics of the current neuroimaging technology (Sec.2), brain networks of interest in this review have a relatively small number of nodes (at most hundreds). For such sizes, the ERGM estimation is in general doable. For further information regarding the computational aspects and a comprehensive overview of the latest technical enhancements of the ergm 4.0 package, we refer to [START_REF] Pavel N Krivitsky | ergm 4: New features for analyzing exponential-family random graph models[END_REF][START_REF] Pavel N Krivitsky | Ergm 4: Computational improvements[END_REF] and references therein.

Interpretation of ERGMs

A crucial, often overlooked, aspect in ERGM is the interpretation of the parameters θ in Eq.1. Let us address the question starting with a simple example [START_REF] Frank | Markov Graphs[END_REF][START_REF] Park | Statistical mechanics of networks[END_REF][START_REF] Albert | Statistical mechanics of complex networks[END_REF]. Consider an ERGM based on a single graph metric x(G) corresponding to the number of edges in the network, i.e. x(G) = i<j a ij . We first consider the forward problem in which there is no observed graph but instead the parameter θ is given and the goal is to evaluate the expected value for x(G).

In such a simple case, the mass probability function can be computed exactly:

P (G|θ) = e θx(G) i<j 1 aij =0 e θaij = e θx(G) (1 + e θ ) ( N 2 ) . (7)
Notably, by defining

p(θ) = e θ 1 + e θ (8)
we can rewrite

P (G|θ) = p x(G) (1 -p) ( N 2 )-x(G)
, which is the probability of a Bernoulli graph with N nodes and x(G) links [START_REF] Erdos | On the evolution of random graphs[END_REF]. The expected value for the number of edges in the graph is then simply ⟨x⟩ = N 2 p.

In practice, the most common scenario is quite the opposite. One does not know the value of θ but rather has an observed graph G * with, say, M edges. In these situations, the goal is to solve the inverse problem by inferring θ * from x(G * ) = M . Taking p * = M/ N 2 and inverting Eq.8:

θ * (p * ) = log p * 1 -p * (9)
so that, consistently, ⟨x⟩ = N 2 p = M . Notably, Eq.9 tells us that if we start from a Bernoulli graph with p * = 0.5, then θ = 0. This is the maximally random case since each dyad corresponds to the toss of a fair coin. But if we start from a denser (p * > 0.5) or sparser (p * < 0.5) graph, then θ * > 0 or θ * < 0, respectively.

In such a simple case we have therefore a complete understanding of both the forward and inverse problems. In real situations, however, the interpretation of the parameters becomes trickier since ERGMs may include several statistics, often exhibiting some degree of dependence. As a consequence, there might be interactions among the model parameters that are in general difficult to compute analytically. Yet, inference is still possible via approximate numerical methods such as the MCMC-MLE. The interpretation of the inferred parameters can be given by quantifying their effect on the likelihood that an edge between two nodes exists or not in the network. For example, consider a single edge toggle between nodes i and j. Let us call P (G aij =1,a \ij |θ * ) the probability of having an edge between the two nodes given the rest of the graph a \ij , and P (G aij =0,a \ij |θ * ) its complementary.

From Eq.1 one easily finds that

P (G aij =1,a \ij |θ * ) P (G aij =0,a \ij |θ * ) = exp θ * • ∆ x ij (G) , (10)
where the

∆ x ij (G) = x(G aij =1,a \ij ) -x(G aij =0,a \ij
) are the so-called change statistics. From Eq. 10 is clear that the probability of a link to exist depends on the magnitude and signs of both θ * and ∆ x ij (G). Let us consider a very simple case where the presence of a new link between i and j increases the value of a given graph statistic i.e. ∆ x ij (G) > 0. If θ * > 0, then this change statistic favors the probability of this edge to exist in the graph (a ij = 1). Instead, if θ * < 0 the same change statistic penalizes the existence of the connection, favouring instead (a ij = 0). Overall, the magnitude of this effect is given by the

term exp [θ * ∆ x ij (G)].
In other words, one can disentangle and quantify the relative influence of each graph statistic on the probability to observe, or not, an edge between two nodes i and j.

Note that the above micro-level interpretation of the parameters is precisely at the core of the MCMC routines used to sample networks from ERGM probability distributions. Markov chains of graphs whose stationary distribution is P (G|θ * ) are obtained by sequentially proposing changes on a graph, evaluating the probability of the new configuration by Eq.10 and accepting it according to the Metropolis-Hastings recipe. We refer to recent reviews in the field for a more detailed description of the MCMC methods used in ERGMs [START_REF] David | Inference in curved exponential family models for networks[END_REF].

Graph statistics for ERGMs

A broad spectrum of graph statistics, or metrics, x(G) have been formulated and included in Eq.1. Here, we present those that have been more frequently adopted in neuroscience.

A first group of graph statistics includes the so-called dyadic independence terms, i.e. combinations of only single-dyad terms. The most intuitive metric is the edges term which measures the number of edges in the graph (Fig. 2a):

x e (G) = i<j a ij . (11)
A straightforward extension is the so-called edge covariate

x dc (G) = i<j a ij ζ ij . ( 12
)
where ζ ij is an attribute associated defined for each dyad, such as the Euclidean distance in a spatial network. ERGM terms based on nodal attributes also belong to this group. If η i is a categorical property of the nodes, then it is possible to define the so-called nodematch term as

x nm (G) = i<j δ ηiηj a ij , ( 13
)
where δ is the Kronecker delta and x nm counts the number of edges whose nodes are labeled by the same categorical attribute. For instance, η i = R, L could indicate whether a node belongs to the right (R) or left (L) hemisphere. More in general, real networks are characterized by the presence of complex connectivity structures leading to dependencies between dyads. Therefore, we now turn into presenting graph statistics that involve products of two or more dyadic variables a ij . Following the seminal work of O. Frank and D. Strauss on Markov Graphs [START_REF] Frank | Markov Graphs[END_REF], we assume that only dyads that share a node can be dependent or, equivalently, that nonincident dyads are conditionally independent.

Under this general definition, several statistics may be introduced based on the count of local connection patterns. For instance, the presence of hubs, documented in many real-world networks [START_REF] Albert | Statistical mechanics of complex networks[END_REF][START_REF] Albert | Statistical mechanics of complex networks[END_REF], can be measured in an ERGM by the k-stars term:

x

(k) st (G) = 1 k! i0 • • • i k a i0i1 . . . a i0i k (14)
with k = 1, . . . , N -1.

The global tendency of networks to form clusters, can be instead measured by the presence of triplets, or triads, of connected nodes. A k-triangles term can be then defined as:

x (k) tr (G) = 1 ξ(k) i0 • • • i 1+k a i0i1 a i0i2 a i1i2 . . . a i0i 1+k a i1i 1+k (15)
where the symmetry factor is ξ(k) = 3! for k = 1 and ξ(k) = 2k! for k = 2, . . . , N -2.

The simplest case x

(1)

tr (G) = 1 6
i0i1i2 a i0i1 a i0i2 a i1i2 counts the number of triangles in the graph.

𝑥 ! = 𝑥 " ($) = 𝑥 &'( ($) = 𝑥 !'( ($) = i 3 2 1 i 3 1 2 i 3 1 2 (a) (b) (d) (c)

3.3.1.

Degeneracy of simple ERGMs. Starting from the earliest numerical investigations, it has been observed that for some basic ERGM specifications, MCMC simulations hardly "mix", meaning they do not reach a stationary distribution within a reasonable amount of time [START_REF] Strauss | On a general class of models for interaction[END_REF][START_REF] Tom | Markov chain monte carlo estimation of exponential random graph models[END_REF]. The lack of mixing hampers the estimation task. This behavior is due to the geometry of ERGM stationary distributions. As it turns out, there exist choices of statistics that result in "degenerate" models, meaning that assign a significant probability only to a few graph configurations, often with a very different structure, such as empty and fully connected graphs. [START_REF] Mark S Handcock | Assessing degeneracy in statistical models of social networks[END_REF][START_REF] Rinaldo | On the geometry of discrete exponential families with application to exponential random graph models[END_REF]. A general, thorough discussion of the causes and consequences of degeneracy in the context of ERGMs can be found in [START_REF] Schweinberger | Instability, sensitivity, and degeneracy of discrete exponential families[END_REF].

For a few simple ERGMs, such as the 2-star model

(H = x e + x (2)
st ) or the Strauss clustering model (H = x e + x (1) tr ), it has been demonstrated through analytic solutions that this phenomenon corresponds to a standard phase separation observed in statistical mechanics between a low-and high-density phase, representing almost empty and full networks, respectively [START_REF] Park | Solution of the two-star model of a network[END_REF][START_REF] Park | Solution for the properties of a clustered network[END_REF]. For such model specifications, there is no parameters combination that is able to reproduce intermediate densities of edges, two-stars and triangles as observed in real-world networks.

These examples demonstrate that degeneracy issues arise when the model includes unstable sufficient statistics §. Such models show excessive sensitivity in the sense that § According to the definition in [START_REF] Schweinberger | Instability, sensitivity, and degeneracy of discrete exponential families[END_REF], x(G) is stable if there exist C > 0 and

Mc > 0 such that Ux ≡ max G∈G x(G) ≤ CM ∀M > Mc, where M = N
2 is the number of degrees of freedom in the network. Note that for instance Ux e = M implies that the edges term is stable while U

x (2) st ∼ N M
small state changes can result in extremely large odds ratios eq.( 10) [START_REF] Schweinberger | Instability, sensitivity, and degeneracy of discrete exponential families[END_REF].

The main conclusion of these theoretical investigations is that inference for simple ERGMs is simply not possible, ill-posed. To solve degeneracy issues, a possible solution is to add structure to the model specification by incorporating additional information such as i) edge-specific or node-specific covariates, ii) block structure, iii) temporal structure, iv) multilevel structure, or v) spatial structure. A recent comprehensive discussion of these cases can be found in [START_REF] Schweinberger | Exponential-family models of random graphs: Inference in finite, super and infinite population scenarios[END_REF].

However, the most common approach to ERGM degeneracy issues in neuroscience applications has been so far to include the so-called "curved statistics" in the model specification. In the next section, we will discuss this approach in detail.

Curved Statistics.

The inclusion of curved statistics into ERGM model specifications has been found to be effective in alleviating degenerate behaviors [START_REF] Schweinberger | Instability, sensitivity, and degeneracy of discrete exponential families[END_REF][START_REF] Schweinberger | Exponential-family models of random graphs: Inference in finite, super and infinite population scenarios[END_REF]. These statistics are linear combinations of the complete distribution of degree or shared partner statistics, where the coefficients are geometrically weighted. Therefore, these statistics not only enumerate patterns but also include additional non-linear constraints on graph structures through the use of geometric coefficients. [START_REF] David | Inference in curved exponential family models for networks[END_REF][START_REF] Tom Ab Snijders | New specifications for exponential random graph models[END_REF].

For instance, the geometrically weighted degree (gwd) term reads:

x gwd (G|τ ) = e τ N -1 k=1 1 -1 -e -τ k x (k) d (G) , (16)
where τ > 0 is a parameter and x

(k) d (G)
is the number of nodes in the graph G whose degree is exactly equal to k (Fig. 2b). Similarly, the geometrically weighted non-edgewise shared partner (gwnsp) statistic:

x gwnsp (G|τ ) = e τ N -2 k=1 1 -1 -e -τ k x (k) nsp (G) , (17)
where τ > 0 and x

(k)
nsp (G) is the number of non connected dyads having k neighbors in common, (Fig. 2c). Finally, the geometrically weighted edgewise shared partner (gwesp) statistics:

x gwesp (G|τ ) = e τ N -2 k=1 1 -1 -e -τ k x (k) esp (G) , (18)
where τ > 0 and x

(k)
esp (G) is now the number of connected dyads that share exactly k neighbors (Fig. 2d).

To provide an interpretation of curved statistics and better understand the underlying modeling assumptions, let's examine the case of gwesp and proceed as in (Sec.3.2). Namely, let G be a graph and consider two connected nodes i, j with k shared partners. Let G ′ be the same graph as G but with an additional common neighbor between the nodes i, j, this implies

x (k) esp (G ′ ) = x (k) esp (G) -1 x (k+1) esp (G ′ ) = x (k+1) esp (G) + 1 , (19)
implies that the 2-star term is unstable.

where x

(k)
esp (G) is the number of connected dyads that share exactly k neighbors. For an ERGM including a gwesp term, we find

P (G ′ |θ) P (G|θ) ∝ exp θ gwesp x gwesp (G ′ |τ ) -x gwesp (G|τ ) ∝ exp θ gwesp (1 -e -τ ) k ≡ ρ . (20)
We start by noting that since τ > 0, we have 1 -e -τ ∈ (0, 1). Moreover: (i) the scaling factor θ gwesp > 0 implies ρ > 1, hence preference for adding shared partners; θ gwesp < 0 on the contrary implies a preference for deleting shared partners. This interpretation is akin to that of all standard ERGM parameters as discussed in (Sec.3.2). (ii) lim k→∞ ρ = 1 i.e. adding/removing shared partners to pairs of connected nodes with already many of them (k ≫ 1) has little effect on the ratio P (G ′ )/P (G). This implies an effective cutoff for the order k of x

(k)
esp statistic that is relevant for the system. (iii) The speed of the geometric decay is controlled by an external parameter τ , small values imply a rapid decay while high values imply a slow decay. Hence, the gwesp statistic is a refined version of the simple sum of triangles. When the parameter θ gwesp is positive, the general tendency for clustering will not push the system into a fully connected state. In fact, its implementation through the gwesp statistic implies a strong advantage for the addition of shared partners for poorly connected nodes, while almost no effect for nodes that already have many connections. This is indeed a more realistic behavior of a triadic closure, which in the case of brain networks it is expected to be constrained by the inner physical and functional resources [START_REF] Bullmore | The economy of brain network organization[END_REF].

The same reasoning applies to other curved statistics, such as gwd and gwnsp. These two statistics are alternatively used to test the hypothesis that shortest paths results from the presence of hubs or two-paths. The former would affect the node degree distribution, the latter the distribution of non-edgewise shared partners. Thus, via the gwd and gwnsp statistics, the presence of short paths in the network will result in positive θ gwd and θ gwnsp values, without the system ending up in a fully connected/empty state.

ERGMs in neuroscience

ERGMs have started to be exploited as statistical framework to study brain networks since 2010 [START_REF] Van Wijk | Comparing brain networks of different size and connectivity density using graph theory[END_REF]. The typical workflow of an ERGM analysis in neuroscience is summarized in Fig. 3. In this section, we present an overview of the most recent results obtained in the last decade. A schematic summary of the existing literature can be found in Tab.2.

While ERGMs are mathematical abstractions and can be applied to any network, most research has so far focused on functional brain connectivity derived from fMRI and EEG data. In the following, we structure the presentation of the current literature around three main scientific questions: i) what are the basic network mechanisms of healthy brain functioning (Sec.4.1), ii) how to statistically compare brain networks between different states (Sec.4.2), and iii) how to explicitly characterize the temporal network evolution following a pathological event (Sec. [START_REF] Bassett | Network neuroscience[END_REF] Table 2: Recent ERGM studies of human brain networks. We summarize each paper by its reference, the kind of experimental method for the data and the average number of nodes of the resulting networks, ERGM-method employed and its R implementation.

Minimal model of brain networks

At the core of many ERGM implementations for neuroscience there is the model first proposed by Simpson et al. in two consecutive works [START_REF] Simpson | Exponential random graph modeling for complex brain networks[END_REF][START_REF] Simpson | An exponential random graph modeling approach to creating group-based representative whole-brain connectivity networks[END_REF], where the H-function reads as

-H(G|θ) = θ 1 x e (G) + θ 2 x gwesp (G|τ ) + θ 3 x gwnsp (G|ν) (21)
The edges term x e (G) characterizes the density of the network. This statistic is typically used to ensure that ERGMs also reproduce the actual number of connections in the observed network. The gwesp term x gwesp (G|τ ) is meant to capture the overall network clustering, a property that is crucial for the segregation of information in the brain [START_REF] Bassett | Small-world brain networks[END_REF][START_REF] Bassett | Small-world brain networks revisited[END_REF]. The gwnsp term x gwnsp (G|τ ) is related to the presence of two paths in the network and can be regarded as a measure of integration of information in the brain [START_REF] Power | Evidence for hubs in human functional brain networks[END_REF]. Hence, Eq.21 defines a parsimonious model of the fundamental properties of brain networks, i.e. connectedness, clustering, and global-efficiency.

In [START_REF] Simpson | Exponential random graph modeling for complex brain networks[END_REF][START_REF] Simpson | An exponential random graph modeling approach to creating group-based representative whole-brain connectivity networks[END_REF] the authors applied this model to fMRI brain networks from 10 healthy subjects, fitted separately per each individual. The same model was also used in [START_REF] Bcl Lehmann | Characterising group-level brain connectivity: a framework using bayesian exponential random graph models[END_REF] to characterize fMRI brain networks from 200 healthy individuals. All these works agreed on assigning negative values to the parameters associated with edges and gwnsp statistics and positive values to gwesp statistics, while the decay parameter for the last two cases was fixed to τ = ν = 0.75 [START_REF] Simpson | Exponential random graph modeling for complex brain networks[END_REF]. The edges parameter θ 1 < 0 confirms the tendency of the network to be sparse. The gwnsp parameter θ 3 < 0 indicated that any tendency for global integration of information in the brain is unlikely to be supported by shortest paths with length equal to 2. The gwesp parameter θ 2 > 0 showed that brain networks are statistically inclined to form clusters, which is a basic hallmark of brain segregation of information.

In subsequent studies, [START_REF] Paul E Stillman | Statistical modeling of the default mode brain network reveals a segregated highway structure[END_REF][START_REF] Paul E Stillman | A consistent organizational structure across multiple functional subnetworks of the human brain[END_REF] investigated whether and how the aforementioned x(G i ) ∼ x(G * ). Goodness of fit (GoF) methods can be used to assess the performance of the estimated model in reproducing other features of the input network [START_REF] David R Hunter | Goodness of fit of social network models[END_REF].

whole-brain statistical network properties were also present in specific subsystems associated with basic behavioral functions. Resting-state fMRI networks were then constructed and fitted separately in 21 healthy subjects and for 8 different subsystems, namely auditory, subcortical, dorsal attention, ventral attention, salience, cinguloopercular task control, fronto-parietal task control, default mode [START_REF] Power | Functional network organization of the human brain[END_REF] (Fig. 4).

To carry out this analysis, the authors introduced the so-called correlation generalized ERGM (cGERGM). This formulation of the model is still based on the ERGM in Eq.21, with a slightly different geometrical weighting strategy of the graph statistics [START_REF] James D Wilson | Stochastic weighted graphs: Flexible model specification and simulation[END_REF][START_REF] Paul E Stillman | A consistent organizational structure across multiple functional subnetworks of the human brain[END_REF]. However, cGERGM differs from the standard approach in two substantial ways. First, cGERGM was designed to handle full weighted correlation matrices, instead of dealing with unweighted sparse networks. Secondly, additional geometric information on the nodes were modeled externally to the ERGM by using a beta regression of the group-averaged link weights on the graph statistics [START_REF] Paul E Stillman | Statistical modeling of the default mode brain network reveals a segregated highway structure[END_REF]. In Asterisks indicate that the estimation was significantly different from 0 for all of the 21 subjects involved in the study. In general, larger networks (cingulo-opercular, salience, default mode, fronto-parietal) yield more robust estimations i.e. significant for a larger number of subjects. The corresponding brain subsystems are shown along the x-axis for illustrative purposes, as defined in [START_REF] Gu | Controllability of structural brain networks[END_REF].

particular, two kinds of geometric information were included: a nodematch term (Eq.13) to measure the tendency of connected nodes to lie in the same hemisphere, and a edge covariate term (Eq.12) to measure the total Euclidean distance between all connected nodes. To emphasize the difference from the ERGM parameters, we shall call β 4 , β 5 the ones controlling for these effects. Across all sub-systems, and almost all subjects, authors found consistent results with parameter values statistically different from 0 especially for the larger subsystems, i.e. cingulo-opercular, salience, default mode, fronto-parietal. The sign of the θ 2 > 0 and θ 3 < 0 parameters associated respectively with gwesp and gwnsp statistics, were in line with those reported for the entire brain network. Instead, the edge statistics was associated with a positive parameter θ 1 > 0. This result is consistent with the fact that brain subsystems are more densely connected internally than between each other [START_REF] Nikola T Markov | Cortical high-density counterstream architectures[END_REF][START_REF] Sporns | Modular brain networks[END_REF].

The estimation of the regression parameters associated with the spatial geometry of the brain structural networks turn out to be much more inconsistent across subsystems and subjects. The significantly positive hemispheric parameter (β 4 > 0) indicates that the nodes in the fronto-parietal and default-mode subnetwork are more likely to be connected within hemispheres. Instead, there is a general weak tendency of all subnetworks, but the auditory one, to penalize long distant connections (β 5 < 0, albeit significant for only between 20% and 62% of the individuals).

The overall emerging picture is that of a modular brain network reflecting a significant clustering behavior together with a strong intrahemispheric connectivity. Note that the negative gwnsp parameters do not necessarily imply the absence of 2-paths in the networks. Instead, they indicate that short paths of length 2 are relatively underrepresented with respect to what could be expected when considering other statistics alone e.g., triangles.

A related question is to what extent integration of information is statistically supported by short paths in the brain network. From a biological perspective, it would be more plausible that information flows through hubs and not via short paths, although both mechanisms are not mutually exclusive [START_REF] Deco | Rethinking segregation and integration: contributions of whole-brain modelling[END_REF]. To address this question, [START_REF] Obando | A statistical model for brain networks inferred from large-scale electrophysiological signals[END_REF] considered EEG resting-state networks in a group of 108 healthy subjects. In particular, they compared two different ERGMs, one including gwesp and gwnsp statistics, and the other including gwesp and gwd statistics, that is

-H(G|θ) = θ 2 x gwesp (G|τ ) + θ 3 x gwd (G|ν); ( 22
)
while holding fixed the number of edges in the network x e (G) ≡ x e (G * ) and τ = ν = 0.75.

Results showed that the second ERGM was able to better reproduce brain network properties not included in the model such as modularity, global-and local-efficiency. Notably, both θ 2 and θ 3 values were in average positive and significant confirming their statistical relevance for brain networks. Taken together, these findings support the hypothesis that while information segregation emerges from triangles, information integration could be actually mediated by hubs rather than shortest paths.

4.1.1.

Group-representative networks. Thus far, ERGMs have been used to characterize brain networks obtained from different individuals, separately. The resulting parameters distribution from many subjects can be then used to assess statistical properties at the group-level via standard tools such as hypothesis testing and regression analysis. Alternatively, assuming that brain networks from different individuals are stochastic realizations of the same system, they can be simultaneously modeled so as to obtain a unified group description. This is the rationale behind the construction of a group representative network (GRN) i.e. a network that summarizes the statistical properties of a given ensemble (Fig. 5).

In network neuroscience, this question has been mainly addressed by aggregating the connectivity matrices of different individuals to obtain a group-average (mean-GRN ) or group-median (median-GRN ) representative network [START_REF] Achard | A resilient, low-frequency, small-world human brain functional network with highly connected association cortical hubs[END_REF][START_REF] Song | Default network and intelligence difference[END_REF]. Alternatively, methods to assess the statistical significance of the pooled connection values for each pair of nodes can be adopted [START_REF] Zalesky | Network-based statistic: identifying differences in brain networks[END_REF]. Despite being computationally straightforward, the major limitation of these methods is that they treat edges independently and cannot directly inform on complex connection properties involving more than 2 nodes. An interesting possibility is therefore to incorporate the intersubject variability directly into the ERGM formulation.

In [START_REF] Simpson | An exponential random graph modeling approach to creating group-based representative whole-brain connectivity networks[END_REF], the authors proposed to generate a GRN starting from the ERGMs in (Eq.21) fitted separately on different subjects. Then, a group representative network is simulated via a new ERGM with the H-function -H(G|θ) = µ 1 x e (G) + µ 2 x gwesp (G|τ ) + µ 3 x gwnsp (G|ν) [START_REF] Rosenbaum | The spatial structure of correlated neuronal variability[END_REF] where τ = ν = 0.75 are fixed and the µ are the group-averaged values of the individual model parameters. This procedure better reproduced several topological properties of the actual brain networks, as compared to standard ERGMs fitted on the groupaveraged network. Note that this is a frequentist approach, the estimation of the group level parameters are point-like averages. We refer to this GRN construction scheme as a mean-ERGM.

The idea of summarizing individual estimates in a vector of hyper-parameters was pushed further in [START_REF] Bcl Lehmann | Characterising group-level brain connectivity: a framework using bayesian exponential random graph models[END_REF]. Accordingly, the model parameter vector for each individual is a realization of a normal probability distribution θ ∼ N (µ, Σ θ ). The hyper-parameters (µ, Σ θ ) are estimated by extending to a multilevel hierarchical setting a Bayesian formulation of the ERGM (BERGM), developed in [START_REF] Caimo | Bayesian inference for exponential random graph models[END_REF]. A GRN can then by obtained by sampling a network from the posterior distribution of parameters. This method is named multi-BERGM.

Identification of discriminant brain network features

ERGMs can also be used to assess statistical differences between brain networks obtained under different experimental conditions or belonging to different groups of subjects.

In [START_REF] Obando | A statistical model for brain networks inferred from large-scale electrophysiological signals[END_REF] the authors analyzed EEG networks with 56 nodes in 108 healthy subjects under two different conditions: eyes closed (EO) and eyes closed (EC) resting-states. To this end, they used the same ERGM defined in Eq.22. In particular, separate ERGMs were estimated from brain networks extracted from each individual, condition and frequency band theta (4 -7 Hz), alpha (8 -13 Hz), beta (14 -29 Hz), gamma (30 -40 Hz). The estimated parameter values were then compared between the two conditions. Results showed that θ 2 values, associated with the gwesp statistics, were significantly higher in EO than EC, for both alpha and beta bands. This result is consistent with the local increase of EEG activity widely reported in literature for those bands [START_REF] Robert | Eeg differences between eyes-closed and eyes-open resting conditions[END_REF][START_REF] Yu | Eeg differences between resting states with eyes open and closed in darkness[END_REF]. Notably, the difference in the beta band could not be detected by a standard network analysis.

ERGM-based analyses have been also adopted to characterize aging. To this end, in [START_REF] Bcl Lehmann | Characterising group-level brain connectivity: a framework using bayesian exponential random graph models[END_REF] the authors compared fMRI brain networks in a group of young (Y: 18 -33 ) and elderly (O: 74 -89) human subjects [START_REF] Bcl Lehmann | Characterising group-level brain connectivity: a framework using bayesian exponential random graph models[END_REF]. Starting from the hierarchical Bayesian ERGM presented in (Sec.4.1.1), a further hierarchical layer is added to account for the population partitioned in two subgroups, yielding two different hyper-parameters ϕ (Y ) = (µ (Y ) , Σ θ ) and ϕ (O) = (µ (O) , Σ θ ) and same covariance structure is assumed for the two groups.

The results showed that the estimates for the parameter associated with the gwnsp statistics were significantly higher in the young group compared to the old one (Fig. 6). These findings indicated that functional brain networks are prone to significant loss of integration of information across the lifespan, consistently with previous findings reporting a decreasing global-efficiency with age [START_REF] Achard | A resilient, low-frequency, small-world human brain functional network with highly connected association cortical hubs[END_REF]. One might then wonder whether similar changes occur in structural brain networks, whose connectivity exhibits much slower changes [START_REF] Gong | Age-and gender-related differences in the cortical anatomical network[END_REF]. Sinke et al. [START_REF] Michel Rt Sinke | Bayesian exponential random graph modeling of whole-brain structural networks across lifespan[END_REF] studied DTI networks from 382 healthy subjects and considered 4 age-groups: 20 -34, 35 -50, 51 -70, > 70 years (Fig. 7). Assuming a low variability between the individuals in the same group, authors constructed a mean-GRN for each group and used a Bayesian ERGM (BERGM) [START_REF] Caimo | Bayesian inference for exponential random graph models[END_REF].

The H-function in this case was slightly different from Eq.21 i.e.

-H(G|θ)

= θ 1 x e (G) + θ 2 x gwesp (G|τ ) + θ 3 x gwnsp (G|ν) + θ 4 x nm (G) , (24)
where x nm (G) is a nodematch term (Eq.13) measuring in this case the number of connections within the same hemisphere.

The signs of the fitted model parameters (θ 1 < 0, θ 2 > 0, θ 3 < 0) confirmed the general tendencies observed in functional brain networks (Sec.4.1). The parameter θ 4 was instead negative for all groups, suggesting that the hemispheric nodematch is relatively low taking into account all the other included statistics. We stress that this does not necessary imply an absolute low intrahemispheric connectivity. Moreover, differently from previous studies [START_REF] Hagmann | White matter maturation reshapes structural connectivity in the late developing human brain[END_REF][START_REF] Willem M Otte | Aging alterations in whole-brain networks during adulthood mapped with the minimum spanning tree indices: the interplay of density, connectivity cost and life-time trajectory[END_REF], only weak changes were found between the parameter values across age-groups. Future studies will elucidate whether this result is actually reflecting a more stable topology of structural brain networks or a result of the group-averaging of the original brain networks.

Predicting states in temporal brain networks

So far, ERGMs have not explicitly taken into account the fact that brain networks are intrinsically dynamic and their topology can change, fluctuate and evolve over multiple time scales. Over short-time scales, brain networks rapidly reconfigure during motor or cognitive tasks to adapt to behavioral efforts or loads [START_REF] De | Persistent patterns of interconnection in time-varying cortical networks estimated from high-resolution eeg recordings in humans during a simple motor act[END_REF][START_REF] Valencia | Dynamic small-world behavior in functional brain networks unveiled by an event-related networks approach[END_REF][START_REF] Bassett | Dynamic reconfiguration of human brain networks during learning[END_REF][START_REF] Giulia Preti | The dynamic functional connectome: State-of-the-art and perspectives[END_REF]. Over longer time scales, brain networks can exhibit significant topological losses, such as in aging or in neurodegeneration [START_REF] Wu | A longitudinal study of structural brain network changes with normal aging[END_REF][START_REF] Dautricourt | Longitudinal changes in hippocampal network connectivity in alzheimer's disease[END_REF], as well as tentative restorative patterns in recovery after stroke or traumatic injuries [START_REF] Siegel | Reemergence of modular brain networks in stroke recovery[END_REF].

In these cases, it is crucial to formulate statistical network models that explicitly include time, so as to capture time-varying connection properties that are predictive of the system behaviour [START_REF] Holme | Temporal networks[END_REF][START_REF] Masuda | A guide to temporal networks[END_REF]. One possibility is to endow the ERGM with a time structure and estimate parameters from a time series, or sequence, of observed graphs G 1 , . . . , G T . One can start by assuming a first-order Markov time dependence, write the probability density function for the entire network sequence as

P (G t , . . . , G t-1 |θ) = T t=2 P (G t |G t-1 , θ) (25)
and generalize Eq.1 in a straightforward fashion:

P (G t |G t-1 , θ) = e θ•x(G t ,G t-1 ) G * ∈G e θ•x(G * ,G t-1 ) , (26)
where the x statistics can be now defined over two consecutive graphs in the timeseries. We shall refer to this model as temporal ERGM (TERGM) [START_REF] Hanneke | Discrete temporal models of social networks[END_REF].

TERGMs have been recently adopted to model brain network reorganization after stroke Obando et al. [START_REF] Obando | Temporal exponential random graph models of longitudinal brain networks after stroke[END_REF]. The study was carried out on longitudinal fMRI restingstate networks measured in 49 patients at 2 weeks, 3 months and 1 year after their first ever unilateral stroke event [START_REF] Siegel | Disruptions of network connectivity predict impairment in multiple behavioral domains after stroke[END_REF][START_REF] Siegel | Reemergence of modular brain networks in stroke recovery[END_REF] (Fig. 8a).

The model (Eq.26) was specified with the following temporal H function:

-H(G t |G t-1 , θ) = θ 1 x e (G t ) + θ 2 x stab (G t , G t-1 ) + + θ 3 x te (G t , G t-1 ) + θ 4 x ttr (G t , G t-1 ) (27)
Besides the standard edges term (Eq.11), the three other temporal statistics were defined as follows.

The stability term x stab simply counts the number of dyads that do not change in the time step t -1 → t and it has a fundamental role for the estimation since it sets the pace of allowed changes in the network sequence. The temporal-edge x te term counts the number of new inter-hemispheric edges appearing at time t, while the temporal-triangle x ttr term counts the number of triangle closures over time within the lesioned hemispheres (Fig. 8b).

These two last statistics are entrusted with the goal of mimicking the brain plasticity process i.e. the recovery of between-hemisphere integration and withinhemisphere segregation [START_REF] Siegel | Reemergence of modular brain networks in stroke recovery[END_REF][START_REF] Grefkes | Recovery from stroke: current concepts and future perspectives[END_REF]. The TERGM analysis resulted in positive estimates (θ 3 > 0, θ 4 > 0) and significantly higher in stroke patients compared to healthy controls (Fig. 8c). This suggested that both x te and x ttr are fundamental temporal processes of post-stroke brain reorganization. Notably, the same conclusions could not be reached when considering a static equivalent formulation of Eq.27 [START_REF] Obando | Temporal exponential random graph models of longitudinal brain networks after stroke[END_REF].

Moreover, the same parameters estimated from the temporal networks in the subacute phase (from 2 weeks to 3 months) significantly predicted the future language Illustration adapted from from [START_REF] Obando | Temporal exponential random graph models of longitudinal brain networks after stroke[END_REF].

score (1 year) of the patients, suggesting their potential use as clinical biomarkers (Fig. 8d). An alternative approach was adopted in [START_REF] Dell | Network analysis in disorders of consciousness: four problems and one proposed solution (exponential random graph models)[END_REF] to characterize the recovery from coma after traumatic brain injuries. Authors evaluated longitudinal fMRI brain networks obtained at 11, 18, 25 and 181 days after the trauma. In particular, they used separable temporal ERGMs (sTERGMs) [START_REF] Pavel | A separable model for dynamic networks[END_REF], which model separately and indipendently the tendencies to form and dissolve instances of different graph statistics.

The network at the time-step t is therefore given by the following:

G t = G + -(G t-1 -G -) (28)
where G t-1 it the graph time t -1, while G + and G -are the so-called formation and dissolution graphs. More in detail, the formation graph is drawn from ∼ P + (G t |G t-1 , θ + ), where P + is the same as in Eq.26 but with the formation statistics x + and parameters θ + . In this formation process, ties can only be added to G t-1 . Analogously, the dissolution graph is obtained by removing edges from G t-1 and the graph is drawn from ∼ P -(G t |G t-1 , θ -), with dissolution statistics x -and parameters θ -. The goal of a sTERMG is to estimate both sets of parameters θ + , θ -.

The approach in this study was to include in the model a huge number of potentially relevant statistics and let sTERMGs select the most relevant ones. Results showed indeed that only a subset of the estimated parameters turned out to be statistically significant. The recovery process was characterized by a tendency to preserve and strengthen connections within pairs of subsystems (e.g., default, frontoparietal, thalamus, visual) and between some of them (eg, default-visual, somatomotor-frontoparietal, ventral-visual). Similarly to the results found in [START_REF] Obando | Temporal exponential random graph models of longitudinal brain networks after stroke[END_REF], a significant propensity to close triangles over time also emerged from this study, thus confirming the importance of the local clustering connections in shaping the brain reorganization process from brain damage-induced coma.

On the interpretation of ERGMs

Up to now, we have not lingered on the foundation and subtleties of ERGMs, albeit they are crucial to understand and interpret the related results. Due to the number and variety of scientists and practitioners that have developed or used these models, there is sometimes confusion in the scientific literature of ERGMs, if not in the notions at least in the lexicon. This might lead to unsupported claims or reckless interpretations of the results. We therefore aim at clarifying in this section the philosophy of the method and its theoretical underpinnings.

Our discussion lies on three fundamental concepts: i) description, ii) prediction, and iii) explanation. We briefly pin-down for clarity the corresponding definitions. Description means a characterization, or a summary, of the state of the system under investigation, based on the choice of suitable variables and observables to be measured. By definition, a description refers to the past or present state of the system. Prediction is the anticipation (forecast) of the future state of a system. Explanation aims at reaching a causal or mechanistic understanding of why a particular phenomenon or event has occurred, based on a theoretical framework or empirical evidence [START_REF] Carl | Studies in the logic of explanation[END_REF][START_REF] Joseph | Explanation, prediction, description, and information theory[END_REF][START_REF] Shmueli | To explain or to predict[END_REF].

We argue that in general ERGMs as those described in (Eq.1) and their derivatives, are certainly descriptive, might be predictive, but they are not explicative. Note that a method can be tremendously efficient in predicting outcomes of an experiment while at the same time providing no explanation on the underlying physics. Analogously, a description of the phenomena under investigation is radically different from a scientific explanation [START_REF] Carl | Studies in the logic of explanation[END_REF].

ERGMs and the maximum entropy principle (MEP)

First, let us explicit the origin of the ERGM formulation from basic information theoretic principles. ERGMs can be indeed seen as one of the many particular instances belonging to the broad class of Maximum Entropy (Max.Ent.) models [START_REF] Edwin | Information theory and statistical mechanics[END_REF][START_REF] Cimini | The statistical physics of real-world networks[END_REF][START_REF] Thomas | Elements of information theory[END_REF]. Similar models have been applied to neural populations [START_REF] Schneidman | Weak pairwise correlations imply strongly correlated network states in a neural population[END_REF], antibody diversity [START_REF] Mora | Maximum entropy models for antibody diversity[END_REF], protein sequences [START_REF] Cocco | Inverse statistical physics of protein sequences: a key issues review[END_REF], economics [START_REF] Sanjay | What is an explanation? statistical physics and economics[END_REF], population genetics [START_REF] Dichio | Statistical genetics in and out of quasi-linkage equilibrium[END_REF] and many more. Let us then derive the ERGMs from the MEP.∥.

∥ Note, however, that this is not the way the ERGMs were originally introduced, which instead was based on the Hammersley-Clifford theorem for Markov graphs [START_REF] Frank | Markov Graphs[END_REF], building on a previous work of J. Besag in the context of spatial models of lattice systems [START_REF] Besag | Spatial Interaction and the Statistical Analysis of Lattice Systems[END_REF].

Consider an observed graph G * ∈ G and assume that all the relevant information about G * is encoded in a set of sufficient statistics x(G * ). According to the MEP, the most unbiased probability density function P (G) consistent with available knowledge is obtained by maximizing the Shannon information entropy S S[P ] = -G∈G P (G) log P (G) [START_REF] Peter | In vivo fiber tractography using dt-mri data[END_REF] under the constraints

G∈G P (G) = 1 , (30)
G∈G x(G)P (G) = x(G * ) : (31)
the first one is just a normalization, the others fix the expected values of the sufficient statistics over the ensamble G to the correspondent values computed for G * . The rationale at the basis of the Max.Ent. recipe is very pragmatic, as already stressed in [START_REF] Edwin | On the rationale of maximum-entropy methods[END_REF], and dates back in similar forms to Bernoulli and Laplace. Let us consider the total number of ways in which the probability distribution of a system can be realized, given that all its possible states are equivalently likely. By a simple combinatorial theorem [START_REF] Edwin | On the rationale of maximum-entropy methods[END_REF], it can be shown that the Max.Ent. distribution is the most likely and also that any other distribution becomes highly atypical as the number of degrees of freedom -in our case, the number of possible edges -becomes large.

In order to solve the constrained maximization problem defined by Eqs.29-31 we introduce a set of Lagrange multipliers λ ∈ R r+1 and maximizing over P the expression: [START_REF] Marcus E Raichle | Behind the scenes of functional brain imaging: a historical and physiological perspective[END_REF] which gives P (G|λ) = e λ0-1+ r ι=1 λιxι(G) . By imposing (Eq.30-31), we find λ 0 = 1 -log Z and λ ι = θ ι for ι = 1, . . . , R, with the θ being defined as the choice of the corresponding multipliers which satisfy Eq.31. The result is precisely Eq.1.

- G∈G P (G) log P (G)+λ 0 G∈G P (G)-1 + R ι=1 λ ι G∈G x ι (G)P (G)-x ι (G * ) ,
The same argument can be used for a temporal series of graphs G = G 1 , . . . , G T . If the information is available in the form of graph statistics x(G * ) evaluated on the observed sequence G * , the path entropy reads as

S[P ] = G 1 ∈G • • • G T ∈G P (G) log P (G) , (33)
and the mass probability function becomes

P (G|θ) = e θ•x(G) G1 ∈G • • • GT ∈G e θ•x( G) , (34)
which is the more general case of the TERGMs presented before (Eq. [START_REF] Fornito | The connectomics of brain disorders[END_REF][START_REF] Jones | White matter integrity, fiber count, and other fallacies: the do's and don'ts of diffusion mri[END_REF]. This generalized variational principle is sometimes referred as maximum caliber (Max.Cal.) [START_REF] Pressé | Principles of maximum entropy and maximum caliber in statistical physics[END_REF].

Description and prediction: the choice of graph statistics

As discussed in (Sec.3.3.1), when formulating ERGMs it is crucial to avoid degeneracy issues which might occur especially with too simplistic model formulations. The possibility of including several potentially relevant effects (model terms) has two advantages. On the one hand, it helps alleviating degeneracy issues, on the other hand it allows to statistically test the joint effect of multiple network properties. Once the model is specified and as long as the estimation process converges properly, an ERGM is always descriptive, the estimated parameters information on the relevance of the associated statistics. Whether an ERGM is also predictive is a less trivial question that we address in the following.

A delicate implicit assumption of ERGMs is that the selected statistics x(G * ) encode all the relevant information of the system, i.e. they are sufficient statistics. This assumption raised major criticism around Max.Ent. and related methods, see e.g.

[132, [START_REF] Auletta | On the relevance of the maximum entropy principle in non-equilibrium statistical mechanics[END_REF][START_REF] Dichio | Statistical genetics in and out of quasi-linkage equilibrium[END_REF]. Indeed, what is deemed to be relevant entirely relies on the modeler's hypothesis and is consequently affected by some arbitrariness. This is particularly true in complex systems science where, almost by definition, we rarely have control on the hidden degrees of freedom [START_REF] Bassett | On the nature and use of models in network neuroscience[END_REF][START_REF] Ladyman | What is a complex system?[END_REF]. It might be that what we aim to measure, or what we are able to measure, simply is not a relevant feature of the system.

In the case of TERGMs, such as those in Eq.25-26, there are two more subtleties. First, the Markov hypothesis in Eq.25 is in general hard to justify a-priori [START_REF] Onsager | Fluctuations and irreversible processes[END_REF][START_REF] Auletta | On the relevance of the maximum entropy principle in non-equilibrium statistical mechanics[END_REF]. Second, when writing Eq.26 we implicitly assume that the selected statistics are relevant throughout the observed time window. This is true if we can assume the stationarity of the process [START_REF] Gardiner | Handbook of Stochastic Methods for Physics, Chemistry and the Natural Sciences[END_REF], whose biological plausibility is however disputable [START_REF] Fang | Nonequilibrium physics in biology[END_REF]. Both these two assumptions are further complicated by the fact that in many real applications the measurements of networks are not equally spaced in time [START_REF] Block | Change we can believe in: Comparing longitudinal network models on consistency, interpretability and predictive power[END_REF].

To address these issues, a-posteriori model-selection approaches like goodness of fit (GoF) or out-of-sample assessment can be adopted to rule out models that are markedly wrong. These methods evaluate the ability of the graphs simulated with a fitted ERGM to reproduce properties of the system other than those explicitly included. When this happens then the estimated ERGM is not only descriptive but also predictive. Thus, when using ERGMs with GOF assessments, it is not only important to use appropriate relevant graph statistics but also to have a precise hypothesis on the network features against which to test the predictive performance.

Explanation: ERGMs and Stat.Mech.

In the previous sections, we have insisted on the role of the modeler in shaping the mass probability distribution by selecting a set of graph statistics. This is consistent with the E.T. Jaynes' subjective interpretation of probability [START_REF] Edwin | Probability theory: The logic of science[END_REF][START_REF] Pressé | Principles of maximum entropy and maximum caliber in statistical physics[END_REF], according to which probabilities are epistemic statements about physical systems. Put differently, they reflect our state of knowledge or, equivalently, of ignorance, about the world.

Such subjective probability is enough for purely inferential purposes. For example, even if Eq.1 is not the true probability density function of the system under investigation, it may nevertheless have good predictive performance (GoF) for a given set of test statistics. Indeed, the ERGM inferential scheme is agnostic about the data-generating process.

This property is somewhat unsatisfactory if the ultimate goal is to explain the physics underlying the observed data. In that case, probabilities should not reflect a subjective state of knowledge but they should have an objective physical meaning. This is, for instance, the case of statistical mechanics where probabilities have an objective irreducible meaning, they do not depend on the observer. The role of the theory is then to explain how the macroscopic properties of bodies, mainly thermodynamic properties, statistically arise from their microscopic structure [START_REF] Peliti | Statistical mechanics in a nutshell[END_REF]. The formal analogy between the ERGM in (Eq.1) and the Gibbs-Boltzmann distribution for physical systems at equilibrium has instilled hopes for a similar microlevel objective interpretation of the ERGM probability density function. However, on a closer inspection we immediately realize that this equivalence is problematic.

These two exponential distributions have in fact very different theoretical underpinnings. Just as an example: at the very foundation of statistical mechanics, and hence of Gibbs-Boltzmann exponential distributions, there is the assumption that the system dynamics obeys a detailed balance [START_REF] Boltzmann | Weitere studien über das wärmegleichgewicht unter gasmolekülen[END_REF]. In its simplest form, detailed balance means that for any two accessible states r, s the transition rate W r→s is equal to its inverse W s→r . This assumption becomes highly disputable when it comes to biological processes. Detailed balance is indeed violated at the molecular level and apparently so also for large-scale brain dynamics even during intrinsic spontaneous activity [START_REF] Christopher W Lynn | Broken detailed balance and entropy production in the human brain[END_REF]. Therefore, in general, that between the ERGM and the Gibbs-Boltzmann distributions should not be considered more than a formal analogy, and their respective interpretations should not be confused.

Conclusion and perspectives

Statistical models are crucial to gain intuition on the connection rules and generative mechanisms of brain networks across multiple time and spatial scales [START_REF] Presigny | Colloquium: Multiscale modeling of brain network organization[END_REF]. Furthermore, they can be used to improve the detection of the underlying network properties in the presence of statistical noise in the data [START_REF] Efron | Bootstrap methods for standard errors, confidence intervals, and other measures of statistical accuracy[END_REF][START_REF] Gfeller | Finding instabilities in the community structure of complex networks[END_REF][START_REF] Karrer | Robustness of community structure in networks[END_REF].

Here, we have discussed a family of models -exponential random graph models (ERGMs) -which allows for simultaneously taking into account different local connection properties, disentangle their interactions and quantify their statistical influence on the global observed network. The number of studies using ERGMs to model brain connectivity networks has rapidly increased in the last decade. We have discussed how they have been used to characterize intrinsic brain connectivity, build group representative networks, discriminate different brain states, as well as model temporal patterns, with focused applications to brain diseases.

In addition, we have placed the ERGM framework in the broader context of maximum entropy models and discussed in detail the subtleties of network analyses based on the Max.Ent. recipe, with particular focus on the meaning of the inferred model parameters. We have discussed the descriptive and predictive power of ERGMs, and warned about its interpretation as an equilibrium theory of the brains' biological processes, for which a non-equilibrium picture is more appropriate.

Despite the increasing number of works using or developing ERGMs to address neuroscience-related questions, the field is still in its infancy and we hope this review will encourage curious physicists, network scientists, and neuroscientists to pursue the research direction we have laid out here. In the years to come, we anticipate a boost in the field fueled by an increasing availability of high-quality experimental data, which is crucial to validate the hypothesized mechanisms. This includes the use of ERGMs to model brain diseases not only from a recovery perspective (e.g., stroke or coma), but also in terms of degeneration (e.g. Parkinson, Alzheimer's), as well as to evaluate the efficacy of treatments for psychiatric conditions (e.g. Schizofrenia) and developmental disorders (e.g. autism).

On a more methodological note, only a small fraction of the ERGM-like galaxy of statistical models has been employed so far in network neuroscience. Among the unexplored alternatives, we spot out the potential use of multilayer ERGMs [START_REF] Pavel N Krivitsky | Exponential-family random graph models for multi-layer networks[END_REF] to capture interactions between different scales of neural organization [START_REF] Presigny | Colloquium: Multiscale modeling of brain network organization[END_REF]. As for temporal analyses, a promising family of models beyond the (T)ERGM paradigm is the related stochastic actor oriented model (SAOM), which simulates time-varying network changes through a microscale process defined at the level of the nodes, i.e. the actors [START_REF] Tom Ab Snijders | Introduction to stochastic actor-based models for network dynamics[END_REF]. Critically, both TERGMS and SAOM only take into account monotonic network changes increasing or decreasing over time. However, modeling brain reorganization might require capturing non-monotonic mechanisms, too. This is particularly true for intrinsic brain functioning and cognitive/motor tasks, which typically exhibit oscillatory time-varying dynamics [START_REF] Giulia Preti | The dynamic functional connectome: State-of-the-art and perspectives[END_REF]. In this direction, the development of models that integrate this dynamic component, such as the varyingcoefficient ERGM [START_REF] Lee | Varying-coefficient models for dynamic networks[END_REF], appear particularly important to overcome such limitations. A distinct mention is for stochastic block models (SBMs), which allows flexible and scalable analysis of complex networks [START_REF] Paul W Holland | Stochastic blockmodels: First steps[END_REF]. In addition to handling community structures, recent SBM extensions also include several types of structures, such as node degree sequences [START_REF] Karrer | Stochastic blockmodels and community structure in networks[END_REF], hierarchical structures [START_REF] Tiago | Hierarchical block structures and high-resolution model selection in large networks[END_REF], edge weights [START_REF] Tiago | Nonparametric weighted stochastic block models[END_REF], and triadic closures [START_REF] Tiago | Disentangling homophily, community structure, and triadic closure in networks[END_REF]. SBM applications to brain networks are still in their infancy (see [START_REF] Faskowitz | Weighted stochastic block models of the human connectome across the life span[END_REF] for a recent example), but they surely represent a future promising research avenue. Finally, a long-term and ambitious goal would be to enrich statistical models by including higher-order interactions, which might reveal previously unappreciated brain organization principles [START_REF] Battiston | Networks beyond pairwise interactions: structure and dynamics[END_REF].

We believe that combining physics-inspired concepts, biological knowledge, and statistical methods is essential for unraveling the intricacies of the human brain's structure and dynamics. This synergy holds the potential to provide crucial insights into the complexity of the brain.

Figure 1 :

 1 Figure 1: Statistical analysis of brain networks inferred from experimental data. (a) Neuroimaging techniques provide images (e.g., diffusion tensor imaging) and signals (e.g., electroencephalography) informing about the anatomical structure and functional dynamics of the brain, respectively. b) Methods from image and signal processing are then used to build the brain network G * , where nodes correspond to different brain regions and links correspond to structural connections or functional interactions. c) The aim of statistical modeling is to describe P (G), the probability distribution in the space of all possible graphs G. It encodes the architectural principles of the observed brain network G * . Typically, the latter has intermediate properties between a totally random G random and a totally ordered G ordered network.

Figure 2 :

 2 Figure 2: Illustration of graph statistics for ERGMs. a) x e counts the total number of edges in the graph; b) x (k) d counts the number of nodes with degree k; c) x (k) nsp is the number of non-connected dyads sharing k common partners; finally, d) x (k) esp counts the connected dyads sharing exactly k common partners.

Figure 3 :

 3 Figure 3: ERGMs workflow for brain networks. a) A graph G * is given, for instance a structural or functional brain network. b) The information stored in the graph G * is condensed in a set of statistics x(•) ∈ R r , whose choice is up to the modeller. The statistics evaluated on G * represent the input of the ERGM inference methods. c) The goal of an ERGM is the estimation of the set of parameters θ * ∈ R r associated with x(•). d) The inferred parameters can be used to simulate n new brain networks G 1 , . . . , G n that statistically reproduce the observed network i.e. 1 n

Figure 4 :

 4 Figure 4: Differences in ERGM parameters across brain subsystems. Data from[START_REF] Paul E Stillman | A consistent organizational structure across multiple functional subnetworks of the human brain[END_REF]. For each subsystem and parameter, the mean value and st.dev. are reported. Asterisks indicate that the estimation was significantly different from 0 for all of the 21 subjects involved in the study. In general, larger networks (cingulo-opercular, salience, default mode, fronto-parietal) yield more robust estimations i.e. significant for a larger number of subjects. The corresponding brain subsystems are shown along the x-axis for illustrative purposes, as defined in[START_REF] Gu | Controllability of structural brain networks[END_REF].

Figure 5 :

 5 Figure 5: Comparison between different GRN construction schemes. Functional (fMRI) brain networks from 100 subjects aged 18 -33. Gray histograms show the distributions of parameter estimates from a Bayesian ERGM (BERGM) across individuals. Vertical lines correspond to group-level parameters. In red, the multi-BERGM (red line). Blue dotted line for mean-BERGM (group-averaged BERGM individual estimates). Green dashed line for a single BERGM estimated on the mean-GRN ; purple dashed line for a BERGM on the median-GRN. Both multi-BERGM and mean-BERGM show greater accuracy than simpler mean/median-GRN in capturing the correct distribution of parameters in the population. Illustration from [95].

Figure 6 :

 6 Figure 6: Comparison of multi-BERGM parameters between young and old individuals. Two groups of subjects are considered: young 18 -33 (Y) and old 74 -89 (O). 100 functional (fMRI) brain networks for each group. The young group is characterized by larger values for the GWNSP parameter while the edges and GWESP posteriors were almost identical between the two. Illustration from [95].

Figure 7 :

 7 Figure 7: Structural brain network changes across lifespan. Top left: mean-GRN for the age range of 20 -34. Other networks correspond to brain networks obtained in distinct older groups of individuals. Red lines represent "lost connections" and green lines represent "new connections". Brain network density changes across the lifespan with +1.6% (age 35-50), -4.4% (age 51-70) and -8.0% (age > 70). Illustration from [90].

Figure 8 :

 8 Figure 8: TERGM analysis of longitudinal brain networks after a stroke event. a) Dynamic brain network of a stroke patient 2 weeks, 3 months and 1 year after the stroke event. Different colours emphasize the belonging to different brain systems. b) Graphical representation of temporal edges x te and temporal triangles x ttr . c) Statistical comparison between θ 3 , θ 4 for stroke patients (red shape) and healthy controls (white shape). Violin plots show the distribution of the values, while box-plots represent median and quartiles. Significant increases in the stroke group are reported for both the parameters considered (p = 0.045). d) Correlation of the TERGM coefficients and future language outcome of stroke patients. Both temporal edges and triangles in the subacute phase (2 weeks, 3 months) are significantly associated with the future language score in the chronic phase (1 year) (p < 0.03).Illustration adapted from from[START_REF] Obando | Temporal exponential random graph models of longitudinal brain networks after stroke[END_REF].

  .3).

	Reference	Data (nodes)	Method/Package
	[88] Simpson et al. (2011)	fMRI(90)	ERGM / ergm [60]
	[89] Simpson et al. (2012)	fMRI(90)	ERGM / ergm [60]
	[90] Sinke et al. (2016)	DTI(96)	BERGM / Bergm [68]
	[91] Obando et al. (2017)	EEG(56)	ERGM / ergm [60]
	[92] Stillman et al. (2017)	fMRI(20)	cGERGM/ gergm [69]
	[93] Dell'Italia et al. (2018)	fMRI(148)	sTERGM / stergm [72]
	[94] Stillman et al. (2019)	fMRI(9 -25)	cGERGM/ gergm [69]
	[95] Lehmann et al. (2021 )	fMRI(90)	BERGM / Bergm [68]
	[96] Obando et al. (2022)	fMRI(81.13

* ) TERGM / btergm

[START_REF] Leifeld | Temporal exponential random graph models with btergm: Estimation and bootstrap confidence intervals[END_REF]

Acknowledgements

We would like to thank the anonymous reviewers, whose comments helped us improving the initial version of the review. We thank Erik Aurell and Mario Chavez for many insightful discussions and for providing useful remarks on the manuscript. We are also thankful to Thibault Rolland and Remy Ben Messaoud for their contribution in the preparation of the figures. FDVF acknowledges support from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program (Grant Agreement No. 864729)