Symmetries and synchronization from whole-neural activity in {\it C. elegans} connectome: Integration of functional and structural networks
Bryant Avila Pedro Augusto David Phillips Tommaso Gili Manuel Zimmer Hernán A. Makse


#1#2Cline#1#2nil
Cline#1-#2#3nil
  multicnt#1
  multispan@ne
  multicnt=nefirstofone   multicnt#2
  multicnt-#1
  multispanne
  height#3








  
  
  
  
   

   Symmetries and synchronization from whole-neural activity in C. elegans connectome: Integration of functional and structural networks
 
 







 





Bryant Avila1*,
Pedro Augusto2,3*,
David Phillips 4,
Tommaso Gili 5,
Manuel Zimmer 2†,
Hernán A. Makse1,6,7†

1 Levich Institute and Physics Department, City College  of New York, New York, NY 10031, USA

2 Department of Neuroscience and Developmental Biology, University of Vienna, Vienna Biocenter (VBC), Vienna, Austria

3 Vienna Biocenter PhD Program, Doctoral School of the University of Vienna and Medical University of Vienna, Vienna, Austria

4 Mechanical Engineering Department, University of New Mexico, Albuquerque, NM
87131, USA

5 Networks Unit, IMT Scuola Alti Studi Lucca, Piazza San Francesco 15, 55100, Lucca, Italy

6 Department of Radiology, Neuroradiology Service, Memorial Sloan Kettering Cancer Center, New York, NY 10065, USA

7 CUNY Neuroscience, Graduate Center, City University of New
York, New York, NY 10031, USA



* Equal contribution 
† Corresponding author, manuel.zimmer@univie.ac.at  
†Corresponding author, hmakse@ccny.cuny.edu 









§ ABSTRACT

Understanding the dynamical behavior of complex systems from their
underlying network architectures is a long-standing question in
complexity theory. Therefore, many metrics have been devised to
extract network features like motifs, centrality, and modularity
measures. It has previously been proposed that network symmetries are
of particular importance since they are expected to underly the
synchronization of a system's units, which is ubiquitously observed in
nervous system activity patterns. However, perfectly symmetrical
structures are difficult to assess in noisy measurements of biological
systems, like neuronal connectomes. Here, we devise a principled
method to infer network symmetries from combined connectome and
neuronal activity data.  Using nervous system-wide population activity
recordings of the C.elegans backward locomotor system, we
infer structures in the connectome called fibration symmetries, which
can explain which group of neurons synchronize their activity. Our
analysis suggests functional building blocks in the animal's motor
periphery, providing new testable hypotheses on how descending
interneuron circuits communicate with the motor periphery to control behavior.  Our approach opens a new door
to exploring the structure-function relations in other complex
systems, like the nervous systems of larger animals.



§ SIGNIFICANCE STATEMENT

Complex biological networks exhibit an intricate relationship between the structural composition of their elements’ interaction and the whole system's functionality. This is particularly true in neuroscience, where the multiscale organization of neurons gives place to a multitude of hierarchical functionalities responsible for the nervous system interacting with the environment in which it is embedded. In this study, we focus on the chemical synaptic network of the backwards-crawling locomotion circuitry of the worm C. elegans to correlate structural clusters of neurons with dynamic patterns of their synchronization. Our findings reveal functional building blocks in the animal's motor periphery, providing new testable hypotheses on how descending interneuron circuits communicate with the motor periphery to control behavior.




§ INTRODUCTION
 

Complex systems depend on an intricate interplay between the dynamical
properties of their building blocks and the network structures that
link them together. For decades, researchers have endeavored to
uncover how the topological attributes of a network could shed light
on its dynamical behavior. While several methods have been proposed to
elucidate this interplay
<cit.>,  
the role of the network's symmetries has garnered considerable
attention in recent years
<cit.>. Identifying
and understanding such symmetries has potentially profound
implications: network theory implies that certain symmetrical
structures termed fibration symmetries in a graph permit
synchronizations
<cit.>,
a ubiquitously observed phenomenon in brain networks
<cit.>. Thus, identifying these symmetries might be
crucial for decoding some of the structure-function relationships in a
neuronal network.

Fibration symmetry refers to a specific type of symmetry where nodes
can be grouped into equivalence classes of balanced colorings based on
their input patterns
<cit.>. Nodes in the same
class receives identical 'input trees', which are hierarchical mappings
of input connectivity patterns, are proposed to share the same dynamics, and thus can
synchronize their behavior. These equivalent nodes are said to belong
to the same 'fiber' and are also 'balanced colorings' of the graph
<cit.>. This is formalized through admissible
ordinary differential equations (ODEs), which describe the evolution
of the system's state <cit.>. Theory ensures that
balanced colored nodes with isomorphic input trees following a set of
admissible ODEs will evolve in synchrony
<cit.>
forming a cluster of synchronization
<cit.>. This synchronization
can be crucial in various biological systems ranging from gene
regulatory networks to brain networks involved in language processing
<cit.>. Thus, identifying
fibration symmetries provide a prediction in which subsets of nodes
potentially exhibit synchronized dynamics.

However, there are technical and conceptual challenges to this
theory. Biological networks like connectomes are inherently
incomplete, and measuring them is limited by annotation errors, 
missing links, and noise, yet identifying fiber symmetries is sensitive
to small variations in a connectivity matrix. More fundamentally,
assuming mathematically perfect symmetries and equivalences of nodes
in a biological network is unrealistic. Here, we propose that network
symmetries represent constraints rather than perfect blueprints for
their implementations in individual animals. We aim to address these
challenges by studying the interplay of symmetry and synchronization
in the locomotor system of the nematode C.elegans.


The nematode worm C.elegans <cit.>exhibits a small nervous system
of just 302 neurons that develop from a stereotypic cell lineage into
118 anatomically and genetically defined cell classes, most of which
are comprised of a bilaterally symmetrical pair of neurons (e.g., AVAL
and AVAR) <cit.>.  Its connectome has been fully
reconstructed with a synaptic resolution
<cit.>, and it is
tractable for large-scale single-cell resolution neuronal calcium
imaging <cit.>.  Both, in
combination, offer unique opportunities to study structure-function
relationships in the nervous system <cit.>.  The worm connectome shares many
features with connectivity data from larger nervous systems, e.g.,
rich-club architecture and modularity 
<cit.>. 
Moreover,
neuronal recordings in C. elegans revealed nervous
system-wide neuronal activity dynamics that exhibit synchronization
patterns among well-defined ensembles of
neurons<cit.>. These
dynamics correspond to motor commands for a set of actions like
forward-crawling, backward-crawling, or
turning<cit.>.

These dynamics are generated in the absence of any known acute
time-varying sensory stimuli
<cit.> and can be
promoted by arousing conditions, such as environmental oxygen or blue
light <cit.>. Moreover, these
dynamics can be observed in immobilized animals; even under such
conditions, they extend to the motor periphery, incorporating motor
neurons, under unconstrained conditions mediating movement execution
<cit.>. These features point toward intrinsic
mechanisms that drive and maintain synchronous activity states. In
accordance, our previous work indicated that rich club architecture
and input similarities are crucial architectural features in the
connectome that permits synchronized neuronal
dynamics<cit.>.


In the present study, we focus on the chemical synaptic network of the
backward-crawling locomotion circuitry of C. elegans to circumvent some of the technical and conceptual constraints posed
by theory and computational limits.  This sub-network of just 21 individual
neurons includes three bilateral pairs of interneurons, AVA, AVE, and AVD,
which are situated in the head of the animal and send descending
chemical synapses to two classes of motor neurons, termed DA1-9 and
VA1-12. These motorneurons form neuromuscular junctions with the
dorsal and ventral body wall muscles, respectively, and are required
to generate the backward crawling gait <cit.>. All
chemical synapses in this sub-network are cholinergic and, with few
exceptions, neurons within the same class share similar morphology and
gene expression patterns <cit.>. Hence, assuming equivalence among nodes within
the same class could be a reasonable simplification.


By leveraging advanced calcium imaging
techniques, we first study neuronal dynamics, aiming to identify
patterns of synchronization within its backward locomotion circuitry. We
then use the observed synchronization dynamics to infer underlying symmetries in the
connectome.

Given the potential variability in connectomes across individual worms,
expecting the universal connectome model to predict the exact neural
synchronization across neurons is unrealistic. This variability
underscores the necessity of adopting tailored approaches to
understand the dynamic behavior from the structure of the
connectome. We approach this problem by reconstructing the underlying
connectome of C. elegans guided by the synchronization dynamics
obtained experimentally. This is done by "repairing" the connectivity
structure of a typical connectome <cit.> by
solving an appropriately modeled integer linear program to
find the minimal number of repairs of the available connectome needed
to ensure the observed synchronization dynamics. We aim for a
fibration-symmetric connectome that reflects the observed
synchronization. This method provides a means to idealize connectomes conceptually.  According to inter-animal connectome differences reported experimentally in
<cit.>, 
the threshold for connectome modifications can be taken to not exceed approximately 50%. However, for our final repair solution we take the stringent threshold of 25% given by animal-to-animal differences found in <cit.>.

To ascertain the robustness and validity of our solutions, we engage
in a rigorous testing procedure. We aim to determine the optimality of
our solution by reshuffling the neuronal labels in the connectome and
revisiting the repair process.  We accept our method's solution only
if the post-relabeling solutions consistently underperform our primary
findings (in over 95% of instances).

Our study bridges the gap between structural and functional neuronal
data, leveraging the power of symmetries in graph theory to shed light
on synchronization dynamics within C. elegans and the
structure of the connectome underlying this synchronization.











§ RESULTS
 


In this section, we delineate the pipeline employed throughout our
study in a condensed fashion and details are elaborated in the
Methods section <ref>.  Figure <ref> provides an
outline of the method to reconstruct the connectome:



  * The
process begins with the experimental setup to record neuron activity in  Fig. <ref>A with a microfluidic device
for immobilization of C. elegans enabling nervous system-wide Ca++ imaging with confocal
fluorescence microscope setup.


  * Time series data 
for activity traces of multiple neurons are recorded simultaneously, as
 shown in Fig. <ref>B. This
procedure is performed on multiple worms under similar
conditions. 


  * To obtain functional pairwise synchrony matrices, various metrics that capture synchrony are applied to the
  recorded time series data(Fig. <ref>C).


  * The synchrony matrices shown in Fig. <ref>D are
  averaged across worms to account for trial-to-trial and/or inter-individual variability. 


  * A standard thresholding process is then applied to obtain the
  functional network (Fig. <ref>E)
  <cit.>. Starting from a
  disconnected graph, we add links between nodes in decreasing order
  of weight of the averaged synchrony matrix. This is done for each
  averaged matrix.


  * Finally, a consensus matrix is calculated
  <cit.> across different methods of synchrony
  (Fig. <ref>F).


  * A hierarchical clustering algorithm is implemented to find a
  partition of synchrony clusters (Fig. <ref>G). Each
  neuron is assigned a color according to its cluster of
  synchrony. These colors are used as inputs to the repair algorithm
  in the next step.


  * We formulate and solve a mixed integer linear program (MILP)
  that finds the minimum number of edges to add or remove from the
  raw connectome (taken from
  <cit.>) to produce an ideal
  fibration-symmetric network that reproduces the coloring in the
  consensus partition obtained in the previous step
  (Fig. <ref>H).


  * The network produced by the solution to the MILP is a fibration symmetric connectome with cluster synchronization that reproduces the experimental data. This network can be collapsed into a smaller representation, the base
  graph, where nodes belonging to the same fiber have isomorphic input
  trees (Fig. <ref>I). The base graph suggests functional building blocks in the connectome.


  * For each consensus
  partition, a permutation p-value test is performed by permuting
  node labels and repairing the structural network 1,000 times. The partition with the lowest p-value is
  chosen as the optimal solution (Fig. <ref>J).
















 §.§ Neuronal activity recordings
 

To simultaneously record calcium activity from interneurons
in the head and motor neurons along the entire ventral cord, we
modify a recently reported whole nervous system imaging pipeline
<cit.>. Briefly, worms expressing the calcium probe
NLS-GCaMP6f in a pan-neuronal fashion, together with NeuroPal
multi-color cell identification labels <cit.>, are
immobilized in a microfluidic device that positions them in the field
of view of a spinning disk confocal microscope. After fast volumetric
GCaMP6f imaging, multi-color stacks are taken for cell class
identification (see Methods section 
<ref> for details).







 §.§ Whole nervous system recording and characterization of motor neuron activity
 

We generated eight datasets from different young adults well-fed
worms. Imaging covered almost the entire body spanning from the head
ganglia, the complete ventral cord, and the tail ganglia with
single-cell resolution (Fig. <ref>A-B). Animals were
recorded for 10 minutes at approximately three volumes per
second. Figure <ref>C shows a multi-neuron time series
with discernible calcium activity patterns of the neurons selected for
this study. These include neurons AVAL/R and AVEL/R, which belong to
the major descending interneurons conveying motor commands for
backward crawling, as well as the downstream backward crawling motor
neurons of the DA and VA classes. Consistent with previous studies
<cit.>, AVA and AVE
activity exhibited discrete transitions in their activity patterns,
characterized as low, rise, high, and fall states. We previously
validated that these states correspond to behavioral states in freely
crawling animals, i.e., low corresponds to forward crawling, rise and
high corresponds to backward crawling, and fall corresponds to turning
<cit.>.



Here, we observe that these features are mirrored in the motor
periphery (Fig. <ref>C). The mean pairwise correlations
among motorneurons and interneurons were high, indicating the strong
coupling between them (Fig. <ref>D). During free backward
crawling, animals generate a posterior-to-anterior traveling body wave
of alternating dorsal-ventral body oscillations. We did not observe
any obvious oscillations or activity patterns alternating between the
DA and VA motorneurons, which innervate the dorsal and ventral body
wall muscles, respectively. Neither did we observe any delay or
sequential activation patterns between motor neurons that innervate
more posterior to anterior muscles (which is indicated by the index
number in their name, n=1 most anterior, n=12 most
posterior). However, all motor neurons follow the rise and fall states
of the AVA descending command interneuron
(Fig. <ref>E). In conclusion, under our experimental
immobilization conditions, A-class motorneurons do not generate
gait-related patterns in their Ca++ activity profiles. However, they appear to respond to descending
motor commands from their presynaptic interneurons reliably. This
feature and the absence of movement, hence the lack of proprioceptive
inputs, makes the present setup particularly suitable for studying
neuronal dynamics dependent on internal neuronal circuit interactions
(see also discussion).





 §.§ Synchronization measures
 

Assessing the synchronization of two or more signals recorded
simultaneously involves various methods tailored to the specific
characteristics and similarities of interest. This results in a
diverse array of techniques aiming to capture synchronization through
different approaches <cit.>. These techniques
are generally classified into four primary groups based on their
focus: time domain versus frequency domain and methods that account
for directional dependencies versus those that do not (see Methods
section <ref> and  
Table
<ref>).










Correlation is a robust
measurement for coherence in neuronal data obtained via calcium
imaging
<cit.>
(Methods section <ref>).






In addition to correlation, to quantify the amount of synchronicity
between two signals, we implement the Level of Synchronicity (LoS)
measure introduced in <cit.>.
This metric determines how closely two
signals are to perfect synchronicity, meaning that
two signals have the same value at the same time, forming a cluster of synchrony:


    V_i(t) = V_j(t)           ∀ i,j∈ C_k.

Here, V_i is the dynamic of a neuron, which, in the present case, is
associated with the neuron's calcium activity. C_k is one of the
non-overlapping sets of neurons partitioning the neural system into
Cluster of Synchronicity (CS) <cit.>.

To quantitatively capture this synchronicity, the LoS evaluates the 
synchronization of two signals over time by considering their instantaneous differences 
and scaling them through a parameter σ defined as:


    LoS_ij = 1/T∑_t^Texp-[V_i(t) - V_j(t)]^2/2σ^2,

where T is the total amount of time steps in which
LoS is measured between the signals of neurons i and j.





The parameter σ serves as a scale to define a benchmark for
closeness between two points in time.

This parameter permits the user to deal with natural variations in
biological systems. No two neurons are perfect copies of each other;
therefore, if identical signals stimulate two similar neurons at rest,
their outputs (membrane potential) will naturally vary in intensity
and phase by small amounts. With this concept in mind, we implement
multiple versions of the LoS, each with a different value for
σ. When the value of σ is zero, each signal will only be
synchronous with its exact copy; as the value of σ increments,
it reaches the state that all signals are synchronous with all
others. This indicates that an ideal value lies between zero and an
upper limit (the largest difference between signals at some time t
should suffice). After many trials, we settle for values ranging from
0.01 up to 0.20

to explore different levels of synchronicity. We proceed with two
classes of synchronization metrics, correlations and LoS, each having
variations through different parameters or relational measurements,
making for a total of 44 metrics.



Matrices of synchronization and correlation are computed for each of
the N=8 worms under study, capturing the synchronous activities
observed within each individual as depicted in
Fig. <ref>C. To obtain a representative overview of
synchronicity across the entire set of worms, these individual
matrices are averaged per type, i.e. for each value of σ
(Fig. <ref>D). However, compiling these averaged matrices
require care, particularly in instances where the activity of
individual neurons are missing for specific worms. Therefore, we
calculate the average by averaging over worms for a particular metric
type and divide each element (neuron-pair functional measurement) in
the summed matrix by the number of times they appear together across
the whole cohort of worms.




 §.§ Identifying  cluster synchronization from the functional network
 


Each element in the averaged matrix contains information on a
neuron pair synchronous relationships used to construct
functional networks composed of all the neurons for the backward
locomotion gait. 








We follow standard thresholding procedures
<cit.> to build
the functional network.  Through this method, we purge the functional
data to only contain the strongest links leading to a functional
network, as seen in Fig. <ref>E.  Using this network, we
identify neuron cliques that are synchronized via community and
cluster detection algorithms.

A group of neurons that are more synchronous to each other
than other neurons belonging to different groups is considered
equivalent to a CS of neurons. We apply cluster and community detection algorithms to extract these functionally synchronous clusters of neurons from the functional network using two methods:


  
  
  * Clique synchronization <cit.>: It decides if a node belongs to a synchronous clique if the average value of the edges of the functional network inside the clique is bigger than any of
    its outside edges (see Methods section <ref>). 
    
    
    
  
  * Louvain community detection <cit.>: It decides the number of clusters and the number of
    nodes in each cluster by optimizing the modularity function as given
    by Eq. (<ref>) in the Methods section <ref>, which is solely cluster dependent. 


The result of applying these methods to a LoS correlation matrix
is seen in Figs. <ref>A-B, which results in a functional network with the node partition as seen in Fig. <ref>C.
The nodes with the same colors belong to the same cluster synchronization, as can be seen since they share thicker edges than
edges between nodes of different colors. 


The clusters obtained using the LoS measure at
different parametric values produce various motor-neuron partitionings. We only keep those that are unique among
the 39 results for the range of σ between 0.01 and 0.20. This
is done for each of the two cluster/community detection methods. The
number of unique partitionings for each method can be seen in the table
within Fig. <ref>E.  


















 §.§ Consensus synchronous clusters
 

We construct 44 metrics to obtain synchronization information: 39
LoS matrices and five types of correlation measures (Pearson,
Spearman, and Kendall coefficient, distance correlation and
covariance). We combine these metrics with the two functional
clustering methods mentioned in the previous section
<ref> (clique synchronization and Louvain) to
create 88 possible partitionings. Then, a consensus is created
following the methods developed in <cit.> to leverage
the information obtained from all the partitionings. To achieve this,
each partitioning, which consists of n clusters of synchronous
neurons, is used to create a co-occurrence matrix, which is simply a
matrix with n diagonal blocks, one for each synchronous group, with
value 1, and with all other values outside these blocks set to zero
(see Fig. <ref>D).  These matrices are summed up and
normalized.  The resulting consensus matrix is depicted in
Fig. <ref>G.  Through this, the sets of neurons
that appear in the same synchronous group more frequently than with
other neurons appear with higher values. This consensus avoids putting
each partitioning in competition with the other and instead compiles
them into a globally agreed-upon result <cit.>.

The consensus matrix, called X (Fig. <ref>G),
has optimal leaf ordering of its rows and columns. Such is produced by
the hierarchical clustering Ward metric (also known as Ward's
minimum variance method), applied to the dissimilarity matrix 1-X
(as seen in Fig. <ref>F). This metric aims to
minimize the total within-cluster variance. The goal is to choose the
successive clustering steps to minimize the increase in the total
within-cluster variance. This metric is particularly effective for
creating clusters that are compact and have a roughly similar number
of elements <cit.>, which was our main reason for
selecting it.

From this consensus, we can obtain various partitionings depending on
where the dendrogram is sliced.  When the dendrogram is sliced at a
value of 1.00 we obtain three major groups as observed in the thick
block structure in Fig.  <ref>G.  Each block can
be further divided into smaller diagonal blocks if we reduced the
threshold. For instance, 7 clusters are observed for a cutoff at 0.55
represented by the smaller blocks in Fig. <ref>G.
In the next sections, we will see that this partitioning is the optimal
solution obtained by the reconstruction algorithm with the least amount of
modifications to convert the baseline Varshney connectome into a
fibration symmetric solution (7 clusters for a cutoff at 0.55 in Table
<ref>).






All these partitions are the input to the mixed integer linear
programming algorithm that we designed to reconstruct the connectome
with a minimal number of modifications, as explained next.






  §.§.§ Balanced coloring partitions, fibrations and
  cluster synchronization



To develop the optimization algorithm to reconstruct the
connectome to satisfy the synchronization found in the previous
section, we need to introduce the concept of 'balanced coloring' and
the fibration of the graph.  Figure <ref>A shows an
example of a graph with a balanced coloring.

Let 𝒞 = C_1, ··· ,C_K be a partition of
the nodes of a network G=(V,E) with V vertices and E edges. We
identify each cluster C_k with a different color and K is the
total number of colors. A 'balanced coloring' is a coloring of the
graph such that each node with color k in cluster C_k is connected
(by edges of the same type) to the same number of nodes with color j
in cluster C_j, for 1≤ k, j ≤ K.  That is, nodes of the same
color receive the same colors from their neighbors
(Fig. <ref>A, center)

Translated into the terminology of dynamical systems, a system of
differential equations for the state variables of each node V_i(t)
can be interpreted as a 'message passing' process of passing colored
messages through the edges of the graph. Since two nodes i and j
with the same color received the same colors (messages) from their
neighbors, we can think, intuitively, that they will synchronize their
activity V_i(t) = V_i(t), forming a synchronous cluster. This
intuition is made mathematically rigorous <cit.>
through the theory of fibrations
<cit.>.

Traditionally, a way to formalize the balanced coloring partitions is
through the automorphisms of the graph forming its symmetry group
<cit.>. In this case, the
orbits of the automorphisms of the graph are the balanced
colors. However, not all balanced colorings are orbits. Many
biological networks contain no automorphisms, yet, they display a
non-trivial balanced coloring partition with many colors
<cit.>.  This is exemplified in
Fig. <ref>A. The graph has no automorphisms (except for
the trivial identity), but has a balanced coloring that reduces the
graph to a base with just two nodes. This balanced coloring is only
captured by the fibration.


The graph fibration formalism introduced in <cit.>
is a more general formalism that captures all balanced coloring
partitions of the graph. Graph fibrations are a particular case of
fibrations between categories introduced previously by Grothendieck
<cit.> and are formal generalizations of
graph automorphisms.

A graph fibration (shown on the right of Fig. <ref>A) is
a morphism of the graph that collapses every cluster of synchronized
balanced colors (called 'fibers') into a single representative node in
the base graph B while conserving the 'lifting property' defined in
the Methods section <ref>.  This transformation leaves
invariant the dynamics in the graph and captures the maximal
symmetries of the network. It is then called a symmetry fibration
<cit.>. Here, we propose that fibers represent potential functional modules or building blocks in a network of neurons.




























 §.§ Fibration symmetry driven repair algorithm
 

The clusters of synchronous neurons found in section
<ref> are used to repair the raw chemical
synapses connectome of the backward locomotion gait of
C.elegans provided by Varshney et al.
<cit.> (plotted in
Fig. <ref>B). We develop an optimization mixed
linear integer program to modify this network by adding or removing
the least amount of edges within the limits of an objective function
<cit.> to match the synchronization obtained
experimentally. We call this algorithm the Symmetry-Driven Repair
Algorithm or SymRep for short.

The goal is to obtain a network with balanced coloring partitioning as
provided by one of the synchronization clusters from
Fig. <ref> (details appear in the Methods section
<ref>). The most important feature of this algorithm in the
context of the present paper is the objective function used to
determine optimal solutions. Many combinations of adding and/or
removing edges can satisfy turning the connectome into a network with
color partitionings as provided by a functional clustering while
respecting the constraints provided by Eq. (<ref>) and
Eq. (<ref>). Of these many solutions, the one that
minimizes the objective function below provides an optimal solution
(Methods section <ref>):

    
    f_α,β(r,a)=α∑_i,j∈ E r_ij + β∑_i,j∈ E^C a_ij
   
where E is the set of edges present in the original connectome and
E^C is the set of non-existing edges in the original connectome that
are permissible to be added. The terms r_ij and a_ij are
binary variables, where the first term, associated with E, takes on
a value of 1 if an edge has been removed. The second term, associated
with E^C, takes on a value of 1 if an edge has been added. If
necessary, one can control the SymRep algorithm to prohibit removing
(adding) connections between physically impossible neurons, which
studies show are always connected (disconnected). This can be done by
simply not including them in E or E^C. The parameters α and
β are penalty weight constants for each of these variables,
respectively. The relative weight between these parameters determines
if the SymRep prefers to repair a network through the addition of
edges rather than removal. In our study, we explore the solutions
found by varying the relative differences between α and
β. We keep the value of α at one and increment the value
of β in steps of one, starting from one and culminating at 10.

We work with the raw chemical synapses connectome of
C.elegans for the backward locomotion gait provided in
<cit.>. This is composed of the connectivity
between motorneurons (VAs and DAs) and 3 main interneuron pairs
(AVAL/R, AVEL/R, AVDL/R). We specifically use two versions of this
connectome.  One is that reported by Varshney et al.
<cit.> shown in
Fig. <ref>B.

We call it the Uncollapsed Varshney connectome. The second, shown in
Fig. <ref>C, is a collapsed version of this
network in which we assume that motorneurons can not distinguish
between signals received from the left or right version of an
interneuron pair (i.e., AVAL or AVAR). This is done by substituting an
interneuron pair for one node (i.e., AVAL and AVAR are collapsed into
AVA) and substituting the two edges a neuron may receive from an
interneuron pair for one edge with a weight equal to 1 if both edges
are present, one if only one edge is present, and 0 if not
present. This version is called the Collapsed Varshney connectome. As
a final note, we assume each interneuron pair belongs to its unique
synchronization cluster.



 §.§ Synchronization-driven reconstruction of the connectome
 

With the 39 unique clusterings found in the previous section, the 19
different combinations of weight penalties and the two different raw
connectomes (original and collapsed versions) used, we run SymRep to
produce 1,482 solutions. All of these solutions satisfy the condition
of balanced coloring. Of the 741 solutions for the collapsed Varshney
connectome, 730 (98.52%) satisfy the stricter condition of minimal
balanced coloring, and 737 (99.46%) fulfill this condition for the
not collapsed version. Of these, we present the most optimal solution
for each type of cluster or community detection algorithm in
combination with the functional measurement group (LoS and
correlations).

We decide on an optimal solution based on the number of modifications
it enacts on the raw connectome with the condition that this number is
the lowest, and at least below  50% of natural variability found from
animal to animal <cit.> to accept the
solution. We pick those that modify the connectome by the least amount
of addition and removal of edges. We consider two cases to measure the
number of modifications. The first one is where adding or removing has
a cost of 1, where the final sum of modifications is divided by the
total number of edges in the raw connectome. The second case is in
which the severity of  
removing an edge is tied to the weight of the
edge in the collapsed connectome before binarization while adding an edge only has
a cost of one.  In the latter, the total sum of modifications
(additions + removals) is divided by the total number of edges in the
raw connectome. Metadata about the solutions can be observed in Table
<ref> with significant p-values for each method.










 §.§ Statistical significance


To test the significance of these functional partitionings and
eventually choose the most significant partition, we proceed to use a
permutation test <cit.>. Specifically, for each of
the partitionings considered in Figure
<ref> we randomly permute the labels
(names) of the neurons, keeping the number of clusters and number of
nodes distributed among these the same. This leads to a shuffling of
synchronous groups keeping their size (amount of neurons in each) the
same. We produce slightly over 1,000 permuted versions of each of
these partitionings and inspect how many lead to networks with these
partitionings as fibers with equal or less number of edge
modifications. This quantity is then divided by 1,000 or by the number
of networks that lead to a minimally balanced coloring solution with
an equal number of fibers as in the partitionings (∼1,000).












 



 §.§ The optimal solution



The consensus matrix allows the detection of fibers at different
scales (number of clusters) through the information contained in its
hierarchical clustering, as seen in
Fig. <ref>F. These fibers are used for repairing
the Varshney chemical synaptic backward gait produces some outstanding
results among the collapsed and non-collapsed connectome. We highlight
the results for the collapsed connectome at the bottom portion of
Figure <ref>.

When it comes to the statistical significance of clusters obtained
through a system of ever-increasing node partitionings to repair a
network, as in our case, there is some expected behavior. For
instance, the trivial partitioning of N neurons into N
non-overlapping clusters will give a p-value of 1 as all permutations
of a trivial partitioning will produce the same partitioning,
resulting in all of them modifying the network by the same amount. The
same can be said for a trivial partitioning with 1 cluster as all
permutations result in the same partitioning. Therefore, as the number
of clusters change between 1 to N, the p-value of statistical
significance is granted to produce a (global) minimum p-value between
1 and N clusters <cit.>.  For a low-noise system
in which node dynamics are faithful to the network structure with
1≤ M≤ N fiber partitionings, it is suggested that as the number
of non-overlapping clusters is reduced from N, the p-value is
expected to improve (be below 1), with the same expected p-value as
the number of clusters is increased above one where the minimum
p-value should be located at M clusters. Deviations from this
depends on how far removed the network under repair is from the
original fiber symmetric version producer of the node mentioned above
dynamics. These predictions in the behavior of p-value are observed in
the repairing of the backward locomotion network in
Fig. <ref>. As the number of clusters
approaches the trivial N-partitioning the amount of repairs needed
decreases because the network we are trying to repair is indeed
originally trivially colored.

Applying all these considerations, we find that the most optimal
partition and reconstructed connectome is obtained for 7 clusters at
consensus clustering cutoff of 0.55 repaired from the collapsed
Varshney connectome.  As shown in
Fig. <ref> and Table
<ref> this solution has the lowest
p-value=0.033, being the only one of the consensus solutions below the
standard 0.05 p-value cutoff for statistical acceptance. Furthermore,
the repair percentage is 21.88%, below the acceptance limit of 25%
variability from animal to animal. Therefore, we choose this
connectome, shown in Fig. <ref>A, as the one
representing the data in the closest possible way.



§ DISCUSSION





Understanding how the structure of the connectome influences the
function of the network is a long-standing problem in system
neuroscience <cit.>. However, addressing this
question is challenging due to low sample number in available connectome datasets and potential 
variability across individuals.  We have addressed the problem of incomplete/missing data in
biological networks by first developing a consensus method to obtain
average clusters of neuron synchrony across worms from whole-body
calcium recording of neural activity. This information drives a
reconstruction optimization algorithm of the connectome consistent with the observed synchronization.  When the modifications to
the connectome are below the experimentally found variation from
animal to animal and offer a
minimal p-value significance, the obtained connectome models can be thought
of as idealized networks consistent with the experimentally measured dynamics of the neuronal circuits. Such a network bridges the gap between structure
and synchronization and can be used to further assess the
functionality of the connectome via perturbations to its structure.

Our research has substantiated the pivotal role of fibration
symmetries in the connectome structure in orchestrating neuronal
synchronization in C. elegans. The refined understanding that
structural connectomics underpins significant aspects of neural
functionality furthers our grasp on the physical bases of neural
synchronization. This synchronization is crucial, as it forms the
foundation for coordinated motor outputs and behavioral responses in
organisms. The use of advanced calcium imaging and graph theory to
correlate these structural motifs with dynamic patterns of neural
activity offers a compelling model for predicting neuronal behavior
based on underlying anatomical data.


The methodological innovations in our study, particularly the use of
integer linear programming (SymRep) to adjust connectomics data based
on functional imaging, demonstrate a novel computational approach to
neuroscientific research. This technique refines existing neural
network models and provides a quantifiable method for aligning
theoretical predictions with empirical observations.

Previous research in C. elegans has suggested local oscillators at the motor
neuron level for rhythm generation in the VNC for reverse
locomotion <cit.>; these patterns have been observed when motor neurons were experimentally deprived of their descending interneuron inputs. Our experimental conditions, Ca++ imaging in immobilized worms with interneurons AVA and AVE generating descending motor commands, did not reveal any oscillations or activation patterns consistent with a rhythmic backward gait pattern i.e., alternations between VA- and DA-motorneurons, or a traveling wave from the posterior to anterior position. Instead, we observed that A-motor neurons faithfully mirror the descending motor command activity in sync. It was shown that C. elegans can adapt their gait flexibly to the physical properties of the environment via proprioception, as shown for forward locomotion <cit.>, and that A-motorneurons have indeed proprioceptive properties themselves  <cit.>. We suggest that proprioception might be an essential driver to unleash the oscillatory or rhythmic properties observed in ref.  <cit.>. In our study, the lack of proprioceptive inputs and absence of complex internal motorneuron dynamics was an advantage, since it enabled us to study the activity of all A-motorneurons solely in response to their synaptic partners, providing a filter on activity that should be mostly explainable by the connectome architecture alone.   

Connectomes between adult worms are remarkably variable, in particular for connections between neurons with few synaptic contacts <cit.> as well as in the posterior ventral cord (VNC) <cit.>. Extrapolating from these data, connections between descending interneurons and motor neurons in the entire VNC might vary by up to 50% i.e., about half of the connections found might not be reproducible across individuals. Our connectome repair algorithm provides a family of statistically significant solutions that require connection additions/removals well below this range (Table 1, Fig. 6). All of these solutions provide circuit models with fiber symmetrical properties consistent with synchronization patterns found in A-motorneurons. It is unlikely, that individual animals indeed perfectly match these idealized symmetric solutions, thus neuronal networks must be endowed with additional properties that ensure stable synchronous dynamics. However, here we propose that fibration symmetries provide constraints, from which individuals should not deviate too much to ensure the robustness of synchronous dynamics. Here we show that symetrization procedures can, however, be useful, for identifying potentially functional circuit modules in the resulting base graph. 

Our analyses led to a base graph of the backward locomotion circuit that suggests functional modules characteristic of differential innervation by the descending interneurons AVA, AVE, and AVD (Fig. <ref>E), which intriguingly appear as a crude topographic map (Fig. <ref>A). Our analyses thereby suggest that these modules could operate as functional units in deferentially transmitting the motor commands to different body parts, perhaps in a manner supporting effective backward locomotion. AVA, AVE, and AVD differ to some degree in their connectivity to other sensory circuits <cit.>; we thus speculate that the modules might enable differential control over how to execute the reversal motor command. Note, that AVD neurons were not active in our recordings. Future studies that activate also AVD, and selective inhibition or interference with AVA, AVE, and AVD will test the relevance of our findings for the control of A-motorneuron activity and backward crawling behavior.   

The present study focused on a smaller neuronal circuit where some of the assumptions that our theory includes might be reasonable simplifications. Anatomically and molecularly, left and right members of interneuron classes AVA, AVE and AVD are nearly indistinguishable <cit.>, justifying bilaterally collapsing the connectome (Fig. <ref> and Table
<ref>). Moreover, most chemical synapses in the circuit are excitatory and cholinergic <cit.>, and VA and DA motorneurons have similar molecular footprints <cit.>, justifying the simplification of node and link equivalency. Future work on larger circuits in the worm and other model organisms should develop our theoretical approach further to realistically allow for more diversity in neuronal and synaptic properties.  



Indeed, the implications of our findings are not restricted to
C.elegans. The principles of neuronal synchronization,
facilitated by connectomic structures observed in this study, likely
have analogs across bilaterian species, including humans
<cit.>. Investigating these principles in more
complex nervous systems could reveal new insights into how brain-wide
synchronization patterns contribute to complex behaviors and cognitive
functions. Furthermore, exploring the impact of structural variations
on the synchronization in disease models could open new therapeutic
avenues, particularly for neurological disorders where dysregulation
of synchronized activity is evident.






§ METHODS
 

In this section, we first explain the experimental setup and the
methods to extract the functional synchronization clusters of neurons
obtained using various combinations of functional measurements and
different methods of extracting clusters and communities. Following
this, we introduce the graph fibration formalism to describe balanced
coloring partitions of graphs and cluster synchronization. We conclude
with the repair algorithm for optimal fibration symmetric network
solutions (SymRep) applied to the synchronization partitionings.






 §.§ Neuronal calcium imaging
 

Whole-nervous system Ca2+ imaging experiments were performed on
transgenic young adult C.elegans hermaphrodites (age was
determined by the number of a maximum of 5 eggs) expressing genetically-encoded
calcium indicator NLS-GCaMP6f in a pan-neuronal fashion and localized
to the cell nuclei <cit.>
together with NeuroPal cell class identification labels
<cit.> (ZIM2001: otIs669; MzmIs52 ; lite-1 (ce314))

Mounting and microfluidic setup: Animals were imaged in
two-layer PDMS microfluidic devices to control the oxygen environment
and with curved channels <cit.> to immobilize and
laterally align animals, enabling reliable positioning of worms across
recordings and fitting them into the field of view of the imaging
system. Such a layout allowed us to cover the head ganglia, ventral
cord, and tail ganglia of the animals.  The worm channel of the
microfluidic device was connected to a syringe that contains NGM
buffer with 1 mM tetramisole to paralyze worms. All components were
connected using Tygon tubing (0.02 in ID, 0.06 in OD; Norton) using
23G Luer-stub adapters (Intramedic). Constant gas flow of 21% O2 and
79% N2 (50ml/min) was delivered using a gas mixer connected to mass
flow controllers (Vögtling Instruments) using lab custom scripts in Micromanager. Adult worms were picked on food-free NGM agar plates in a
drop of NGM with 1mM tetramisole and aspirated into the worm channel.
Animals were first habituated for 10min and afterward imaged at 21%
O2 for 10 min.

Microscope setup: High-resolution data of neuronal activity in
the head and tail ganglia were acquired with an inverted spinning disk
confocal microscope (Zeiss Axio Observer.Z1 with attached Yokogawa
CSU-X1) using a sCMOS camera (pco.edge 4.2 with Camera Link HS connection) and a 40x 1.2 LD LCI Plan- Apochromat
water-immersion objective (Zeiss). Moreover, we used 0.5x demagnification relay optics at the port of the spinning disk unit to increase intensity/pixel and the total area that can be imaged on the camera sensor both by the factor of 2, allowing imaging of a larger field of view encompassing the full worm nervous system. We use a custom-made GUI in
Micromanager to control the different elements of the
microscope. Exposure time was 20 ms, with 2μm steps between Z-planes
operated by a Piezo stage (P-736 PInano, Physik Instrumente GmbH); the
total plane number of z-planes varied between 16-20, leading to a
volume acquisition rate of up to 3 Hz.


Neural trace extraction: As described in detail in our previous
work in Kato et al. <cit.>, neuronal activity
traces were obtained by tracking the intensity maxima in each volume
over time and calculating the single-cell fluorescence intensities. F0
was calculated for every neuron as the mean fluorescence intensity
across each trial. After background subtraction, DF/ F0 was calculated
for each neuron, following bleach correction by linear detrending and
exponential fitting. See reference <cit.> for details.

Neuronal identification: Identification of each neuron was done
based on a neuron dictionary termed NeuroPAL <cit.>
(Fig. <ref>A).  In each recording, we aim to detect at least 25 crucial neurons for our analysis: two reversal interneurons AVA and AVE, and all the reversal motor neurons DA01-DA09 and VA01-VA12.  Some neurons in the head and tail
tip were lost in individual recordings, depending on animal size, when it exceeded the imaging
region. Other reasons for missing neurons might be a failure in
segmentation and tracking due to low signal/noise ratios based on low
expression levels and/or low calcium levels. Neuronal traces were
curated after each recording, and obviously erroneous traces were
removed.
Average correlation matrices were generated by calculating the mean pairwise correlations between the activity time series of identified active neurons (up to n=8, depending on the number of pairwise observations, but at least n=3). The identified active neuron numbers are as follows: 72 for recording 1, 76 for recording 2, 98 for recording 3, 79 for recording 4 and 5, 69 for recording 6, 63 for recording 7, and 62 for recording 8


NeuroPAL labeling: The identities of neurons were determined via
NeuroPAL using the following procedure. We obtained image stacks from
each recorded animal after GCaMP6f imaging.  For each plane, we
acquired spectrally isolated images sequentially of CyOFP1,
mNeptune2.5, and mTagBFP2. We excited CyOFP1 using the 488nm laser at
5% intensity with a 585/40 bandpass emission filter. Afterward,
mNeptune2.5 and TagRFP-T were recorded using a 561nm laser at 20% intensity (percentage of max) with a 655LP filter and a 570LP filter, respectively. mTagBFP2 was isolated using a 405nm laser at 30% intensities with a 447/60 bandpass filter. These imaging
conditions allow for the color-specific labeling of
C. elegans neurons. Neuronal identities were manually annotated
according to the NeuroPal guidelines. <cit.>.



 §.§ The zoo of synchronization measures
  

Determining the synchronization between two or more simultaneously
recorded signals is a task with multiple approaches depending on the
characteristics that one is interested in measuring and determining
their level of similarity. Naturally, all this leads to a zoo of
methods that try to capture the synchronization between signals
through different features <cit.>. These
methods tend to be broken into four main camps of measurement, which
are the intersection of the time domain vs the frequency domain and
those that consider directional dependencies and that do not; some of
these examples can be observed in Table
<ref>.



Depending on the data and the aspects of the signals one is interested
in, one may use a few examples in Table
<ref> and even other methods not
mentioned in the table. Selecting the appropriate metrics can
sometimes be a difficult problem, and one must compare these against
some reference, producing competition between the selected
metrics. Below, we indicate which metrics we have decided to
use. Further, down in the paper, it is shown that it is best not to
put these metrics in competition one against the other to find the
best metric that captures synchronization. Instead, a better approach
is to use the information that these provide in a consensus manner to
decipher the synchronization in our system <cit.>.




 §.§ Clique synchronization and Louvain method
 


From section <ref>, the Clique Synchronization
method developed in <cit.> accepts a clique if all
nodes composing a clique satisfy the conditions outlined in section
<ref>. We define a cluster of N neurons to
be synchronous as a fully connected clique consisting of N nodes
that meet the conditions:

    ∑_i<j^1,Nσ(x_i (t),x_j(t)) ≥N(N-1)/2σ(x_k (t),x_k'(t))   ∀ k= 1, …, N  and  k' ∈ℳ_k,

where ℳ_k denotes the set of nearest neighbors of node k
(for k=1, …, N) that are not part of the clique in
question. Here, σ(x_i (t),x_j(t)) represents the value for
synchronization, measured by the LoS metric or any correlation
metric used in this paper.







The Louvain method <cit.> is a popular algorithm for
module and community detection. At every execution, it assigns every
neuron to its community, as per the method used here, it proceeds to
randomly group neurons into bigger clusters, only accepting mergers if
the modularity value given by Eq. (<ref>) increases and
halting once this equation can no longer increase in value
<cit.>. Due to the stochastic nature of this process,
each execution may lead to a slightly different result. Due to this,
the Louvain method is executed 1,000 times for a given network
retaining the partitioning with the highest measured modularity. The
modularity measure is defined as:
    
    Q=1/2m∑_i,j[A_ij-k_ik_j/2m]δ(c_i,c_j),

where A_ij is the adjacency matrix of the network, m is the
number of edges in the network, k_i is the number of in-degree
edges attached to node i and c_i are numeric labeled given to
each synchronization cluster <cit.>.




 §.§ Fibration symmetry of graphs and cluster synchronization

  

In the study of complex networks, understanding the interplay between
structure and dynamics play a crucial role. Specifically, the concept
of synchronization, where nodes in a network adjust their behavior in
accordance with each other is fundamental in various disciplines,
including physics, biology, and engineering. This section delves into
the intricate relationship between graph topology, characterized by
fibration symmetry, and the emergence of synchronization through the
lens of admissible ordinary differential equations (ODEs)
<cit.>.


A graph G=(V_G,E_G) consists of a set of vertices V_G and a set of
edges E_G where each edge connects a pair of nodes u to v
represented as (u,v) ≡ e^u → v_G. The set of edges
from node u to node v is written as E_G(u,v), that is, the set
of edges e ∈ E_G that have as a source node u brought by the
function s(e)=u and have as a target node v brought by the
function t(e)=v. The set of edges in graph G that have as a target
the node v is denoted by E_G(-,v) <cit.>.

A graph G can be mapped into another graph W through the structure
preserving process of a morphism Ψ:G ⟶ W using a
pair of functions Ψ_V:V_G⟶ V_W and
Ψ_E:E_G⟶ E_W. These functions map vertices to
vertices and edges to edges, respectively, from graph G to W,
while being commutative with the source and target mapping
functions. In mathematical terms, these functions must obey
<cit.>:

    s_W∘Ψ_E =
    Ψ_V∘ s_G ,

and

    t_W∘Ψ_E = Ψ_V∘
      t_G .


If the mapping functions Ψ_V and Ψ_E are surjective, that is
to say, they map multiple elements (vertices or edges) in their domain
(graph G) to one representative element in their range (graph W),
then morphism Ψ is called an epimorphism.



A (surjective) graph fibration is a particular type of epimorphic
morphism that maps elements that belong to the same category into one
representative element while preserving the lifting property which
guarantee dynamical invariance.  The original definition of fibration
was defined between categories by Grothendieck and others
<cit.>. Boldi and Vigna
<cit.> worked out the definition of fibrations
between graphs, which is the main theoretical framework of the present
work.  Morone et al.  <cit.> then showed the
application of fibration to biological networks to understand their
building blocks <cit.>, and synchronization
<cit.>.

A graph fibration is a morphism from
graph G to the base  B:

    φ: G ⟶ B

that satisfy the lifting property <cit.>.  This
means that for every edge e = (u,v) ∈ E_G, there is an edge e' =
(φ(u), φ(v)) ∈ E_B, and for every vertex v ∈ V_G there
exists a unique edge e^v∈ E_G(-,v) satisfying (w,φ(v))
∈ E_B.

The importance of the lifting property is that the dynamics between
the graphs G and B is preserved. This means that if we add a set
of admissible ODEs to the graph with dynamical variables x_i(t),
then the dynamical evolution of the system in graph G is the same as
the dynamical evolution of the graph B.

    (G) =   (B) .

This property defines the clusters of synchrony as follows.

When a function φ: G → B acts as a fibration, we
refer to G as the total graph and B as the base graph of this
mapping function. In this context, G is said to be fibred over
B. The set of vertices in G that φ sends to a specific
vertex x in B is called the fiber over x, denoted by
φ^-1(x). Fibers are also balanced color clusters.

Vertices belonging to the same fiber have isomorphic input trees as
defined in <cit.>. The input tree for a node v,
denoted T(v), is a rooted tree centered at node v that captures
all the paths in the graph leading to v.  The first layer of the
tree is the node's in-neighborhood, called its input set.  Each
subsequent layer is then iteratively defined as the input set's input
set. A visual example of an input tree can be seen in
Fig. <ref>I, where all blue nodes in the total space
graph have the same input tree structure.

In terms of dynamics, we can think of a message passing process where
the information travels through the input tree and arrives at the
rooted node v. If another node u has an isomorphic input tree with
v, meaning that T(v) ∼ T(u), then these two nodes receive the
same messages through the network, although from different
pathways. Therefore u and v synchronize their dynamics x_v(t) =
x_u(t). This statement has been put in rigorous mathematical terms by
DeVille and Lerman <cit.>.


A fiber is called trivial if it contains exactly one vertex, that is
if |φ^-1(x)| = 1 and called nontrivial if
|φ^-1(x)|>1. Throughout this paper, the function |·|
denotes the number of elements in the mathematical structure it
encloses unless specified otherwise.



Vertices belonging to the same fiber, say u and v, have an
equivalence relation ≃ called an in-isomorphism, which can be
seen as a particular case of a discrete homeomorphism
<cit.>, where a discrete bijective function
(one-to-one association of elements in two domains) ψ: G(-,u)
→ G(-,v) is applied to the discrete topological space of
graphs.

Additionally, the in- in in-isomorphism is the added
condition for the bijective function to obey s(e) ≃ s(ψ(e))
for all e ∈ G(-,u). In short, in-isomorphism holds the notion that
the vertices in a graph G can be converted into one another while
conversing the adjacency connectivity from the same (or equivalent)
source vertices.


A fiber in a graph G is associated with the cluster containing the
vertices of the fiber u, v, w,...∈ C_i where n_i represents the
number of vertices in cluster C_i with its vertices noted as
∪_l=1^n_i v^i_l = C_i. The union of these clusters
contains the entirety of vertices C_1 ∪ C_2, ∪...C_K∈ V_G
where K is the number of fibers in graph G. The overlap of
different clusters is always empty C_i∩ C_j=∅ for i≠
j  ∀ i,j∈ K_G where K_G={x∈ℤ | 1 ≤ x
≤ K } such that ∑_c=1^K n_c = |V_G|.

These fibers are the equitable partition or colored-balanced
partitioning of the graph as the number of in-degree edges of the
in-isomorphic vertices in a cluster C_i receives from another
cluster C_j only depends on the choice of the clusters.  Thus, both
fibers and balanced colorings describe, in two different ways, the
same synchronous dynamics of clusters of nodes.

We can represent this as:


    C_i↢ C_j≡{|∪_l=1^n_j E_G(v_l^j,
    v^i)|=|∪_l=1^n_j E_G(v_l^j, v^k)|}      ∀ v^i,v^k∈ C_i

where the number of edges |·| every vertex in C_i receives
from a cluster C_j must the same <cit.>.



The particular type of fibration we work with in this paper is the
minimal fibration, which leads to a minimal equitable
(colored balanced) graph partitioning. The word minimal enforces the
value K representing the number of colored balanced partitionings of
a graph to be the lowest it can be. This leads to the situation in
which no two clusters receive the same in-degree edges relative to a
third cluster, which leads to the constraint

    ∑_k=1^K |(C_i↢ C_k)-(C_j↢
      C_k)|>0      ∀ i,j∈ K_G.


The surjective minimal graph fibration is called a symmetry fibration
<cit.> since it collapses the graph into its
minimal base capturing the minimal number of fibers (or balanced
colors) of the graph, and thus collecting the maximal symmetries of
the graph.

An implementation of the algorithm to find minimal balanced colorings
in a graph in the form of an R package is available at
<https://github.com/makselab/fibrationSymmetries> and
<https://osf.io/z793h/>.



 §.§ Fiber symmetries constrain admissible ODEs into synchronization


The fibration symmetry of a graph imposes a structural constraint that
can facilitate synchronization. Specifically, if the graph has a
fibration symmetry, the nodes in each fiber can be expected to
synchronize with each other <cit.> due to the
consistent interaction patterns enforced by symmetry. This can be
understood through the synchronization of a system of ordinarily
differential equations (ODEs) 'admissible' to a graph
<cit.>.

Admissible means the structure of a graph G imposes the coupling
terms between the |V_G| number of ODEs, where each equation
characterizes the state of a vertex in the graph. Consider a dynamical
system on a graph where each vertex v_i has a state x_i that
evolves according to an ODE:

    dx_i/dt = F(x_i,t) + ∑_j ∈∂_i A_ijH(x_i,x_j,t) .


Here, F represents the intrinsic dynamics of each vertex, H
embodies the interaction between vertices, and A_ij are the
elements of the adjacency matrix A of the graph, where A_ij=1
indicates the presence of an edge and A_ij=0 indicates the absence
of an edge from vertex j to i. ∂_i denotes the set of neighbors
of vertex i. As an example, the Kuramoto model, taken along with its
master stability function can be used to study synchronization in a
network
<cit.>.

Cluster synchronization as defined by Pecora et al.
<cit.> in this context refers to the situation in
which the dynamical states of a group of vertices for a given fiber
converge to a common trajectory, i.e., x_i(t) →
c_ℓ(t)  ∀  v^ℓ_i∈ C_ℓ, given the
appropriate initial conditions and after any transients have died out.
























The symmetries in the graph fibration allow the system of ODEs
admissible to the graph G to be reduced from a |V_G| number of
equations to a K number of equations for each fiber
cluster. Equation (<ref>) can be reduced to a system of equations
for each of the fibers:

    dc_i/dt = F(c_i,t) + ∑_j ∈ V(c_i) Q_ijH(c_i,c_j,t).


Here, Q is the adjacency matrix of the base of the graph obtained by
the fibration.

The fibration formalism extends the automorphism symmetry groups of
the graph <cit.>. An automorphism is a permutation
symmetry of the nodes of the graph that leaves invariant the adjacency
of nodes. That is, the nodes permuted by the automorphism have the
same in- and out-neighbors before and after the application of the
automorphism. The analogous of fibers in fibrations are the orbits of
the automorphisms. An orbit a node is the set of nodes obtained by the
application of all the automorphisms of the graph.

Both orbits and fibers are balanced colorings of the graph. All orbits
are fibers, but not all fibers are orbits. Thus, the fibers of the
fibration capture more balanced colorings than the orbits of
automorphisms. In short, fibration symmetries are rigorous extensions
of automorphisms, and all automorphisms are fibration symmetries but
the opposite is not always true.

This situation is exemplified in Fig. <ref>A. This graph
has no automorphism. That is, there is no permutation of any nodes
that leaves invariant the adjacency matrix. Therefore it has a trivial
orbital (and balanced coloring) partition where each node is its own
orbit or its own color. However, there is another minimal balanced
coloring partition that is shown in the two colors in the center of
Fig. <ref>A. This is the balanced coloring partition
captured by the fibration, which is then applied to reduce the
network to the base on the left of Fig. <ref>A.

A final fundamental difference between automorphism and fibration
symmetry is that the former is a global symmetry of the network. That
is, it is a permutation that constrains the global adjacency of nodes
to be the same. However, the fibration is a local symmetry leaving
invariant the input trees of the nodes, which are local views of the
network of the nodes collapsed by the fibration.  This fundamental
property makes the fibration applicable to a larger set of complex
networks than the more restricted automorphism.





























The interplay between fibration symmetry and synchronization in
networks via admissible ODEs provide profound insights into the
dynamics of complex systems. The structural constraints imposed by
fibration symmetry can dictate the synchronization patterns, where the
quotation or base graph provides a simplified version of a larger
networks leading to predictable and potentially controllable behavior
in networks. This understanding has significant implications for the
design and analysis of complex systems in various domains, ranging
from biological networks (as in this paper) to engineered distributed
systems
<cit.>



 §.§ Symmetry Driven Repair Algorithm
 



In this section, we present a set of equations and inequalities based
on the concepts discussed previously to construct our integer linear
programming model for our problem, considering a specific graph, node
clustering, parameters, and constraints
<cit.>. We name this as the
Symmetry-Driven Repair Algorithm or SymRep for short.

For this model, we consider a directed graph where we denote n=|V| and m=|E| as the number of nodes and directed edges, respectively. We also define

    E^C = { (i,j) : i,j ∈ V, ij ∉E}

as the set of node pairs between which no directed edge exists in G: these pairs indicate potential edges that could be added to the graph G. We define  as a coloring of G,  represents the sets dividing V, i.e., the various clusters of nodes.

The three types of decision variables in the model are:




r_ij for (i,j) ∈ E     is a binary indicator for when edge
  (i,j) is removed; 
 a_ij for (i,j) ∈ E^C     is a binary
  indicator for when a new edge (i,j) is added; and 
 s_ijR for
  i ∈ P, j ∈ Q, R, P, Q ∈     is a binary indicator that i
  and j has an imbalance in the color set R.





The objective function aims to minimize the weighted sum of edges
removed and added, defined as

    f_α,β(r,a) = α∑_ij ∈ E r_ij + β∑_ij ∈ E^C a_ij.

Constants α, β are parameters that adjust the importance
between edge removal and addition in the objective.  The main
constraint ensures that  represents a balanced coloring of the
graph G, as defined by


    ∑_ip ∈ E: i ∈ S (1 - r_ip) + ∑_ip ∈ E^C:i ∈ S a_ip    = 
    ∑_iq ∈ E: i ∈ S (1 - r_iq) + ∑_iq ∈ E^C: i ∈ S a_iq; 
             p,q ∈ T; S,T ∈.


This aims to ensure the nodes in the same division receive equal
influence from all other divisions; this is the SymRep equivalent of
Eq. (<ref>). An optional constraint ensures that the
in-degree of all nodes is at least one.

    ∑_ip ∈ E (1 - r_ip) + ∑_ip ∈ E^C a_ip≥ 1,        p ∈ V.

This helps avoid scenarios where nodes in a fiber are disconnected
from the rest of the network. The following constraints are necessary
but not sufficient for minimally balanced colorings
Eq. (<ref>).

    ∑_ip ∈ E: i ∈ R (1 - r_ip) + ∑_ip ∈ E^C :i ∈ R
        a_ip - (∑_iq ∈ E: i ∈ R (1 - r_iq) + ∑_iq
          ∈ E^C :i ∈ R a_iq) ≥ s_pqR + ns_qpR; 
     p
        ∈ S; q ∈ T; R,S,T ∈



    ∑_iq ∈ E: i ∈ R (1 - r_iq) + ∑_iq ∈ E^C :i ∈ R
        a_iq - (∑_ip ∈ E: i ∈ R (1 - r_ip) + ∑_ip
          ∈ E^C :i ∈ R a_ip) ≥ s_qpR + ns_pqR; 
     p
        ∈ S; q ∈ T; R,S,T ∈,



    s_pqR + s_qpR≤ 1;    p ∈ S; q ∈ T; R,S,T ∈,




    ∑_R ∈ (s_pqR + s_qpR)    ≥ 1;     p ∈ S; q ∈ T;
        S,T ∈
    ∑_O ∈-(S∪
          T) (s_pqO + s_qpO) + ∑_I ∈ (S∪ T) |s_pqI -
        s_qpI|    ≥ 1;     p ∈ S; q ∈ T; S,T ∈



These inequalities define the limits to ensure that two different
divisions do not have the same influence, termed the unbalancing
constraint. The entire model is defined as follows.

    min{f_α,β: (<ref>),
    (<ref>), (<ref>),
    (<ref>),(<ref>),
    (<ref>), r_ij, a_kℓ, s_pqR∈{0,1}, ij
    ∈ E,kℓ∈ E^C, p ∈ P, q ∈ Q, P≠Q,R ∈}.

where Eq. (<ref>) within the equation above is a
reference to only select one of its
sub-equations. Eq. (<ref>) restricts feasible
solutions and potentially better minimal balanced solutions relative
to Eq. (<ref>) but is stronger with respect to the
minimal balanced coloring property. In our implementation, we always
try to find a solution to our networks by first implementing
Eq. (<ref>) after which if it fails to produce a
minimal colored solution, then Eq. (<ref>)
is used.  We solve the integer linear programs with the solver
Gurobi. The implementation is available at
<https://github.com/makselab/PseudoBalancedColoring> and
<https://osf.io/26u3g>.




§ ACKNOWLEDGMENTS

We acknowledge discussions with Ian Leifer and Luis Alvarez. We thanks Ev Yemini for the NeuroPal animals and advice. 



§ FUNDING

Funding was provided by NIBIB and NIMH through the NIH BRAIN Initiative Grant  R01 EB028157 and ONR Grant N00014-22-1-2835. The work in this manuscript was partially supported by the Simons Foundation (#543069) (M.Z.) and the International Research Scholar Program by the Wellcome Trust and Howard Hughes Medical Institute (#208565/A/17/Z) (M.Z.).



§ AUTHORS CONTRIBUTIONS

B.A. and T.G. performed the analysis of the data. P.A. performed Ca++ imaging experiments, neuron identifications, and the analyses shown in Fig. 2. D.P. performed the algorithmic developments. M.Z. led the experimental work and provided advice for the theory part. H.A.M. led the modeling and theoretical work. B.A., P.A., D.P., T.G., M.Z., and H.A.M. wrote the manuscript. M.Z. and H.A.M. led the study. 



§ DATA AND CODE AVAILABILITY

All data and code are available at <https://github.com/MakseLab>,
<https://osf.io> and <http://kcorebrain.com>.  In particular,
an implementation of the algorithm to find minimal balanced colorings
in a graph in the form of an R package is available at
<https://github.com/makselab/fibrationSymmetries> and
<https://osf.io/z793h/>. An implementation of the MILP solver is
available at <https://github.com/makselab/PseudoBalancedColoring>
and <https://osf.io/26u3g>.
    


naturemag



§ SUPPORTING INFORMATION



























