Vito Dichio

Fabrizio De

Vico Fallani

The exploration-exploitation paradigm for networked biological systems

The stochastic exploration of the configuration space and the exploitation of functional states underlie many biological processes. The evolutionary dynamics stands out as a remarkable example. Here, we introduce a novel formalism that mimics evolution and encodes a general explorationexploitation dynamics for biological networks. We apply it to the brain wiring problem, focusing on the maturation of that of the nematode C. elegans. We demonstrate that a parsimonious maxent description of the adult brain combined with our framework is able to track down the entire developmental trajectory.

Introduction. Modelling and analyzing the dynamics of biological systems is notoriously challenging. Critically, they are often stochastic in nature as they involve and possibly exploit some degree of randomness. At the same time, biological dynamics are also shaped by functional constraints that determine which outcomes are viable. The constraints emerge from the need for biological systems to perform specific tasks, and act on the system as a whole, not on specific components. Stochastic events that violate these constraints are unlikely to persist, while those that align with them are more likely to become integrated. If the details of such explorationexploitation (EE) dynamics are context-dependent, general principles can still be formulated [START_REF] Bialek | Biophysics: searching for principles[END_REF]. This entails addressing several challenging questions. For instance, how do biological systems explore the space of possible configurations? How do they identify the optimal states that satisfy specific functional demands?

A possible solution is offered by Nature itself. In evolutionary dynamics, a population primarily evolves under the combined action of mutations and recombinations (exploration) and natural selection (exploitation). The latter is based on the notion of fitness: those individuals that are more apt to the environment will have a higher reproductive success (high fitness) and survive to the next generations, while the others will go extinct, Supplemental material Section I (SM-I) [START_REF] Crow | An introduction to population genetics theory[END_REF]. We argue that evolutionary dynamics is a particular instance of the aforementioned EE dynamics and build upon it to construct a general EE formalism for networked biological systems.

We use it to tackle the brain wiring problem and model the developmental dynamics of the C. elegans connectome, recently obtained by serial-section electron microscopy [START_REF] Dichio | The exploration-exploitation paradigm: a biophysical approach[END_REF][START_REF] Hassan | [END_REF][5].

Theoretical framework. Let us begin by clarifying the terminology. (a) Exploration refers to the act of stochastically searching the configuration space. (b) Exploitation refers to the harnessing the discovered configurations to optimise the system function. The resulting optimisation problem is defined once we specify (b.i) how the optimal states are encoded and (b.ii) how the system approaches them.

Formally, let us consider a biological system represented as a simple graph (or network) G ∈ G over N nodes, unweighted, undirected, with no self loops. It can be identified with a finite, binary, symmetric and with zero-diagonal adjacency matrix G = {a ij }, where a ij ∈ {0, 1} indicates the absence or presence of an edge within the pair of nodes, or dyad, (ij). There are L = N (N -1)/2 dyads, hence L possible edges. Let P (G, t) be the probability of the graph G at time t.

(a) Exploration. Each dyad mutates its state in the time interval ∆t with rate µ ≥ 0. A simple exploration scheme is to randomly create or dissolve edges, e.g., an edge is added if none existed or removed if present. The effect on the graph distribution is P (G, t+∆t) = P (G, t) + ∆tµ i<j [P (M ij G, t) -P (G, t)] , [START_REF] Bialek | Biophysics: searching for principles[END_REF] where M ij is the operator that mutates the dyad a ij of the graph G. The exploration rate µ is here constant and uniform across dyads.

(b) Exploitation. A functional metric F (G) : G → R serves the purpose of representing the concept of biological function, with optimal states defined as maxima of F (b.i). In the time interval ∆t, we formally define exploitation as follows:

P (G, t + ∆t) = e ∆tφF (G) ⟨e ∆tφF ⟩ t P (G, t) , (2)
where ⟨•⟩ t stands for the ensemble average at time t i.e. ⟨e ∆tφF ⟩ t = G e ∆tφF (G) P (G, t) and φ ≥ 0, the exploitation rate, is an overall scaling. Therefore, the way in which the dynamics approach the most functional (highest F ) configurations is by exponentially increasing the probability of those graphs that have higher F values than the ensemble average at time t (b.ii). We will refer to the ratio ρ = φ/µ as the functional pressure: ρ ∼ 0 implies a dynamic dominated by randomness, similar to a random walk in the graph space G, while ρ → ∞ corresponds to the limit of perfect exploitation, where only the most functional graph configurations have non-negligible probabilities.

In SM-II we introduce four simple models, namely the cases of no exploitation, edge penalty, edge covariate and distance-like F metric. We show that these cases can be treated analytically and offer a formal understanding of the intuitions (i) that the optimal states implied by a F metric are not strictly attainable as long as µ ̸ = 0, and (ii) that the functional pressure ρ controls not only the rate of approach to the maxima of F but also the final stationary state.

The above framework closely mimics a Darwinian evolution driven by mutations and fitness-based natural selection, in the infinite population limit [6,7]. Eq.(1-2) can be regarded as an algorithm, inspired by evolution, that (a) uses random choices to (b) direct an exploitative search for solving an optimisation problem. In this sense, it is similar to a genetic algorithm [START_REF] Goldberg | Genetic algorithms in search, optimization, and machine learning[END_REF], see SM-I. We further explore the parallel with the evolutionary dynamics to design simulations based on eq.(1-2), concisely described in SM-III.

F is shaped by the environment and has the role of mapping the functional requirements of the biological system onto the configuration space. As a consequence, the particular form of F depends on the specific context and system. In general, there are no requirements on the properties of F , which could be regarded as a black box which returns a real number for each possible input (graph). In this work, however, we will study a whitebox F metric which admits a mathematical formulation. In particular, we will describe the state of a graph G by a set of sufficient statistics x(G) ∈ R r and in this latter space the F : R r → R will be formally defined.

C. elegans brain maturation. In the following, we will consider the so-called brain wiring problem [START_REF] Hassan | [END_REF], i.e., how the structural complexity of a natural brain arises during the development of an organism. Answering the brain wiring problem is an open challenge in neuroscience and entails tackling at least two kind of questions: (i) what drives the brain maturation (structural principles) (ii) which is the driving algorithm (dynamical principles).

Here, we will formulate them in terms of a F -metric and EE dynamics, respectively. In particular, the latter is consistent with three essential and general features of the brain wiring dynamics, which are (a) functionally robust -the adult brains are capable of supporting the functions that sustain the life of an organism;

(b) not hardwired -genetically encoded developmental algorithms give rise to similar yet non-identical structures, resulting in high inter-individual variability [START_REF] Hassan | [END_REF];

(c) self-referential -the updating rules evolve in time, as a function of the state and therefore of the history of the system [START_REF] Goldenfeld | [END_REF].

To tackle the brain wiring problem, a natural choice is to consider that of the nematode C. elegans [10], SM-V.A. This is the only organism for which a comprehensive map of neuronal connections within a brain has been reconstructed across development [5]. The dataset consists of 8 fully reconstructed brains of the hermaphrodite C. elegans, obtained from different isogenic individuals at different developmental ages, including one at birth and two adults (t ∼ 45 h after birth), SM-V.B. We consider the unweighted and undirected networks of chemical synapses between sensory, inter, motor and modulatory neurons (161 -180 nodes, 617 -1669 edges). This choice of representation is motivated by the statistical properties of the adult C. elegans connectome [11,12], along with the effort to devise a simplified growth model, a critical discussion can be found in SM-V.B.

The developmental principles that guide the C.elegans brain maturation are not entirely known. On the one hand, approximately 43% of the synaptic connections between neurons are not conserved among genetically identical individuals, suggesting a prominent role of stochasticity in the brain wiring [START_REF] Hassan | [END_REF]5]. Conversely, the diverse range of behaviors exhibited by adult C.elegans [13] demands functional selection. Our EE framework captures these two tendencies simultaneously. An overview of the approach is illustrated in fig. [START_REF] Bialek | Biophysics: searching for principles[END_REF]. A micro-level interpretation of eq.(1-2) for the wiring dynamics of the individual neurons is extensively discussed in SM-V.C.

A preliminary step of our modelling approach is the characterization of the worm brain by a set of sufficient statistics x(G). Based on recent evidence [14,15], we consider a parsimonious representation in which:

x(G) = k>0 w (k) τ d x (k) d (G) k>0 w (k) τesp x (k) esp (G) , (3)
where

x (k) d and x (k)
esp are the number of nodes with degree k and the number of connected dyads sharing exactly k partners, respectively. The coefficients are w

(k) α = e α {1 -1 -e -α k }, with α = τ d , τ esp > 0 decay param- eters.
In other words, these statistics are linear combinations of the degree and edgewise shared partner distributions. They yield a model that is both realistic and computationally tractable [16], SM-IV. More specifically, the first statistic is called geometrically weighted degree (gwd) and encodes the information, e.g., on the presence/absence of hub-nodes in the graph, well documented in the case of the C.elegans [5,14]. The second statistic is called geometrically weighted edgewise shared partner (gwesp) and is a proxy for a triadic-closure phenomenon in the graph, i.e., pairs of nodes that have links to one or more common neighbours have a higher chance of being connected to each other. The latter could in turn result from a tendency of the network to segregate into densely connected modules. [5,14].

Given a choice of statistics as in eq.( 3), we can characterize any observed graph G * within the inferential framework of exponential random graph (ERG) models [16,17], SM-IV. Accordingly, the ensemble G is endowed with a maximum entropy probability distribution

P ERGM (G|θ) = e -H(G,θ) / G∈G e -H( G,θ) (4)
where

H(G, θ) = -θ • x(G)
is the Hamiltonian. Given an observed graph G * , the vector of parameters θ * can be inferred as approximate solution the maximum likelihood estimation problem θ * = arg max θ log P (G * |θ). The inferred parameters quantify the contribution of the associated statistics to the structure of the observed graph. For example, if the statistic x α is non-negative, positive θ * α imply the existence of a bias towards graphs with higherthan-random values of x α -given the rest of the model β̸ =α θ * β x β [16], SM-IV. We are now in the position to propose the following F metric for the C. elegans brain maturation

F (G) = θ * • x(G) , (5)
where x(G) are defined in eq.( 3). The parameters θ * are obtained from the ERG inference in the adult stage, so that the correct (functional) balance of model statistics can be achieved at the end of the developmental process. The topography of the functional landscape results from the combined effect of physical, genetic and functional constraints.

In particular, we consider the average estimated parameters from the two adult worms

G * T = (G * T,1 , G * T,2
) and obtain θ * gwd = 0.44, θ * gwesp = 0.58, SM-VI.A. EE dynamics based on eq.( 5) will favor both the emergence of hubs and of a triadic closure behaviour since, by virtue of the positive values of the linear parameters, higher values of the statistics in eq.( 3) will imply higher F values. This is in line with experimental observations that, during development, hub neurons at birth get more inputs and that the overall modularity of the C. elegansbrain network increases [5].

We can now proceed to model the developmental dynamics by setting appropriate boundary conditions, and the EE parameters of the dynamics: the exploration rate and the functional pressure. As the argument goes, µ and ρ are characteristic of the specific instance (C. elegans) of the biological process (brain wiring), they are genetically encoded and therefore result from the evolutionary history of the species.

We set the graph G * 0 corresponding to the network at birth as the starting point of the dynamics, P (G = G * 0 , 0) = 1. In fact, (i) as reported in [5], the brain morphology at birth serves as the structural foundation upon which the adult connectivity unfolds. Moreover, (ii) an implicit assumption of the EE graph dynamics is the functional homogeneity, i.e., that the same F metric holds true throughout the whole dynamics. This assumption is likely to be violated before hatching (birth), during the embryonic stage, where a different growth regime of the nervous system has been observed [18].

Throughout development, the removal of synaptic connections happens rarely [5]. Accordingly, we modify the mutation scheme described in eq.( 1) by restricting the removal of edges. By Occam's razor, we assume a constant exploration rate

µ * = 1 T L i<j āij (G * T )-a ij (G * 0 ) = 1.43×10 -3 h -1 , (6)
where T = 45h is the adult age, L = N (N -1)/2 is the number of dyads of the adult brain graphs (N = 180), i<j āij (G * T ) is the average number of edges between the two adult worms and i<j a ij (G * 0 ) is the number of edges at birth. At the end of the developmental dynamics, µ * is assumed to drop to zero.

We are left with only one free parameter, i.e., the functional pressure ρ = φ/µ > 0. Nonzero ρ describe the scenario in which new connections emerge primarily in locations where they result in an enhanced system function. Its value must be biologically regulated to ensure the development of adequately specialized functional circuits prior to reaching the adult stage. Therefore, we use the corresponding degree of freedom to inform the EE graph dynamics about the age of the adulthood. In particular, we set ρ * = min ρ δ mah T , where the quantity to be minimized is the Mahalanobis distance [19], at time T , between the two-dimensional ensemble distribution of the graph statistics and the average experimental values, SM-VI.B.

In fig.

(2a) we show that that δ mah T is a convex function of the functional pressure ρ. Both insufficient and excessive ρ lead to the ensemble distribution diverging from the experimental values. The minimization procedure yields ρ * = 9.017 × 10 2 (R 2 = .98).

Notably, in estimating the dynamical parameters we have relied solely on the C. elegans brain graphs at birth and in the adult stage. We can ask if and how the estimation would change when considering the whole available data, which include the developmental time-points at 5 h, 8 h, 16 h, 23 h, 27 h and 45 h after birth. A linear fit of the growth of the number of edges based on the whole time series yields µ * * = (1.39 ± 0.08) × 10 -3 h -1 . This estimation is compatible with the value µ * in eq.( 6), fig. (2b). As for the functional pressure, we can define an equivalent minimization problem where the Mahalanobis distance is summed over all experimental time points, yielding a value ρ * * that is statistically consistent with ρ * , fig.(2a), SM-VI.B.

This hints that (i) given the model eq.( 5), the assumption of functional homogeneity for the worm brain wiring dynamics holds true and (ii) the EE graph dynamics, informed about birth and adulthood, capture the entire developmental trajectory.

To further investigate this result, we can fix µ * , ρ * and look at individual simulations of the brain growth, as in fig. (2c). For comparison, we also plot the results of a null model where ρ = 0, i.e., a random graph growth with exploration rate µ * . As could be expected, the adult stage is correctly reached in terms of the model statistics eq.( 3) and, by consequence, of the F -metric. Notably, however, our simulations approximate en passant the other observed developmental ages, which we have used nowhere in inferring the parameters. This opens up the possibility of using our framework to make predictions also for stages of brain maturation for which no data are available.

The model described so far is a simple, low-dimensional model of the underlying biological dynamics. Almost by construction, a choice as simple as eq.( 3) is unlikely to capture the finer-scale topological details of the observed graphs. Analogously, an exploration dynamics as simple as eq.( 6) cannot capture transient dynamic patterns. Yet, we can meaningfully ask to what extent the EE graph dynamics based on the features eq.( 3) reproduces other network properties not included in the model formulation, SM-VI.C. In fig. (2d) we show that our model retrieves the propensity of the C. elegans brain networks to exhibit relatively high efficiency (like random graphs) and clustering (unlike random graphs) [START_REF] Latora | Complex networks: principles, methods and applications[END_REF].

Discussion. To summarize, we have presented a parsimonious, interpretable framework for the dynamics of networked biological systems. It is built upon the dynamical principle of the exploration-exploitation (EE) paradigm, which is general. It serves as theoretical scaffolding for formulating specific dynamical models, which must be tailored to the biological system. We have used it here to model the growth of the C. elegans connectome, from birth to adulthood. Notably, our results suggest that the knowledge of the birth and adult age is sufficient for the EE graph dynamics to describe the whole developmental trajectory. We speculate that the same may be true for the connectomes of other living systems [START_REF] Hildebrand | [END_REF][22][23], for which no such data as the developmental trajectory are available to date. This hypothesis is poised for experimental validation in the near future.

Our model should be regarded as a first step towards a more detailed understanding of the brain maturation. To this end, the framework here presented supports straightforward extensions to more complex exploration schemes, accounting for non-uniform synapse addition, directed flow of synaptic information, neuron-specific information, homophily effects, and physical or functional constraints [5,14]. A detailed discussion of the possible model extensions can be found in SM-VII.A. Beyond structural connectivity, it would be interesting to study under the same lens the C. elegans brain functional connectivity, recently mapped for the adult stage [24], where there exists a closer correlation between the notion of biological function and the topology of the graph.

Zooming out, our framework can be broadly used to study the dynamics of complex systems arising from the interplay between (i) the variability fueled by a stochastic search of the configuration space and (ii) the statedependent optimisation of an objective function -we propose three examples in SM-VII.B. Importantly, as showcased here, this can be done by introducing only a very limited number of interpretable parameters [25] [START_REF]Additional discussion can be found in Supplemental Material at[END_REF].

In this document, we present additional information to supplement the discussion in the main manuscript, which we will refer to as MS. All references (equations, figures, tables, listings) in the Supplemental Material are labeled with roman numbers.

I. DARWIN AND THE OTHERS

The theoretical framework presented in MS is constructed by abstracting the core algorithm of Darwinian evolution to a general exploration-exploitation (EE) dynamic. To elucidate this analogy, we here provide a glimpse of the main features of the former.

Darwinian evolution, or simply Darwinism, is the widely-accepted theory of biological evolution introduced by the English naturalist Charles Darwin in his seminal work, On the Origin of Species (1859). The problem considered is that of a population (a group of organisms of a species) in the same environment, that reproduce across successive generations. The key ingredient of evolution is inheritance: offspring inherit traits (phenotype) from their parents through genetic information passed down from generation to generation. The different forms in which each gene may exist are called alleles, the pool of alleles of an organism is its genotype. At the population level, two opposing principles cooperate in defining the evolutionary dynamics.

(a) Genetic variation. Stochastic events drive the generation of genetic diversity within a population by introducing variability in the genetic makeup of individuals. This diversity arises primarily through two mechanisms: mutations, which involve random alterations in the genotype, and recombinations, which involve the exchange and shuffling of genetic material between homologous genotypes.

(b) Natural selection. It acts upon this genetic variation by favoring individuals with traits that confer a reproductive advantage in a given environment. As a consequence of selection, populations gradually become better adapted to the environment (adaptation).

General, up-to-date discussions of these concepts can be found in [START_REF] Hamilton | Population genetics[END_REF][START_REF] Manrubia | [END_REF]. Among the existing theoretical approaches to evolutionary dynamics, we mention in particular the one by Neher & Shraiman, i.e., the framework of statistical genetics, proposed in [6] and recently reviewed in [7]. Here, genotypes are modelled as Ising spin-chains g = {σ 1 , . . . , σ L }, where σ i ∈ {-1, 1} and L is the fixed genome length. The evolutionary dynamic is framed in terms of a master equation for P (g), the probability distribution in the genotype space, which changes under the effect of mutations (spin swaps), recombinations (reshuffling of a pair of genomic chains) and selection. The latter is based on a fitness function F (g) ∈ R, which quantifies the aptness of an individual to the environment: at any time, those individuals that are more apt to the environment (higher F ) than the others in the population will have a higher chance to pass their genetic information to the next generation. Note that this implies that selection is state-dependent, since what matters is not the aptness of an individual the environment per se but with respect to that of all others in the population.

In formulating the exploration-exploitation (EE) dynamic, we borrowed the formal structure of the approach above described. The key logical step is to recognize that the evolutionary dynamic is a particular instance of a general EE dynamic in which (i) the configuration space is that of all possible genotypes (ii) exploration is realized by genetic mutations and recombinations and (iii) the exploitation is driven by natural selection. If the interpretations we attach to the concepts of exploration and exploitation are context-dependent, the algorithm coded by evolution to solve its exploration-exploitation problem can be adapted and used elsewhere.

This same line of ideas underlies the family of genetic algorithms, a broad class of computational methods, designed in a way to mimic the evolutionary process to solve complex optimisation problems [START_REF] Goldberg | Genetic algorithms in search, optimization, and machine learning[END_REF]. An evolutionary process is replicated in silico in order to optimise a target objective function. To this end, M individuals (instances of the system) of a population are iteratively mutated, recombined and selected, until the appropriate termination criterion is met. Albeit powerful, this approach is inherently algorithmic, meaning that its interpretation as a direct representation of a real-world process (other than the evolutionary dynamic) is problematic. In particular, outside the context of evolutionary (population) dynamics, it is unclear what the recombinations should correspond to.

In order to retain the potential for interpretation, in eq.(1-2) we implement an EE dynamic where exploration is solely achieved by random mutations, while no recombination-like mechanism exists. In other words, when borrowing from the Darwin's algorithm, we limit ourselves to those mechanisms that are not specific of the evolutionary case [START_REF] Dichio | The exploration-exploitation paradigm: a biophysical approach[END_REF]. Both mutations and fitness-based selection indeed offer a straightforward abstraction -a dictionary can be found in tab.(I) -and an interpretation in different contexts. In sec. V C the case of the C. elegans brain maturation is discussed.

II. THE MATH OF SIMPLE SCENARIOS

Let us start by combining the EE graph dynamic eq.(1-2) in a single formula:

P (G, t + ∆t) = P (G, t) + ∆tµ i<j [P (M ij G, t) -P (G, t)] + e ∆tφF (G) ⟨e ∆tφF ⟩ t -1 P (G, t) . (i)
where all symbols hold as defined in MS. In order to simplify the calculations, it is convenient to define the spin variables σ ij = ±1 which are related to the dyadic variables a ij in MSby the relations:

σ ij = 2a ij -1 , a ij = 1 + σ ij 2 . (ii)
It is important to note that in the context of graphs, the dyads, not the nodes, assume the role analogous to spins in classical statistical mechanics. In a similar fashion, we can define the average magnetization m ∈ [-1, 1] as the spin analogous of the average graph density d ∈ [0, 1] as

m = 1 L i<j ⟨σ ij ⟩ = 2d -1 , d = 1 L i<j ⟨a ij ⟩ = 1 + m 2 . (iii)
The action of an operator M ij that mutates the state of the dyad a ij from 1 → 0 or 0 → 1 (edge toggle) can compactly be expressed as M ij : σ ij → -σ ij , a spin flip. We will consider in this section the regime of weak exploitation in which ∆tφF (G) ≪ 1. A continuous-time description of eq.(i) can then be given:

d dt P (G, t) = µ i<j [P (M ij G, t) -P (G, t)] + φ[F (G) -⟨F ⟩ t ]P (G, t) , (iv)
where we have used e ±x ∼ 1 ± x for x ≪ 1. Eq.(iv) allows to compute the dynamic of the expected value of any graph-observable O(G) : G → R. In fact:

d dt ⟨O⟩ t = d dt G O(G)P (G, t) = G O(G) d dt P (G, t) , (v)
where we have introduced the shorthand

G = σ11=±1 σ12=±1 • • • σ LL =±1 .
Let us now examine some simple cases that are analytically tractable.

A. No exploitation

In the trivial case in which F = const (no exploitation), the dynamic of the system is a random walk in the graph space, whose speed is tuned by the rate µ. We can easily compute the expected value of any spin variable as:

d dt ⟨σ ij ⟩ t = G σ ij µ k<l [P (M kl G, t) -P (G, t)] = µ G σ ij P (M ij G, t) - G σ ij P (G, t) = -2µ⟨σ ij ⟩ t , (vi)
where in the last step we have used

G σ ij P (M ij G, t) = G -σ ij P (G, t).
Eq.(vi) is the differential equation of an exponential decay with characteristic time (2µ) -1 , yielding the solution ⟨σ ij ⟩ t = e -2µt ⟨σ ij ⟩ t0 . In words, after a sufficiently long time, any memory of the initial condition is lost. Each local spin variable converges to a state that is statistically analogous to the outcome of a fair coin toss. Since the dynamics of each variable σ ij are independent, the magnetization follows the same exponential decay m t = e -2µt m t0 . In terms of graph density, this implies that starting from any initial condition, the average state of the system melts down in an Erdős-Rényi random graph where each edge exists with probability p = 1/2.

B. Edge penalty

A simple exploitation scheme is the one in which the only graph sufficient statistic is the number of edges and each existing edge results in a fixed penalty, i.e.,

F (G) = - 1 L i<j a ij , (vii)
where L = N (N -1)/2. The case of a fixed benefit per edge is equivalent to eq.(vii), modulo a minus sign. Using eq.(ii), we can express the exploitation component of eq.(iv) in term of spin variables as

[F (G) -⟨F ⟩ t ]P (G, t) = - 1 2L i<j (σ ij -⟨σ ij ⟩ t ) P (G, t) , (viii)
and again compute the dynamic of the expected value for any spin variable as

d dt ⟨σ ij ⟩ t = -2µ⟨σ ij ⟩ t - φ 2L k<l ⟨σ ij σ kl ⟩ t -⟨σ ij ⟩ t ⟨σ kl ⟩ t (ix)
Differently from the previous case, the dynamics of the average spin variables are now coupled by the sum of covariances C ij;kl = ⟨σ ij σ kl ⟩ -⟨σ ij ⟩⟨σ kl ⟩ in the last term. The calculation simplifies if we restrict to the case in which the spincovariance matrix has approximately a diagonal form, i.e.,

C ij;ij = 1 -⟨σ ij ⟩ 2 ∼ O(1) and C ij;kl ∼ O(ϵ) for ij ̸ = kl (decoupling approximation). Eq.(ix) becomes d dt ⟨σ ij ⟩ t ∼ -2µ⟨σ ij ⟩ t - φ 2L 1 -⟨σ ij ⟩ 2 t , (x)
valid as long as Lϵ ≪ 1, where ϵ is the order of magnitude of the off-diagonal spin covariances. In practice, this means either small graph sizes (small L) or mild functional pressures ρ = φ/µ, for which the dynamic approaches that of the previous case of no exploitation (small ϵ). Under the decoupling approximation the same exact equation as eq.(x) can be written for the magnetization m t and integrated explicitly by partial fractions, yielding the solution:

m t = m 2 1 + m 1 /m 2 -1 1 + m1-m0 m0-m2 e 2µt √ 1+(2L/ρ) -2 , ( xi
)
where m 0 = m t0 and m 1/2 = 2L/ρ ± 1 + (2L/ρ) 2 are the fixed points of the dynamic. Eq.(xi) is an exponential relaxation dynamic to the stable fixed point m 2 , hence m t → m 2 for t → ∞. The magnitude of the asymptotic value m 2 depends exclusively on the ratio between the number of dyads L and the functional pressure ρ. For fixed ρ, the rapidity of the approach to it additionally depends on the specific value of the mutation rate µ, fig.(I). The edge count being the only sufficient statistic of this problem, eq.(xi) completely determines the dynamic. We note that for ρ → 0 we have m 2 → 0 similarly to the previous section, while for ρ → ∞ we have m 2 → -1, corresponding to an empty graph. For any ρ < ∞ the asymptotic value m 2 strikes a balance between the strengths of the exploration and exploitation mechanisms. The value of the asymptotic state d2 solely depends on L, ρ (held fixed) while the rapidity of the approach to it additionally depends on the exploration rate µ.

C. Edge covariate

A straightforward generalization of the case discussed in this section is that to the case in which the F is the sum of an edge covariate, i.e.,

F (G) = 1 L i<j a ij γ ij , (xii)
where Γ = {γ ij } is a N × N real-valued matrix. Each entry of Γ contributes to the F metric only when the corresponding edge exists. For instance, in the case of spatially embedded graphs, the Γ matrix can represent the distance between each possible pair of nodes. A matrix Γ can also quantify homophily effects between nodal attributes, such as the membership to the same community. Following the same steps as above, we get an equivalent expression to eq.(x):

d dt ⟨σ ij ⟩ t ∼ -2µ⟨σ ij ⟩ t + φ 2L γ ij 1 -⟨σ ij ⟩ 2 t , (xiii)
A solution of the same form as eq.(xi) can be written for each spin average ⟨σ ij ⟩. However, each eq.(xiii) depends now on the dyad-specific value γ ij , by consequence the same is true for the fixed points of the dyad dynamics.

D. Distance-like F

We finally discuss an exploitation scheme driven by a distance-like F metric. In this case, an optimal state is explicitly encoded as the global maximum of F , which corresponds to a zero distance between a set of graph statistics x(G) and the target statistics x * , i.e.,

F (G) ∝ -∥x(G) -x * ∥ 2 2 , (xiv)
where ∥•∥ 2 indicates an l 2 -norm. In particular, we consider the simple such case in which the x(G) = ij a ij is one-dimensional and corresponds to the edge count of a graph, the optimal (target) number being E * :

F (G) = - 1 L 2 i<j a ij -E * 2 , ( xv
)
where 0 < E * < L. We can follow the same steps as in sec. II B. The expression analogous to eq.(ix) is

d dt ⟨σ ij ⟩ t = -2µ⟨σ ij ⟩ t - φ L 2 L 2 -E * k<l ⟨σ ij σ kl ⟩ t -⟨σ ij ⟩ t ⟨σ kl ⟩ t + + 1 4 k<l,m<n (kl)̸ =(mn) (⟨σ ij σ kl σ mn ⟩ t -⟨σ ij ⟩ t ⟨σ kl σ mn ⟩ t ) . (xvi)
We generalize the decoupling approximation in sec. II B and define it as the regime in which all moments factorize:

⟨σ (1) σ (2) . . . σ (k) ⟩ ∼ ⟨σ (1) ⟩⟨σ (2) ⟩ . . . ⟨σ (k) ⟩ , (xvii)
where the left-hand side of the equation contains no repeated spin variables. Eq.(xvi) becomes:

d dt ⟨σ ij ⟩ t = -2µ⟨σ ij ⟩ t - φ L 2 1 -⟨σ ij ⟩ 2 t L 2 -E * + 1 2 k<l (kl)̸ =(ij) ⟨σ kl ⟩ t , (xviii)
valid again for Lϵ ≪ 1, where ϵ is the order of magnitude of the spin-covariances. Finally, under the hypothesis of using the same initial condition for the dynamics of all dyads, then we can approximate the last sum ∼ (L -1)⟨σ ij ⟩ t . All dynamical equations are then completely decoupled and an equation of the same form as eq.(xvi) can be written for the magnetization

m t ṁt = -2µm t - φ L 2 (1 -m 2 t ) L -1 2 m t + L 2 -E * . (xix)
The integration of the above expression with boundary condition m t0 = m 0 completely solves the one-dimensional system dynamic and can be done numerically, fig. (II). We find the same qualitative behaviour as in sec. II B, with the solution approaching an asymptotic value that strikes a balance between the influx of mutations -driving the system towards d = 1/2 -and the strength of the exploitation -steering the system towards E * /L. The rapidity of the approach to the asymptotic density value d ∞ depends on µ. where f is specified by eq.(xix). The dynamic has in this case 3 fixed points, labelled as m1, m3 (red, unstable) and m2 (green, stable), found by solving numerically the equation f (m) = 0. Main: dynamic of the average graph density dt (numerical solution). Similarly to the case discussed in sec. II B, the value of the asymptotic state d2 solely depends on L, ρ (held fixed) while the rapidity of the approach to it depends on the exploration rate µ.

0 1 0 1 1 0 1 0 … 2 1 4 1 1 1 0 1 1 0 0 … 0 1 1 0 1 0 0 1 … … … ! ! 𝑛 ! = 𝑀 𝐿 𝒫(𝑡) = (𝐆(t), 𝐧(t)) 0 1 0 1 1 0 1 0 … 5 1 2 0 1 1 0 1 0 1 0 … 0 0 0 0 1 0 0 1 … … … 0 1 0 1 1 0 1 0 … 1 0 6 0 1 1 0 1 0 1 0 … 0 0 0 0 1 0 0 1 … … … ∆𝑡 = 𝜈 "# 𝜇 𝜑 0 0 0 1 0 … 0 0 1 1 1 … 1 0 0 0 1 … … 1 1 0 0 0 … 1 0 1 0 0 … 𝑀 $ (

III. POPULATION-BASED SIMULATIONS

We further exploit the analogy with evolutionary models mentioned in sec. I to design population-based simulations for eq.(i) [29,30]. The current implementation has been coded in Python 3.9.7 and is available on Github [START_REF] Dichio | EE-graph-dyn[END_REF].

In this study, we consider graphs that are unweighted (a ij = 0, 1), undirected (a ij = a ji ), with no self-loops (a ii = 0). Therefore, if N is the (fixed) number of nodes, each graph can be represented by a string of L = N (N -1)/2 binary values 0100 . . . 01 and the configuration space is the L-dimensional hypercube {0, 1} L . Instead of tracking the dynamic of a single instance of the system (individual) in the configuration space, we design a population-based framework where we track simultaneously the evolution of many samples of the probability distribution. More specifically, if G ∈ G is a graph configuration and n ∈ N is the number of times G is observed in the population, we call the couple (G, n) a clone. At any time t, a population P(t) = (G(t), n(t)) is defined as the set of existing clones. The total number of individuals α n α (t) = M is held fixed while the total number of clones M c (t) ≤ M fluctuates. At each simulation step ∆t the population undergoes the processes of mutations (exploration) and selection (exploitation) as follows:

• Exploration. Each possible dyad of each existing graph in the population mutates with probability 1e -∆tµ ∼ ∆tµ for ∆tµ ≪ 1. The mutation (exploration) rate is uniform across dyads. Using spin-variables eq.(ii), in the case of bit-flip mutations we have σ ij → -σ ij . If edges are only allowed to appear,

σ ij → |σ ij | = 1.
• Exploitation. Prior to selection, the F metric of all existing graph configurations G α is computed. The counts n are then updated by performing M independent draws from a multinomial distribution where each G α is chosen with probability

p α = n α e ∆tφF (Gα) / β n β e ∆tφF (G β ) , α ∈ 1, . . . , M c (t) . (xx)
An illustration of the workflow of a simulation step can be found in fig.(III). A pseudo-code for the simulations is the following:

P(0) = (G(0), n(0)) t = 0 while t < T do Exploration: σ ij → -σ ij with probability ∆tµ ∀(i, j), ∀G α update P * = (G * , n * ) compute F α = F (G * α ) ∀ G * α ∈ G * Exploitation: M draws ∼ multinom. distrib. with p α = n * α e ∆tφF (G * α ) / β n * β e ∆tφF (G * β ) ⇒ new counts n * * set P(t) = (G * , n * * ) t = t + ∆t end while
Such a simulation framework has hence two internal free parameters: the total number of samples M and the simulation step ∆t (or, equivalently, its inverse ν = ∆t -1 ). Tab.(??) summarizes all the parameters needed for a single run. A distinction is made between (i) the structural (struct) parameters N, T , which set the dimension of the configuration space and the length of the time window (ii) internal simulation (sim) parameters M, ν, mentioned above and (iii) the dynamical (dyn) parameters µ, ρ, which tune the system dynamic. The simulation time has an obvious linear scaling with the inverse time interval ν, since one run of simulation involves repeating the above steps a number νT of times, where T is the final time. In fig.(IVa) we show that a there is a liner dependence from the number of individuals in a population M , too. Each simulation steps in fact involves mutating M individuals independently and (in the worst case) evaluating the F metric of M different individuals. In fig.

(IVb) we show that, in the limit of large M , a linear dependence is also found on L. This latter result however is strongly dependent on the choice of the F metric, which represents the potential bottleneck of the simulation framework here presented. Additional details can be found in the documentation available at [START_REF] Dichio | EE-graph-dyn[END_REF].

A. Distance-like F , simulated

Let us briefly re-consider the distance-like F metric eq.(xv), discussed in sec. II D. We can set-up the corresponding simulations, use them to showcase the essential features of a single simulation run and compare it to the analytical approximation eq.(xviii).

At each given point in time t, we can use the graph population P(t) = (G(t), n(t)) to approximate the distribution of any any observable O : G → R by where δ is the Dirac-delta function ( dO δ(O) = 1). By consequence, the expected value of O at time t is xxi) is subject to two opposing forces. On the one hand, the influx of random dyadic mutations drives the distribution towards the F of the maximally random state, an Erdős-Rényi random graph with p = 0.5 [31]. On the other hand, the exploitation term increases(decreases) the probability of those F which are above(below) ⟨F ⟩ t , resulting in a net movement of the distribution towards higher F values. This is analogous to a mechanism of adaptation in evolutionary dynamics.

P (O, t) = 1 M Mc(t) α=1 n α (t) δ[O -O(G α (t)] , (xxi)
⟨O⟩ t ∼ 1 M Mc(t) α=1 n α (t) O(G α (t)) . (xxii
In fig.

(Vb) we show the resulting simulated dynamic of d t = ⟨x⟩ t /L and the corresponding confidence interval -whose calculations are based on eq.(xxii). A stationary value is reached where the strengths of the explorationexploitation drivers are balanced. The value of this stationary point and the speed at which it is approached depend on the functional pressure ρ, fig.(Vc). In particular, the higher ρ the closer will be the asymptotic solution to the optimal value, fig. ∞ obtained numerically under decoupling approximation, eq.(xviii). We show evidence that d dec ∞ (decoupling approximation) always lies closer to the target value E * /L than d ∞ , the difference between the two vanishing for increasing exploitation strength φ. We conclude that, for fixed N , the decoupling approximation agrees qualitatively everywhere in the parameter space with simulations, see fig.(Vb), quantitatively for large φ.

IV. EXPONENTIAL RANDOM GRAPH (ERG) MODELS IN A NUTSHELL

Exponential random graph (ERG) models in this work serve the purpose of providing a parsimonious characterization of a given graph G * . They belong to broad class of maximum entropy (maxent) inference methods [32], and are designed for graph data. We here briefly illustrate the essential features of the method, further details can be found in [16,17] and references therein.

(xviii). The d dec

∞ is always closer to the target than d∞, the difference between the two vanishing for increasing φ (where both approach E * /L), as expected.

Let G * ∈ G be a single realization of the system under investigation. Suppose the existence of r sufficient statistics x : G → R r for G * , i.e., that characterize it. According to the maximum entropy principle (MEP), the most unbiased probability density function P (G) consistent with the available information x(G * ) is obtained by maximizing the Shannon information entropy S

S[P ] = - G∈G P (G) log P (G) (xxiii)
while at the same time imposing the normalization G∈G P (G) = 1 and the soft constrains:

G∈G x(G)P (G) = x(G * ) . (xxiv)
This maximization problem is easily solved by introducing a set of Lagrange multipliers θ ∈ R r , and the result is the ERG probability distribution:

P (G|θ) = e θ•x(G) G∈G e θ•x( G) , (xxv)
where the θ are set so to satisfy eq.(xxiv). The parameters θ can be interpreted by quantifying their effect on the likelihood that an edge exists between any pair of nodes i, j in the graph. Let P (a ij = 1|G \ij , θ * ) be the probability of an edge within the dyad (i, j), given the rest of the graph G \ij . We have log

P (a ij = 1|G \ij , θ * ) P (a ij = 0|G \ij , θ * ) = θ * • ∆x(G ij ) , (xxvi)
where ∆x(G ij ) is the vector of change statistics.

Defining G +ij = {a ij = 1, G \ij } and G -ij = {a ij = 0, G \ij }, the α-th element of the change statistics is ∆x α (G ij ) = x α (G +ij ) -x α (G -ij ).
The dyadic interpretation (xxvi) of the parameter θ α is the change in the log probability of a graph, resulting from switching from G -ij to G +ij (i) per unit increase of the corresponding statistic ∆x α (G ij ) = 1, and (ii) holding fixed the cumulative effect of the other statistics β̸ =α θ β ∆x β (G ij ). The interpretation of the ERG parameters is therefore based on a characterisation of the ensamble distribution P (G|θ * ). This because, by construction eq.(xxiv), the properties of the latter reflect those of the original graph. Large positive (negative) values of the parameter θ α indicate an over-(under-) representation in the original graph of the corresponding x α , with respect to the null expectation -i.e., ⟨x α ⟩ of an ERG model with θ α = 0 and unaltered θ \α [START_REF] Dichio | The exploration-exploitation paradigm: a biophysical approach[END_REF]16].

Given a graph G * , the parameters θ of an ERG model can be inferred from data. In theory, this is achieved by maximizing the log-likelihood

θ * = arg max θ log P (G * |θ) . (xxvii)
In practice, the latter computation is hampered by the evaluation of the partition function, i.e., the denominator of eq.(xxv). Nevertheless, the inference problem can be solved numerically by resorting to standard approximation schemes for exponential probability distributions. In particular, Markov chain Monte Carlo maximimum likelihood estimation (MCMC-MLE) iteratively explores the parameter space R r looking for the set of parameters θ * that solves the maximization problem eq.(xxvii). This is done by exploiting eq.(xxvi) to sample graphs from eq.(xxv) with MCMC routines and bypassing the computation of the partition function, an idea dating back to the 90s [33]. A number of packages have been developed in recent years to perform such inference task, the implementation used in this work is the ergm package of the statnet suite [START_REF] Krivitsky | Statnet: Tools for the statistical modeling of network data[END_REF], written in R language. The inferred parameters θ serve the purpose of assigning weights to each statistic and quantifying their contribution towards the generation of the observed graph.

A. ERG statistics

The starting point of an ERG model -as well as of all maxent inferential methods -is the choice of statistics, which is up to the modeller. Choosing a set of statistics entails making an hypothesis about the effects that are relevant for the system. For a given choice of the statistics, the ensemble eq.(xxv) -whose parameters are inferred from datadefines a minimal model in which only the information represented by x is accounted for, and is otherwise maximally unbiased. A huge variety of statistics x : G → R has been proposed and might be included [START_REF] Krivitsky | Statnet: Tools for the statistical modeling of network data[END_REF]. In this work, we have used in particular two kind of statistics. The first, edge covariates, have already been introduced in eq.(xii). We here briefly describe the second, the so-called curved statistics, fig. (VIa).

Let us consider for instance the geometrically weighted degree gwd. It is based on the graph degree distribution and is able to account, e.g., for an over/under representation of hubs in a graph, while at the same time avoiding the so-called degeneracy problem of simplistic ERG models [START_REF] Schweinberger | [END_REF]. Formally, where τ > 0 is a decay parameter and x (k) d (G) is the number of nodes in the graph G with degree k. It is a linear combination of the degree distribution, where the linear coefficients

x gwd (G|τ ) = e τ N -1 k=1 1 -1 -e -τ k x (k) d (G) , (xxviii)
w (k) τ = e τ 1 -1 -e -τ k (xxix)
are based on the geometric series (1e -τ ) k . In order to interpret the effect accounted by eq.(xxviii), we can reason as follows [16,36]. Let us consider a single edge between the nodes i, j. As a consequence of this addition, the degrees of both the incident nodes increase, let us focus on one such increase k → k + 1, using eq.(xxv) one finds:

log P (a ij = 1|G \ij , θ * ) P (a ij = 0|G \ij , θ * ) ∝ θ gwd (1 -e -τ ) k . (xxx)
The role of the ERG parameter θ gwd is the same as described in eq.(xxvi). Let us consider the case θ gwd > 0, which corresponds to the scenario we found in MS. Eq.(xxx) implies that the relative advantage in adding an edge to a node with degree k decreases geometrically with k, fig.(VIb). The rapidity of the decay is controlled by τ , a large τ resulting in a slow decay. It should be noted that, without the geometric decay, an MCMC dynamic based on eq.(xxx) would end up in a fully connected graph (degeneracy), making the estimation impossible. Therefore, for θ gwd > 0, a bias towards high-degree nodes is present, but thanks to the geometric weights it attenuates and eventually disappears for increasing k -where w (k) τ reaches a plateau and (1e -τ ) k ∼ 0. Thus, the system is not steered into an unrealistic fully-connected graph state.

In addition to the gwd, we also introduce the geometrically weighted shared partner distribution:

x gwesp (G|τ ) = e τ N -2 k=1 1 -1 -e -τ k x (k) esp (G) , (xxxi) where x (k)
esp (G) is the number of pairs of connected nodes that share exactly k neighbors. The latter statistic is employed to capture an over/under representation of triads in the network, for θ gwesp ≷ 0. As in the previous case the use of geometric weights helps to prevent this tendency from implying an unrealistic fully connected graph.

V. C.ELEGANS CONNECTOME(S) A. Overview

The C. elegans connectome -i.e., the wiring diagram of neural connections within the worm nervous system -, has garnered significant attention within the field of neuroscience, due to its complete mapping and simplicity [12,37]. It provides valuable insights into the relationship between neural networks and behavior, making it an ideal model system for understanding fundamental principles of brain organization and function.

The nervous system of an adult, hermaphrodite C. elegans was first described in the 1986 by White et al. [10] and consists of only 302 cells (neurons), classified into different types based on their morphology, function, and connectivity. These include sensory neurons that detect various environmental cues, interneurons that relay information within the nervous system, motor neurons that control muscle movement and modulatory neurons that release neuromodulators, playing a crucial role in modulating behavior and neural activity. The neurons of the C. elegans are spatially arranged in a stereotypical fashion, the adult body length being ∼ 1 mm [38]. A high bilateral symmetry is found, most of the neurons occurring in pairs located along the left and right sides of the body. Neurons communicate through different types of cellular junctions. Gap junctions are specialized channels that connect the cytoplasm of adjacent cells, facilitating the rapid exchange of ions, molecules, and electrical signals. Chemical synapses serve as specialised junctions that enable the unidirectional transmission of chemical signals, known as neurotransmitters, from a presynaptic neuron to a postsynaptic neuron or target cell.

This work focused specifically on the maturation of the C. elegans brain [5,14,39]. By this, we refer to the process by which the nervous system develops and matures over the worm's life cycle -the latter includes one embryonic stage (before hatching), four larval stages and one adult stage. At a molecular scale, it involves intricate cellular and molecular events that occur during embryonic and post-embryonic stages. In MS, we developed a model for synaptogenesis, i.e., the process by which synapses are formed and established in the developing nervous system. It involves the growth and differentiation of neuronal processes, such as axons and dendrites, and the precise alignment and interaction of presynaptic and postsynaptic components.

A number of online resources are available that provide comprehensive information and tools for studying the C. elegans nervous system, e.g., Wormatlas.

B. Dataset and pre-processing

In MS, we modelled the data recently published by Witvliet et al. (2021) [5]. In the latter, serial-section electron microscopy was used to reconstruct the brains of eight isogenic C. elegans (wild-type N2, hermaphrodites), imaged at different ages and during different post-embryonic stages, see fig. (VII). The estimation of the age was based on the cell division pattern of the worms at the moment they are selected for imaging. The reconstructed brains, consisting of the nerve ring and ventral ganglion, include 161 of the total 222 neurons at birth and 180 of the total 302 neurons in adulthood, each cell was unambiguously identified by a code and classified as either sensory, modulatory, or interneuron, tab.(III). Chemical synapses and gap junctions were manually annotated. The former were fully mapped while the latter were mapped only partially and, for this reason, they were not considered in our analysis. A directed synaptic connection was defined as a pair of neurons (presynaptic, postsynaptic) connected by at least one chemical synapse. The eight networks of directed synaptic connections were the starting point of our analysis, fig. (VII).

In MS, we used an ERG model, sec. IV, to infer the parameters of an F -landscape, fig.

(1). This required formulating the F metric in terms of graph motifs. In doing this, a natural starting point is looking at the subnetwork distributions of simple motifs, as they are the building block of higher order structures. Previous graph analyses of the adult C. elegans network of synapses [11,12] studied the triad census, i.e., the counts of all directed connection patterns involving three nodes -there exist 16 possible motifs, see fig. (VIII). These investigations found that the patterns involving the same number of empty dyads have similar statistical properties, i.e., they are under-or over-represented with respect to randomized null models. For instance, the triplets of nodes with connections between each pair (regardless of the directions) occur with greater frequency than it would be expected by chance -this is true for the motifs 030T, 120D, 120U, 120C, 210, 300 while 030C represents an exception. Based on this evidence, we simplified the analysis by considering the network of undirected synaptic connections of the C. elegans, fig. (VIII). As a coarsegraining procedure, this implies a loss of information. However, it allows to formulate an F metric eq.( 5) in terms of only a few sufficient statistics, avoiding the combinatorial proliferation of possible patterns for higher-order motifs in directed networks. This should be regarded as a first step towards more complex models of the C. elegans brain maturation, which will account for the directed nature of the information flow.

In tab.(IV), we provide a comprehensive summary of the properties exhibited by the eight undirected graphs analyzed in this study. We observe an increased neuronal connectivity over time, which is further reflected in the growth of the number of simple motifs such as 2-paths and triangles. The rise in the clustering coefficient is consistent with the increasing formation and consolidation of functional circuits or modules [5]. It further suggests the deployment of a biological strategy for increasing the system's redundancy, thereby bolstering the resilience against potential failures. Moreover, we observe a notable increase in efficiency alongside a reduction in the shortest path lengths. This observation aligns well with the documented emergence of highly connected hubs [5], leading to a more tightly woven neural network.

C. The EE dynamics as brain wiring dynamics: an interpretation

In this section, we illustrate how the EE dynamics eq.(1-2) could be interpreted from the point of view of the individual developing system. In so doing, we aim to show how the EE dynamic is consistent with the view of brain wiring dynamics proposed in [START_REF] Hassan | [END_REF], as essentially driven by a set of (genetically encoded) rules that allow neurons to make local decisions about whom to connect with. Note that we take the network representation of the brain as TABLE III. List of the 180 neurons of the adult C. elegans brain (hermaphrodite, N2), as reported in [5]. Interneurons in red, modulatory in yellow, motor in blue, sensory in pink. We have marked with an asterisk * those neurons that were not present at birth. Each neuron in the worm nervous system is uniquely identified by a code, which consists in two or three letters (or, occasionally, numbers), followed by the position in worm's body D/V (dorsal/ventral), R/L (right/left) [10]. The left-right symmetry increases over time and reaches the ∼ 90% in the adult brain. [11,12]. Right:

In an undirected graph, 4 unique connectivity patterns are possible among three nodes. We have employed a color-coding system to illustrate the correspondence between patterns when approximating a directed graph with an undirected one. The biological process we consider in this work -i.e., the formation of neuronal circuits -entails a multitude of complex molecular and cellular events, not fully understood [START_REF] Colón-Ramos | Current topics in developmental biology[END_REF]. In the case of the C. elegans, synapses are mainly found between physically proximal neuronal processes -the latter term generally refers to any projection from the cell . Postsynaptic neuronal processes -namely, dendrites (black lines) and their spine-like protrusions (thin black lines)sprout from a postsynaptic neuron (white circle). Occasionally, they establish synaptic connections with physically proximal presynaptic sites (red dashed lines). The presynaptic processes of the postsynaptic neuron and the postsynaptic processes of the presynaptic neuron have not been depicted. (b) A simple scenario. One presynaptic neuron A and four postsynaptic neurons B, C, D, E are represented. The thick black line represent the axon elongating from A. At time t a synaptic connection exists between the nodes AB (black line connecting B to the axon shaft). Additional connections exist between the neurons BC and DE (summarized by gray dashed line). After a time interval ∆t, postsynaptic neuronal processes extend from the neurons C, D towards the axon, that might develop into new connections (red dotted lines). On the contrary, neuron E does not exhibit the emergence of any such postsynaptic process. (c) Graph representation of the scenario in (b). We hypothesise an F metric that simply counts the number of triangles in the undirected graph representation of the neural system. The two potential connections between AD and AC at time t + ∆t can be represented as two different graph configurations, G1, G2, associated to different F values. G2, by virtue of its higher F , will be observed with higher probability. (d) Example of two time steps of the EE dynamic (∆t = 1 for simplicity). Each square represents a graph. In blue, we indicate the F values.

In black, the unconditioned probabilities computed at each time as exp[F (Gi)]/ j exp[F (Gj)] where the sum runs over all graphs at that time (column). In brown, the probabilities conditioned on the previous time-point. They can be computed either as above, restraining the sum to those graphs that come from the same parent graph at the previous time, as done in eq.(xxxii), or starting from the unconditioned probabilities and using P (Gi, t + 1|Gj, t) = P (Gi, t + 1 ∩ Gj, t)/P (Gj, t), where P (Gj, t) = k P (G k , t + 1 ∩ Gj, t). In bold-red we highlight the most likely developmental pattern.

body: dendrites, postsynaptic or axons, presynaptic. More in detail, most of the neurons have a simple morphology, their processes (one or two per each neuron) run in parallel bundles along the worm's body [10]. Across development, they extend under the influence of molecular guidance cues. Presynaptic sites appear as en passant swellings along the axon shaft. Postsynaptic neuronal processes are either dendrites or spine-like protrusions that grow out from dendrites [5]. Occasionally, dendrites or protrusions establish new synaptic connections, fig.

(IXa).

To demonstrate how the synapse formation process could be encompassed by the EE dynamics, let us take into account the simple scenario illustrated in fig. (IXb). Let us suppose that, at a given time t, synaptic connections exist between neurons AB, DE, and BC. In the time interval ∆t, postsynaptic neuronal processes from both neurons C and D emerge (or extend sufficiently), potentially leading to the formation of new synaptic connections with A. Contrarily, no such process appears on the neuron E, which therefore has no chance of establishing a connection with A. Provided that ∆t is sufficiently small, only one of the two possible synaptic connections AC or AD is likely to be observed. Although the choice is inherently stochastic, a connection offering a higher functional advantage will be associated to a higher probability. Under the hypothesis that the notion of biological function is represented by the number of triangles in the undirected graph representation of the system, fig.(IXc), we anticipate the preferential formation of the AC connection, as it implies the formation of the ABC triangle. In this sense, the "less functional" mutation AD is penalised.

In this interpretation, a mutation event does not correspond to the establishment of a physical connection, but to its precondition (formation/elongation of neuronal process), therefore to a potential connection -note that this marks a difference with respect to the evolutionary dynamics where different genetic mutations correspond to different physical individuals. The exploitation on the other hand acts on configurations that have been allowed by mutations. If G is the graph configuration at time t and G are the potential configurations allowed by mutations at time t + ∆t, then

P (G i , t + ∆t|G, t) = e ∆tF (Gi) / Gj ∈ G e ∆tF (Gj ) , (xxxii)
where we have taken φ = 1 for simplicity. Note that the above probabilities are conditioned on the configuration G observed at time t.

The one above illustrated is an EE dynamic of a single, developing system. From the standpoint of the individual neuron, the process of synaptogenesis consists of a series of stochastic decisions about which other neuron to connect with. These decisions are biased towards those connections that lead to higher functional gains, which in turn are evaluated based on the information available to the neuron at any given time. Note that for this interpretation to be possible, it is essential to define the biological function F in terms of local network structures -e.g., the triangles in the previous example -because neurons only have knowledge of their neighbourhood, not of global system properties.

The EE dynamic contains a comprehensive information about the possible dynamical trajectories, being defined for the probability distribution rather than for a single evolving system. In particular, it allows to compute the (unconditioned) probability of all possible configuration that might have appeared at time t, including those that might have taken very unlikely developmental paths. The difference is illustrated in fig. (IXd). This provides a justification for employing an EE dynamics to encapsulate the inter-individual variations observed in brain wiring.

VI. C. ELEGANS MODEL FIT

The model of the C. elegansbrain growth illustrated in MS has only six parameters, four describe the F -landscape, two the EE dynamics. Their values must be deduced from the data, as detailed in sec. VI A-VI B. Only the birth and adult snapshots of the C. elegans brain are employed for this purpose. The assessment of how well a model's predicted outcomes aligns with the actual observed data is discussed in sec. VI C .

A. F -landscape

The inference of the F -landscape was carried out through exponential random graph (ERG) models, sec. IV, using the two adult snapshots of the C. elegans network of chemical synapses. The ERG inference is defined for a single graph, let us call G * T,i a generic adult snapshot. The model described in eq.(3, 5) can be summarized in the following ERG Hamiltonian H:

-H(G * T,i ) = θ gwd e τ d N -1 k=1 1 -1 -e -τ d k x (k) d (G * T,i ) + θ gwesp e τesp N -2 k=1 1 -1 -e -τesp k x (k) esp (G * T,i ) , (xxxiii)
where we have used the explicit formulation of the statistics gwd eq.(xxviii) and gwesp eq.(xxxi). The above expression has four parameters, namely θ gwd , θ gwesp and τ d , τ esp . Note that from a theoretical point of view, they are fundamentally different. The former are the linear weights of the Hamiltonian that originate from the maxent origin of the ERG probability distribution, sec. IV. Conversely, the latter dictate the formal definition of the graph statistics, therefore they can be ascribed to the problem of model selection. In practice, however, they can be estimated simultaneously by the ergm package [START_REF] Krivitsky | Statnet: Tools for the statistical modeling of network data[END_REF]. We used the following specification:

# ERGM formula fit <-ergm ( formula = G ~gwdegree ( fixed = F ) + gwesp ( fixed = F ) , constraints = ~edges , control = snctrl ( init = c (1 ,1 ,1 ,1) ) )

Listing 1. Specifying an ERG model based on the Hamiltonian eq.(xxxiii) using the ergm package. G is the graph to be used for inference, fixed=F implies that the decay parameters of the curved statistics are to be estimated. The model is constrained to graphs that have the same number of edges as the G. As initial guess of the four parameters to be estimated, we provide (1,1,1,1). For the estimation, we have used Rstudio v2022.12.0.353, R v4.0.4 and ergm v4.3.2. The scripts ara available in the Github folder [START_REF] Dichio | EE-graph-dyn[END_REF].

Note that in the specification of the formula, we constrained the space of graphs explored by the ergm numerical routines to those graphs that have the same number of nodes as the input graph. The reason for this is that the number of edges of the graphs in the EE dynamics is controlled by the exploration rate and therefore it is not a degree of freedom at the disposal of the ERG inference. The results of the estimations are summarized in tab.(V). In the ERG literature, this procedure is referred as mean-ERG [16]. Here, we are essentially limited by the availability of only two individual connectomes.

B. Developmental trajectory

The theoretical picture of the C. elegans brain wiring dynamics discussed in MS is that of a stochastic, statedependent dynamics in the functional landscape F (G) as in eq.( 5), with parameters eq.(xxxiv). The long-term dynamics of this model is quite simple: if we waited long enough, we would end up with (almost) full graphs, independently of our starting point. However, we do not wait long enough. Instead, the edge formation is controlled by the exploration rate µ, and set to reach the appropriate number of edges at adulthood. Therefore, although the maximisation of the F metric learned in the adult stage drives the dynamics, the adult stage itself is not intended to be a maximum of the F landscape.

In MS, we have already discussed the choice of the boundary conditions and the inference of the exploration rate µ * (growth only). This leaves us with only one degree of freedom, the functional pressure ρ, which we discuss in sec. VI B.2. Running a single simulation of the developmental trajectory of the C. elegans requires a few additional specifications.

A single run

If the theory is completely specified by the above parameters, the simulations described in sec. III have two additional degrees of freedom.

The first one, M is the number of samples drawn from the graph probability distribution at each simulation step. Small number of samples will result in high level of noise for the statistics computed as in eq.(xxii). Conversely, a tangible constraint arises from the computational load, which scales linearly with M . In the simulations described in this section, we employed M = 1024 and empirically verified that the outcomes are not significantly altered when simulations are run with M = 2048, M = 4096. The second internal degree of freedom is the inverse time-step ν = ∆t -1 . We fixed it by requiring that in one simulation step a single edge addition was to be observed in each graph within our population. This means choosing ∆t = (Lµ * ) -1 ∼ 4.34 × 10 -2 h, hence ν ∼ 23. Once again, the computational time scales linearly with ν.

Our EE simulations currently do not support a node dynamic (appearance/disappearance). Consequently, to set the initial conditions, we embedded the birth connectome, consisting of 161 nodes, within a larger network that matched the adult connectome's node count of 180.

A single simulation with the configuration here presented takes approximately 1.5 hours. The scripts are available in the Github folder [START_REF] Dichio | EE-graph-dyn[END_REF].

Optimal functional pressure

In MS, we proposed to fix the degree of freedom of the functional pressure ρ by requiring that our simulations optimally reproduce the experimental graph statistics at the adult age. The EE simulated dynamics allows to compute at each time t a whole distribution of graph statistics, while experimental data consists of one or two isolated point per each t. A natural generalization of the Euclidean distance to the distance between a multivariate distribution Q on R r and point y * ∈ R r is the Mahalanobis distance δ mah Q [19], defined as

δ mah Q = (⟨y⟩ Q -y * ) ⊤ Σ -1 Q (⟨y⟩ Q -y * ) (xxxv)
where ⟨y⟩ Q and Σ Q are the mean and covariance matrix of Q. δ mah T used in MS to find the optimal ρ * corresponds to eq.(xxxv) where Q is the two-dimensional distribution of graph statistics eq.( 3) at the adult age T = 45 h (we use the time T to label the distribution), ⟨y⟩ Q = ⟨x⟩ T and Σ Q = Σ T are the ensemble average and covariance matrix, while y * = x(G * T ) are the experimental graph statistics, averaged over the two adult connectomes, i.e.,

y * i = [x i (G * T,1 ) + x i (G * T,2
)]/2. The Mahalanobis distance, as defined in eq.(xxxv), has two intriguing characteristics that render it particularly appropriate for the task at hand. (i) It takes into account the covariance structure of the multivariate distribution. This is crucial, since the statistics eq.( 3) are clearly not independent, fig.(Xa). In particular, let us consider the whitening transformation [START_REF] Kessy | [END_REF] y

→ z = Σ -1 2 Q y ∀y ∈ R r (xxxvi)
where Σ

-1 2 Q
is the inverse principal square root of the covariance matrix Σ Q . It is easy to show that (a) the transformed variables have unit diagonal (white) covariance matrix and (b) eq.(xxxv) corresponds to the Euclidean distance of the transformed variables. In words, the Mahalanobis distance accounts for covariance structure by computing the Euclidean distance of the whitened (standardized) data.

(ii) It is scale invariant, i.e., it is invariant under affine transformations y → Ay + b ∀y ∈ R r , where A is an r × r matrix and b ∈ R r . For instance, it would not change if we scaled the statistics by their corresponding ERG model parameter. Therefore, eq.(xxxv) provides a common ground for comparing models defined by different sets of θ -note that the same is not be true if instead we consider a distance function based on the F metric.

In order determine an optimal functional pressure, we started by scanning uniformly possible values of ρ -100 runs for each. We computed δ mah T for each simulation, averages and standard deviations are depicted in fig. (2a). As expected, for mild functional pressures, the simulations do not attain the desired value at adulthood, hence the Mahalanobis distance increases for low values of ρ. Furthermore, it also increases for too high values of ρ, implying an overestimation of the functional constraints deployed by the biological system.

We fitted the data with a quadratic curve δmah T (ρ) = aρ 2 + bρ + c and took the position of the minimum -b/2a as optimal functional pressure ρ * . This last methodological step brought the inference process to completion. Finally, in order to showcase the results of a single simulation of the C. elegans brain maturation in fig. (2c-2d), we ran 500 simulations with the inferred parameters µ * , ρ * and again selected as best simulation the one that minimized δ mah T . Similarly to eq.(xxxv), one can compute the distance δ mah t of the simulated distribution of statistics from the corresponding experimental values at all observed time points t * = 5 h, 8 h, 16 h, 23 h, 27 h and 45 h. Note that by assumption δ mah t=0 = 0 and that, except for the adult age, only one experimental graph is available, therefore y * = x(G * t ). The total error over the whole (observed) time series can be simply defined by

δ mah t * = t∈t * , t>0 δ mah t . (xxxvii)
Using the latter, we repeated the same steps described in the previous paragraph. The minimum of the quadratic fit of δmah t * (ρ) was found to be located at ρ * * = 7.001 × 10 2 . In fig.

(2a), we plotted δmah T (ρ * * ), i.e., the average value of the Mahalanobis distance defined using the adult age exclusively that we would have obtained if we had chosen ρ * * . This value fell within one error bar from δmah 3). At any given time t our simulation allow to extract a two-dimensional distribution in the space gwd-gwesp. We show those corresponding to four of the seven experimental time-stamps. Note that the statistics are (anti-)correlated, in computing a distance from the distribution, the covariance structure must be taken into account. Here we used the same simulation as in fig. (2b-c (blue line and shaded area, respectively). The data are fitted with a parabolic curve δmah t * (R 2 = .99) and its minimum ρ * * is highlighted (blue diamond). On the same curve, we show the position of δmah t * (ρ * ), the value of the all-time Mahalanobis distance that corresponds to the functional pressure ρ * we selected in MS, based on the adult stage exclusively (orange diamond). The two overlap within the error bars. This plot is complementary to fig.(2a) in MS.

C. Feature generalization

By feature generalization, we refer to the ability of a model to reproduce other features of the data than those included in the model formulation. We here briefly provide definitions of the metrics used in this work to this purpose, see also [START_REF] Latora | Complex networks: principles, methods and applications[END_REF] for a general reference. The package NetworkX v2.6.3 with Python 3.9.7 has been used, the scripts are available in the Github folder [START_REF] Dichio | EE-graph-dyn[END_REF].

Let us consider an unweighted, undirected graph G with N nodes, and adjacency matrix A = {a ij }. Let k i = j a ij be the degree of the the i-th node.

• Clustering. A first family of graph metrics aims at computing the extent to which nodes in a graph tend to cluster together. The most straightforward manifestation of clustering is a higher-than-random probability that two nodes that are connected to a common node are also connected to each other. Such a behaviour can be quantified by the following two metrics:

• Transitivity, fig.(XIa). It evaluates the ratio between the number of existing triangles -triples (i, j, k) with a ij = a jk = a ik = 1 -and the number of connected triples, i.e., the triples (i, j, k) with a ij = a jk = 1. The symmetry factor accounts for the fact that each triangles has 3 connected triples. Formally,

T = 3 × # triangles(G) # connected triples(G) = tr(A 3 ) i̸ =j (A 2 ) ij . (xxxviii)
• Average clustering coefficient, fig. (2d). It is defined as: . Given a statistic y, at each experimental time point t ∈ t * we compute its experimental value y(G * t ) (circles), the average ⟨y(G)⟩t (diamonds) based on our population-based simulations eq.(xxii), and the same average with simulations based on a null growth model with ρ = 0 (squares). In both cases our simulations closely track the experimental values. (b) Cumulative degree distributions at adulthood, T = 45 h. Those corresponding to the two adult C. elegans brain graphs are represented as green solid lines. Dashed orange line for the simulated cumulative degree distribution, bin-wise average over the distribution of graphs we obtain from simulations. The Kolmogorov-Smirnov statistics between the latter and the two empirical distributions are DKS = 0.08, 0.10, respectively. Gray dashed line for the corresponding values from a null model with ρ = 0, DKS = 0.24, 0.18. Inset: zoom in the high-degree tail, loglog plot. The black line is a linear fit of the simulated high-degree tail of the distribution.

C = 1 N N i=1 C i , C i = # connected pairs of neighbors of i # pairs of neighbors of i = j,l a ij a jl a il k i (k i -1) . (xxxix
where k i = j a ij is the degree of the node i. In words, C is the average of the local clustering coefficient C i . The latter is computed by considering the subnetwork induced by the node i and its first neighbors and quantifies the relative number of neighbors of i that are also themselves neighbors.

Albeit both T and C take values in the interval [0, 1] and reach 1 in the case of perfect transitivity, they are not equivalent. The clustering coefficient is more influenced by low-degree nodes (since it averages over all nodes), while transitivity is more influenced by high-degree nodes (since it considers all possible triangles in the network). Therefore, transitivity might be more representative of the overall network structure, while the clustering coefficient could provide more insight into local structures or subnetworks within the network.

• Efficiency. A second family of metrics aims at quantifying the effectiveness of information or resource exchange over the network. As a general idea, the closer are two nodes in the graph, the more efficiently information will be exchanged among them.

• Global efficiency, fig. (2d). A convenient way to quantify a graph's efficiency is to compute the harmonic mean of geodesic lengths. More in detail, let D be the matrix whose elements d ij represent the shortest path (geodesic) from the node i to j -by definition; d ij = ∞ for disconnected nodes. The global efficiency E g is defined as

E g = 1 N (N -1) i̸ =j 1 d ij . (xl)
A network with high global efficiency is typically characterized by short paths between any given pair of nodes, meaning information or resources can be disseminated rapidly across the network. This property is often observed in random or in 'small-world' networks.

• Local efficiency, fig.(XIa) .

A complementary definition of efficiency can be given by averaging over all nodes a local notion of efficiency. The local efficiency E l can be defined as

E l = 1 N i E (i) g , (xli)
where

E (i)
g is the global efficiency of the subgraph induced by the node i and its neighbors. By looking at the efficiency of each node's immediate subnetwork, local efficiency provides insight into the network's robustness or resilience to failures or attacks. If a network has high local efficiency, the removal of a node would not significantly disrupt communication between its neighboring nodes.

• Node degrees. The degree distribution of a network is a critical aspect to consider when studying network structures, as it encapsulates fundamental information about the network's structure, robustness to failures, and information spreading dynamics, fig.(XIb).

Given a network, one may look directly at its degree distribution P deg (k) = n k /N where n k is the number of nodes with degree k, or its cumulative distribution

P (c) deg (k) = P deg (j ≥ k) = 1 N j≥k n j , n j = # nodes with degree j . (xlii)
Cumulative degree distributions are frequently employed in network analysis, especially for networks with heavytailed degree distributions typical of many real-world systems. In fact, they offer an enhanced visualization of the degree structure, reducing noise in the distribution's tail and facilitating the identification of long-tail or powerlaw behaviors indicative of scale-free networks. The distance between two cumulative distributions P (c) , P (c) can be quantified e.g. by computing the Kolmogorov-Smirnov statistic:

D KS = sup k |P (c) (k) -P (c) (k)| . (xliii)

VII. A ROADMAP

In this section, we point more in detail at relevant lines of development of our work.

A. More detailed models of a worm brain

The EE model of the C. elegans brain wiring dynamics designed in MS lends itself to a number of extensions for achieving more detailed depictions of the underlying process. To start with, when formulating the biological function F , simplicity can be traded for realism in several ways.

In the first place, the system's representation. For instance, a natural extension of our approach would be to take into account the directed nature of chemical synapses, by considering directed dyads (i → j), hence directed graphs. From a technical point of view, this would entail storing 2L = N (N -1) bits per graph -i.e., twice the memory L needed for an undirected graph. The statistics eq.( 3) in MS admit generalisations to the directed case, which however require non-trivial modelling choices, see [START_REF] Krivitsky | Statnet: Tools for the statistical modeling of network data[END_REF]. Directed graphs would allow modelling the existence of a feedforward bias for synaptogenesis that increases with time, i.e. new synaptic connections mainly appear in the direction: sensory neurons → (modulatory neurons →) interneurons → motor neurons [5]. As usual, this should be specified by including a suitable term in the formulation of F .

A general and useful ERG-like term to model effects is the edge covariate:

x(G) = 1 2L i̸ =j a ij γ ij , (xliv)
here in the case of directed graphs, cf eq.(xii). The latter is a sum over all existing (directed) edges of the edgeattribute γ ij . The latter can be used to specify a number of different effects. For this to be possible, the information about γ ij of each possible dyad must be available. Three examples are the following.

• Feedworward bias. We define γ ij = 1 if the synaptic connections goes from sensory to motor neurons, γ ij = -1 if it goes in the opposite direction, γ ij = 0 if it connects neurons of the same time. Then, eq.(xliv) quantifies the feedforward bias.

• Homophilies. This refers to the fact that neurons in the worm's brain are more likely to form a connection if they share common attribute, see [14] for an overview. For example, there is a bilaterally symmetric pairing homophily, i.e., neurons are more likely to be connected if they belong to a bilaterally symmetric pair, tab. III.

To account for such effect, we could use eq.(xliv) with γ ij = δ(α i , α j ), where δ is the Kroenecker delta, and α i = L, R if the i-th neuron belongs to the left or right part of the worm body, respectively.

• Spatial embedding. The worm brain network is embedded in the physical space [18]. Long connections are anticipated to incur higher penalties due to increased wiring costs. The simplest way to take this into account would be to use eq.(xliv) with γ ij corresponding to the euclidean distance between the cell bodies of each pair of neurons.

Furthermore, the biological function F can be endowed with additional terms that can account for the presence of stable neuronal circuits in the wiring [5]. For example, suppose that the (directed) triangle i → j → k → i is of particular biological importance to the worm's nervous system and is therefore expected to be disproportionately conserved across individuals. A term of the form x(G) = 1a ij a jk a ki could be added to F . In the case of a negative associated parameter θ, this term would reduce the likelihood of those graphs that do not have the above circuit.

According to the spirit of the approach here illustrated, the idea would be to include these effects (some, all, more) in the formulation of the F metric and let the ERG estimation "judge" which are the most relevant, as quantified by the associated parameter θ. Large magnitudes would indicated relevancy, while θ ∼ 0 would indicate marginal or no contributions to the observed network. More detailed models will allow for reproducing more detailed features of the data. A systematic exploration of the space of possible models along the lines outlined above would holds considerable interest and we envision this as a promising avenue for future investigations.

Zooming out to the EE framework, several possible refinements could improve the model realism. A straightforward extension would implement the appearance of new neurons at different stages during the developmental dynamics, as mentioned in sec. VI B 1. A less straightforward generalisation would be to abandon the hypothesis of time homogeneity of the process. For instance, this would mean to use time dependent EE parameters µ(t), ρ(t), or even a time dependent functional landscape -or, functional seascape -F (G, t). Here we are currently limited by the availability of only a few snapshots of the worm brain, which limits the complexity of the hypotheses that can be meaningfully made and tested. However, given the rapid advances in experimental imaging techniques, this family of generalisations of the current approach may be within reach in the near future.

Finally, if and when data will be available, our model could be readily extended to different kinds of brain connectivity, as for example the network of gap junctions within the C. elegans brain, or its functional connectivity [24]. Also, an exciting line of investigations would entail developing a similar model for different natural brains, whose reconstruction are currently at different levels of completion: larval zebrafish, Drosophila, mouse (...) [START_REF] Hildebrand | [END_REF]42,43]

B. Beyond the brain wiring problem

The EE dynamics arise from the interplay between a stochastic exploration of the configuration space and a statedependent exploitation of the more functional states, regardless of the data representation. It encodes the interplay between randomness and functional robustness of the dynamics. The brain wiring dynamics, we have argued, fits this theoretical picture. However, EE dynamics can be used to describe also other types of dynamics, that can be placed within the same conceptual framework. In the following, we anticipate some concrete possible examples (biological and non-biological) that could be studied under the EE lens.

• Human brain diseases. Within connectomics applications, different questions can be asked, at different scales of neuronal organisation. EE dynamics could be used to understand how the human brain -as reconstructed from electrophysiological or neuroimaging experiments -longitudinally reorganises after a traumatic event. For example, during the recovery from stroke, the brain tends to spontaneously explore local rewiring strategies to reestablish functional connections both within the affected hemisphere and with the unaffected hemisphere [44,45]. Although at a different scale of neural organisation -here nodes are brain regions, edges measure structural or functional connectivity -a similar approach to that developed in MS could be used to model individual recovery/reorganisation dynamics, with the same rationale.

• Social contacts. Social bonds between individuals are essential for many living organisms. Modeling the resulting network formation is therefore crucial to understand, predict, and eventually influence the overall system behavior. For example, the formation of a social network on digital platforms such as Twitter occurs when users are exposed to various potential connections, often through suggested posts in their feeds or recommended profiles. Some lead to the establishment of new social bonds, typically in the form of 'follow' relationships.

The probability of this process is influenced by several factors, including the degree of homophily among users and the presence of mutual acquaintances [START_REF] Mislove | Proceedings of the 7th ACM SIGCOMM conference on Internet measurement[END_REF]. When representing the social network as collections of nodes (users) and directed edges (indicating 'follow' relationships), these dynamics can be conceptualized as resulting from two key drivers: (1) a non-uniform, stochastic exploration of the graph space and (2) a state-dependent exploitation process driven by a ERG-like F accounting for homophilies (edge covariates) and graph-topological effects -e.g., the presence of hub nodes (popular individuals) and clustering (well-connected communities).

• Spatial navigation. Beyond the network representation, the EE paradigm could be used also to study the problem of spatial navigation in complex environments. For example, the slime mold Physarum polycephalum uses a form of reactive navigation to explore its environment. In the presence of attractants like food, its cytoplasm flows stochastically toward the food source, while repellents like light or salts push it away. Moreover, as the slime mold forages, it leaves behind a trail of nonliving extracellular slime, which it later avoids. This simple, noneuronal organism is able to solve complex navigation tasks, such as the U-shaped trap problem [START_REF] Reid | [END_REF]. In terms of an EE dynamics, this can be modeled as (i) a stochastic exploration of the physical space, and (ii) a state-dependent exploitation driven by a F metric which enforces hard constraints (physical barrier and previous path), soft constraints (repellents) and favors those displacements that minimise the distance from the food source.

FIG. 1 .

 1 FIG. 1. EE dynamics for the C. elegans brain maturation. We consider the eight snapshots of the worm network of chemical synapses, at different developmental ages t = 5 h, 8 h, 16 h, 23 h, 27 h and 45 h (two adults), from [5] -some are omitted for visual clarity. The birth configuration is fixed as the starting point of the dynamics. The two adult snapshots are used to infer (i) the topography of the functional landscape F (G) (bottom left), encoded in the set θ * -ERG inference starting from the empirical statistics x * ≡ x(G * T ); (ii) the EE parameters, i.e., the exploration rate µ * and the functional pressure ρ * , see SM-VI.

FIG. 2 .

 2 FIG. 2. Tracking down the C. elegans brain maturation. (a) We run 100 simulations ∀ρ ∈ {200 + 20 i, 0 ≤ i ≤ 60}. For each ρ, we compute the mean and standard deviation of δ mah T (red line and shaded area, respectively). We fit the data with a quadratic curve (dash-dotted line) and take the abscissa ρ * of its minimum (red diamond) as estimation of the functional pressure. On the same curve, we show the value (blue diamond) corresponding to the abscissa ρ * * we get by minimizing the sum of the Mahalanobis distances over all experimental time-points, SM-VI.B. The two overlap within the error bars. (b) The exploration rate is calculated by using the first and last time point (average). The shaded area corresponds to the estimation by a linear fit over the whole time series. (c) One simulation run with µ * , ρ * . Main: the trajectory in the space of statistics (gwd, gwesp). Experimental data (circles) are closely tracked by our simulations (dashed line, diamonds highlighting the observed time points t * ). The trajectory of a null model with ρ = 0 is also shown (dotted line, squares) Inset: The simulated and experimental time course of F in time. (d) Feature generalisation. The temporal trajectory of the average clustering coefficient and global efficiency, markers as described in (c). See also SM-VI.C.

FIG

  FIG.I. Analytical solution of the exploration-exploitation dynamic based on eq.(ix), decoupling approximation. Here N = 10 (L = 45), μ = 1/L, ρ = 2 × 10 2 , d0 = 0.9. Inset: plot of ṁ = f (m), where f is specified by eq.(x). The dynamic has two fixed points, corresponding to the solutions of f (m) = 0, labelled as m1 (red, unstable) and m2 (green, stable). Shaded area corresponds to the inaccessible regions for m ∈ [-1, 1]. The initial condition m0 (white diamond) and the one-dimensional vector field (green arrows) are indicated. Main: dynamic of the average graph density dt as obtained from eq.(xi) using eq.(iii). The value of the asymptotic state d2 solely depends on L, ρ (held fixed) while the rapidity of the approach to it additionally depends on the exploration rate µ.

  FIG.II. Analytical solution of the explorationexploitation dynamic based on eq.(xix), decoupling approximation. Here N = 10 (L = 45), μ = 1/L, ρ = 5 × 10 2 , d0 = 0.9, E * = 3. Inset: plot of ṁ = f (m), where f is specified by eq.(xix). The dynamic has in this case 3 fixed points, labelled as m1, m3 (red, unstable) and m2 (green, stable), found by solving numerically the equation f (m) = 0. Main: dynamic of the average graph density dt (numerical solution). Similarly to the case discussed in sec. II B, the value of the asymptotic state d2 solely depends on L, ρ (held fixed) while the rapidity of the approach to it depends on the exploration rate µ.

  FIG. III. A simulation step of population-based simulations for eq.(i). A population P(t) at time t consists of matrix G(t)with dimensions Mc(t) × L and a vector n(t) of length Mc(t) (top left). Each row of G(t) corresponds to a unique graph configuration, and in particular to the vectorized lower triangular matrix of the corresponding adjacency matrix (top right). First, the population undergoes mutations, with rate µ per each dyad (green, bottom right), here bit-flips. Consequently, both the matrix G and the counts n change. As new configurations are found, the number of clones Mc increases. Second, selection is implemented, tuned by the parameter φ (violet, bottom left). Selection only affects the vector of counts n, those graphs that end up with n = 0 are removed from G before the next time-step. The total number of individuals M and time interval ∆t are free internal parameters of the simulations.

  FIG. IV. Simulation time as a function of the number of individuals M and L = N (N -1)/2, where N is the number of nodes. Each dot is the average value obtained from 10 simulation runs, error bar are indicated (when visible). Here T = 100, ν = 1, µ = 1.0 × 10 -4 , ρ = 5, an F metric as in sec. II B is used. (a) The simulation time is linearly dependent on the number of individuals M . Dotted lines are the result of a linear fit with the curve y = a * M , where a is a parameter. A log-log scale has been used. (b) The simulation time for this problem has a mild exponential dependence on L. Dotted lines are the result of a fit with the curve y = a * L b , where a, b are parameters. In the limit of large M , we find β ∼ 1, indicating an approximately linear dependence on L. The dependence on L (N ) however, strongly depends on the complexity of the F metric.

  Fig.(Va) illustrates the EE dynamic in the space of the F values. At any given time t the F -distribution eq.(xxi)is subject to two opposing forces. On the one hand, the influx of random dyadic mutations drives the distribution towards the F of the maximally random state, an Erdős-Rényi random graph with p = 0.5[31]. On the other hand, the exploitation term increases(decreases) the probability of those F which are above(below) ⟨F ⟩ t , resulting in a net movement of the distribution towards higher F values. This is analogous to a mechanism of adaptation in evolutionary dynamics.In fig.(Vb) we show the resulting simulated dynamic of d t = ⟨x⟩ t /L and the corresponding confidence interval -whose calculations are based on eq.(xxii). A stationary value is reached where the strengths of the explorationexploitation drivers are balanced. The value of this stationary point and the speed at which it is approached depend on the functional pressure ρ, fig.(Vc). In particular, the higher ρ the closer will be the asymptotic solution to the optimal value, fig.(Vd). Finally, in fig.(IIe) we compare the asymptotic values d ∞ resulting from simulations with the corresponding values d dec ∞ obtained numerically under decoupling approximation, eq.(xviii). We show evidence that d dec ∞ (decoupling approximation) always lies closer to the target value E * /L than d ∞ , the difference between the two vanishing for increasing exploitation strength φ. We conclude that, for fixed N , the decoupling approximation agrees qualitatively everywhere in the parameter space with simulations, see fig.(Vb), quantitatively for large φ.

  Fig.(Va) illustrates the EE dynamic in the space of the F values. At any given time t the F -distribution eq.(xxi)is subject to two opposing forces. On the one hand, the influx of random dyadic mutations drives the distribution towards the F of the maximally random state, an Erdős-Rényi random graph with p = 0.5[31]. On the other hand, the exploitation term increases(decreases) the probability of those F which are above(below) ⟨F ⟩ t , resulting in a net movement of the distribution towards higher F values. This is analogous to a mechanism of adaptation in evolutionary dynamics.In fig.(Vb) we show the resulting simulated dynamic of d t = ⟨x⟩ t /L and the corresponding confidence interval -whose calculations are based on eq.(xxii). A stationary value is reached where the strengths of the explorationexploitation drivers are balanced. The value of this stationary point and the speed at which it is approached depend on the functional pressure ρ, fig.(Vc). In particular, the higher ρ the closer will be the asymptotic solution to the optimal value, fig.(Vd). Finally, in fig.(IIe) we compare the asymptotic values d ∞ resulting from simulations with the corresponding values d dec ∞ obtained numerically under decoupling approximation, eq.(xviii). We show evidence that d dec ∞ (decoupling approximation) always lies closer to the target value E * /L than d ∞ , the difference between the two vanishing for increasing exploitation strength φ. We conclude that, for fixed N , the decoupling approximation agrees qualitatively everywhere in the parameter space with simulations, see fig.(Vb), quantitatively for large φ.

  FIG. V. Exploration-exploitation (EE) dynamic as defined by eq.(i, xv), N = 10, μ = 1/L, ρ = 5 × 10 3 ( φ = ρμ), T = 50, ∆t = 1/ν = 10 -2 , E * = 3, ∆t = 1, M = 4096. (a) Distribution of F values at different t, the corresponding ⟨F ⟩t are indicated. At each t, exploitation increases the probability of those graphs with higher F than the ensemble average (violet arrows). As a result, the distribution moves towards higher F values, until the max is reached (red dotted line). (b) The average graph density ⟨x⟩/L evolves towards the target density E * /L, until a stationary value ⟨x⟩∞/L is reached -shaded area indicates the 95% confidence interval, black dashed line indicates the analytical solution from eq.(xix). (c) The value ⟨x⟩∞/L and the rapidity of the approach to it depends on the functional pressure ρ, it approaches E * /L for increasing exploitation rate φ and decreasing exploration rate µ. (d) The distance from the target increases for increasing exploration rate µ and decreasing selection rate φ. (e) Difference between d∞ as in (d) and the corresponding value d dec ∞ decoupling approximation eq.(xviii). The d dec∞ is always closer to the target than d∞, the difference between the two vanishing for increasing φ (where both approach E * /L), as expected.

  FIG. VI. Geometrically weighted statistics. (a) An illustration of the gwd and gwesp statistics, as defined in eq.(xxviii-xxxi). Both have the structure of a linear combination of a distribution (degrees, edgewise shared partners), the weights w (k) τ are defined in eq.(xxix). (b) Weights w (k) τ(solid lines) and geometric decaying factors (dashed lines), as they appear in eq.(xxx). Consider as an example the case of xgwesp, in the case θgwesp > 0. Adding a new partner to a pair of connected nodes k → k + 1 always leads to an increase in the probability associated to the graph, eq.(xxx), since eq.(xxix) is a monotonically increasing function of k. However, this increase diminishes and eventually vanishes for increasing k.

  FIG. VII. C. elegans brain networks at different developmental ages. The embryonic stage (E, dotted line) terminates with the hatching (birth, 0 h, white circle). The post-embryonic stages include four larval ages (L1 -L4, dashed line) and ends with the onset of adulthood (red circle). The adult stage (A, solid line) lasts about 2 -3 days. The dataset Witvliet et al. (2021) [5] consists of eight snapshots (microscope icons), including one at birth 0 h, three L1 ∼ 5 h, 8 h and 16 h, one L2 ∼ 23 h, one L3 ∼ 27 h and two adults, both ∼ 45 h. In each network, nodes are individual neuron cells, colored by cell-type: interneurons in red, modulatory in yellow, motor in blue, sensory in pink -see tab.(III) for a complete list. A directed edge (connection) is placed between two neurons if at least one synapse exists between them.

15

 15

F

  ! , 𝑡 + ∆𝑡|𝐺, 𝑡 < 𝑃 𝐺 " , 𝑡 + ∆𝑡|𝐺,

T

  (ρ * ). In fig.(Xb), we corroborate this result by showing the plot complementary to fig.(2a), for eq.(xxxvii).

  FIG. X. (a)Population-based simulations, in the space of model statistics eq.(3). At any given time t our simulation allow to extract a two-dimensional distribution in the space gwd-gwesp. We show those corresponding to four of the seven experimental time-stamps. Note that the statistics are (anti-)correlated, in computing a distance from the distribution, the covariance structure must be taken into account. Here we used the same simulation as in fig.(2b-c) -these are the underlying distributions of the average values plotted in fig.(2b). (b) We run 100 simulations ∀ρ ∈ {200 + 20 i, 0 ≤ i ≤ 60}. For each ρ, we compute the mean and standard deviation of δ mah

  FIG. X. (a)Population-based simulations, in the space of model statistics eq.(3). At any given time t our simulation allow to extract a two-dimensional distribution in the space gwd-gwesp. We show those corresponding to four of the seven experimental time-stamps. Note that the statistics are (anti-)correlated, in computing a distance from the distribution, the covariance structure must be taken into account. Here we used the same simulation as in fig.(2b-c) -these are the underlying distributions of the average values plotted in fig.(2b). (b) We run 100 simulations ∀ρ ∈ {200 + 20 i, 0 ≤ i ≤ 60}. For each ρ, we compute the mean and standard deviation of δ mah t *

)

  FIG. XI.Feature generalization assessment of the EE dynamics. (a) Transitivity (brown) and local efficiency (yellow). Given a statistic y, at each experimental time point t ∈ t * we compute its experimental value y(G * t ) (circles), the average ⟨y(G)⟩t (diamonds) based on our population-based simulations eq.(xxii), and the same average with simulations based on a null growth model with ρ = 0 (squares). In both cases our simulations closely track the experimental values. (b) Cumulative degree distributions at adulthood, T = 45 h. Those corresponding to the two adult C. elegans brain graphs are represented as green solid lines. Dashed orange line for the simulated cumulative degree distribution, bin-wise average over the distribution of graphs we obtain from simulations. The Kolmogorov-Smirnov statistics between the latter and the two empirical distributions are DKS = 0.08, 0.10, respectively. Gray dashed line for the corresponding values from a null model with ρ = 0, DKS = 0.24, 0.18. Inset: zoom in the high-degree tail, loglog plot. The black line is a linear fit of the simulated high-degree tail of the distribution.

TABLE I .

 I Translation dictionary between exploration-exploitation EE graph dynamic (left) and evolutionary dynamic (right). We group by colour those terms that refer to the structure of the configuration and state spaces (top), to the the structure of the dynamics (middle) and to the dynamic parameters (bottom). No equivalent mechanism to genetic recombinations exists in the EE graph dynamic. *See also sec. III.

	EE graph dyn.	Notation	Evolutionary dyn.
	graph space number of dyads	G L	genotype space genome length
	graph (unweighted, undirected)	G/g	genotype (biallelic)
	graph statistics	x	phenotype (traits)
	dyad (bit-like)	aij/σij	locus (spin-like)
	time window # samples *	T M	# generation population size
	biological function (F metric)	F	fitness function
	exploitation r.	φ	natural selection r.
	exploration r.	µ	mutation r.
	×	r	recombination r.

TABLE II .

 II Parameters of simulations for EE graph dynamics. Our computational framework has six degrees pf freedom, which we group by color: structural parameters (top), internal degrees of freedom (middle) and parameters of the dynamics (bottom).

	Parameter Description
	N	number of nodes
	T	time window span
	M	population size
	ν	inverse time step ∆t -1
	µ	exploration rate
	ρ	functional pressure φ/µ

TABLE IV .

 IV Properties of the C. elegans networks of undirected synaptic connections. Each row corresponds to a graph, the first (birth) and the last two (adulthood) are highlighted. We compute the number of nodes, edges, two-star -or connected triples -, triangles, the average shortest path (av.sh.path) -or average geodesic distance -, the global efficiency (glob.eff.) and the average clustering coefficient (clust.coeff.). See also sec. VI C for the definitions.t[h] nodes edges two-star triangles av.sh.path glob.eff. clust.coeff.

	0	161	617	5976	346	2.993	0.380	0.208
	5	162	782	9273	601	2.712	0.416	0.232
	8	162	788	9299	614	2.712	0.416	0.245
	16 168	907	11838	830	2.617	0.428	0.246
	23 173	1166	18449	1406	2.430	0.459	0.262
	27 174	1175	18866	1433	2.429	0.458	0.274
	45 180	1633	34124	2889	2.217	0.498	0.286
	45 180	1669	35677	3003	2.206	0.501	0.292

correct and focus on its growth (edge addition) during development. Therefore, a model of the brain wiring process at the molecular level (axon guidance, growth factors, etc.) is neither possible nor desired here.

  FIG. IX. Individual-based interpretation of the EE dynamic. (a) Establishment of a synapse, schematics. Presynaptic sites, presented as en passant swellings (black circles), appear on the axon shaft (thick black line) of a presynaptic neuron (gray circle)

				2
				.11
				.67
			0	0
			.27	.01
			.27	.09
		0		1
		1		.04
				.24
	𝐺 ! , 𝑡 + ∆𝑡	𝐺 " , 𝑡 + ∆𝑡		4
				.80
			1	.95
			.73
			.73	1
				.04
				.05
	𝑡 + ∆𝑡	𝑡 𝐹 ! = 0	𝑡 + ∆𝑡 𝐹 !"∆! = .5	𝑡 + 2∆𝑡 𝐹 !"$∆! = 1.6

TABLE V .

 V Estimated parameters of the ERG model eq.(xxxiii) from the two adult C. elegans (undirected) networks of chemical synapses. The maxent parameters θ * are both significant and positive for all networks. The parameters λ * controlling for the geometric decays of the model statistics are significant -and positive by construction. The third row represents the mean-ERG, final parameter estimations for the EE dynamic are boldfaced.An output procedure is needed to combine the estimations from the group of adult individuals in our dataset. Considered that each individual yields independent estimation of the ERG parameters, a simple choice is to take the average estimations, i.e.,

	θ * gwd	λ * gwd	θ * gwesp	λ * gwesp
	G * T,1 0.45 ± 0.20 1.91 ± 0.46 0.626 ± 0.056 1.432 ± 0.067
	G * T,2 0.43 ± 0.20 1.97 ± 0.48 0.529 ± 0.048 1.542 ± 0.075
	ERG 0.44 ± 0.14 1.94 ± 0.33 0.578 ± 0.037 1.487 ± 0.050
	θ * gwd = 0.44, τ * d = 1.94, θ * gwesp = 0.578, τ * esp = 1.487	(xxxiv)

We thank E Aurell, E Mauri, D Battaglia and M Josserand for the many useful discussions. FDVF ac-knowledges support from the European Research Council (ERC), Grant Agreement No. 864729.

Supplemental Material

The exploration-exploitation paradigm for networked biological systems