Synchronization Patterns in Modular Neuronal Networks: A Case Study of C. elegans

We investigate synchronization patterns and chimera-like states in the modular multilayer topology of the connectome of Caenorhabditis elegans. In the special case of a designed network with two layers, one with electrical intra-community links and one with chemical inter-community links, chimera-like states are known to exist. Aiming at a more biological approach based on the actual connectivity data, we consider a network consisting of two synaptic (electrical and chemical) and one extrasynaptic (wireless) layers. Analyzing the structure and properties of this layered network using Multilayer-Louvain community detection, we identify modules whose nodes are more strongly coupled with each other than with the rest of the network. Based on this topology, we study the dynamics of coupled Hindmarsh-Rose neurons. Emerging synchronization patterns are quantified using the pairwise Euclidean distances between the values of all oscillators, locally within each community and globally across the network. We find a tendency of the wireless coupling to moderate the average coherence of the system: for stronger wireless coupling, the levels of synchronization decrease both locally and globally, and chimera-like states are not favored. By introducing an alternative method to define meaningful communities based on the dynamical correlations of the nodes, we obtain a structure that is dominated by two large communities. This promotes the emergence of chimera-like states and allows to relate the dynamics of the corresponding neurons to biological neuronal functions such as motor activities.


INTRODUCTION
Synchronization phenomena are widely studied across fields, from classical mechanics [1] to complex dynamical systems [2][3][4][5] and music [6,7]. Surprising phenomena in nature, for instance, the synchronized flashing of fireflies [8] or the unexpected motion of bridges due to the emergence of synchronized walking [9] have sparked the interest in synchronization patterns. However, some more peculiar patterns of synchronization can also be observed in complex systems. These include the surprising coexistence of coherent and incoherent parts of coupled identical oscillators, a hybrid state that became known as chimera.
Previous works on the effect of nontrivial topologies on chimera states have involved scale-free and random networks [51,52], hierarchical (fractal) schemes [53], modular structures [54], and "reflecting" connectivities [55]. Our aim is to contribute in this direction and take it one step further by considering a multilayer structure. In recent years, the study of multilayer networks has become highly popular owing to their significant relevance in several complex systems [56][57][58]. In the context of neuronal networks, such a multilayer approach is ideal for addressing the relationship between structure and function, an essential question in theoretical neuroscience [59]. Concerning chimera states, studies on multilayer/multiplex networks are limited and mainly deal with artificial coupling schemes [23,[60][61][62][63]. For example, the case of two populations with various coupling schemes has been systematically studied in reference [22]. In the present work, we focus on the possibility to observe chimera-like patterns in a multiplex structure of a realworld system, namely the neuronal network of Caenorhabditis elegans (C. elegans). Our main focus is to demonstrate the existence and emergence of synchronization patterns in a multilayer network obtained from the connection of this real organism. In short, we will show that chimera-like states can be hard to identify in real-world networks and suggest an alternative approach to dynamically define communities, whose dynamics can be related to biological functions to the involved neurons. The C. elegans multilayer connectome is used as a case study, but the proposed approach can be easily generalized to other networks.
The nematode C. elegans has been studied for many decades as a standard model organism for many processes of biological interest and beyond [64]. Particularly for neurobiology, the structure and connectivity of its nervous system has been deducted from reconstructions of electron micrographs of serial sections [65,66]. Its nervous system includes sensory organs in the head and can produce highly plastic behavior, e.g., disassociative and associative learning and memory as response to taste, smell, temperature, touch and slightly to light, even though the nematode has no eyes [67]. A number of molecular mechanisms is involved in learning and memory, mediated through the same neurotransmitters as in humans and every species with a nervous system. In fact, neurons in C. elegans are very similar to those of humans, and their synapses are also classified as electrical or chemical.
It has been found that some of C. elegans' neurotransmitters, specifically the four monoamines dopamine, octopamine, serotonin and tyramine, act at both neurons and muscles to affect egg laying, pharyngeal pumping, locomotion and learning, or in general, modulate behavior in response to changing environmental cues [68]. Not only in C. elegans but also in many animals, one important route of neuromodulation is through monoamine signaling, and it is well known that this extrasynaptic communication is critical to some human brain functions. In both humans and C. elegans, many neurons expressing aminergic receptors are not post-synaptic to releasing neurons, indicating that a significant amount of monoamine signaling occurs outside the wired connectome. This defines a wireless connectivity network between neurons [69] for C. elegans. In general, this wireless or extrasynaptic communication is known as volume transmission in neuroscience [70][71][72], and is characterized by three-dimensional signal diffusion in the extracellular fluid, for distances larger than the synaptic cleft. This leads to multiple extracellular pathways connecting intercommunicating cells that are not well characterized from a structural perspective. For further details about our implementation of this network please refer to the Methods 5.1.
In our multilayer approach to modeling the neuronal connectivity of C. elegans, the network's nodes represent neurons connected by either electrical, chemical, or wireless pathways, which defines three layers. The electrical network is undirected, while the chemical and wireless networks are directed. Taking a closer look at the three layers' hubs in Figure 1, a strong overlap exists between hubs in the electrical and in the chemical network, however not so in the wireless network. This is in agreement with results from Bentley et al. [69], where a multilayer analysis of the three networks delineates topological overlaps between the two synaptic networks and discrepancies between the wireless and the synaptic networks. These topological differences of the different layers can be further explored by the degree distributions depicted in Figure 1, where it is clear that all of them present a heavy-tailed degree distribution with a maximum degree/average degree of 40/4 (electrical layer), 53/16 (chemical layer), and 137/15 (wireless layer). Since there are not many monoamine-releasing neurons, we observe many nodes with zero in-and out-degree. Note that the values of the wireless degree distribution need to be interpreted as maximum possible values. The underlying parameters that could lead to an equivalent of link weights are disregarded due to their complexity. Biologically, a source neuron in the wireless network influences a potential target neuron depending on parameters related to the diffusion processes of neurotransmitters throughout the body. Such parameters can be the physical distance between the presynaptic source and the postsynaptic target or the chemical concentration of the released neurotransmitters [70][71][72].
In order to investigate the synchronization patterns of the neuronal network at hand, we will first investigate the dynamics based on a previous approach from Hizanidis et al. [54]. This modeling approach is then extended to get closer to the biological connectome of C. elegans. The synchronization of the network is computed based on the pairwise Euclidean distances of the dynamical variables of nodes (cf. Methods 5.3 and Kemeth et al. [73]). Furthermore, we introduce an alternative way of finding meaningful communities in the neuronal network and relate the observed synchronization patterns-including chimera-like states -to biological functions of the involved neurons.

COUPLING BY DESIGN
Before investigating the C. elegans network using the actual connectivity data, we discuss results following the modeling approach of Hizanidis et al. [54], where first, the communities are computed from the aggregated connectome irrespective of the link type. Then, all intra-community links are assigned to the electrical layer and all inter-community links to the chemical layer. In other words, the designed network is modular or multilayer, where the neurons of each module and their intra-community links occupy a different layer. We follow this approach in order to test the applicability of the Euclidean distance method in evaluating the synchronization of the network. Figure 2 shows the six communities obtained from the Walktrap algorithm [74] applied to the aggregated connectome of C. elegans. The resulting topology is exactly the same as in Hizanidis et al. [54]: electrical intra-community links and chemical inter-community connections. Therefore, numerical integration of the coupled Hindmarsh-Rose (Equation 1) using this topology leads to similar time series of the dynamical variable as in Hizanidis et al. [54]. See Methods 5.2 for details on the Hindmarsh-Rose model. In Figure 3, the level of synchronization of every community γ 1 to γ 6 , the global level of synchronization γ of the whole network, the chimera-like index χ γ as well as the metastability index λ γ are shown for a large range of electrical and chemical couplings. Note that the global level of synchronization γ is computed from the pairwise Euclidean distances of the three dimensional coordinates defined by the dynamical parameters p, q and n, between any pair of nodes in the network (cf. Methods 5.3). Thus, γ , in contrast to χ γ and λ γ is independent of the community structure.
It can be seen that for an electrical coupling strength of g el = 0.4 and no chemical coupling (g ch = 0) the chimera-like index based on Euclidean distances attains a value of χ γ ≈ 0.25. Furthermore, the metastability index for these coupling strengths λ γ ≈ 0.12 is only half as large as the chimera-like index. This serves as justification why we may call the state a "chimera-like" state, in contrast to a "metastable" state where the metastability index prevails. Figure 4A depicts a space-time plot of the p-variable of the Hindmarsh-Rose model, which corresponds to the membrane potential (cf. Methods 5.2) for disconnected communities with high electrical (intra-community) coupling. One can see that every community operates in the synchronized regime. For smaller electrical coupling, Figure 4B shows the space-time plot of the p-variable with a pattern of mixed synchrony that resembles a chimera-like state, that is, varying synchronization across communities. This state is achieved for g el = 0.4 and g ch = 0, that is, that communities are not connected. It can be seen that especially the nodes in the larger communities 2 and 4 are much more synchronized than the nodes in the small communities. The chimera-like index is still large in a close neighborhood of this point (cf. Figure 3). However, when the chemical coupling is increased, the intra-community synchronization weakens and inhibits the emergence of chimeralike states. In summary, different levels of synchronization can be achieved by means of reducing the inter-community coupling strength. This raises to the hypothesis that, when observing the designed model, the main driver of chimera-like behavior is in fact the relative size of the communities, since larger communities do not need a high intra-community coupling strength to reach a high level of synchronization. Figure 4C shows the mean order parameter depending on the size of the community. As expected, the order parameter grows with the number of nodes in the community. While the mean order parameter of the three small communities (1, 3, and 5) is always below 0.6, the largest community reaches a value of γ 4 ≈ 0.95.
The actual size of the community affects the order parameter only in an indirect way. What seems to be more important is the higher mean degree of nodes in the larger communities. In Figure 4D, it can be seen that there is a strong correlation between the number of nodes and the mean node degree of each community. In other words, since the nodes in large communities have more neighbors than the nodes in small communities, it is easier for them to synchronize. Note that this does not imply causality between node degree and level of synchronization. To corroborate this, more general investigations, e.g., on randomized versions of the network would be insightful. For our purpose however, it suffices to know that different community sizes influence the level of synchronization significantly.

COUPLING BY BIOLOGY
In order to investigate the synchronization of the C. elegans network even further, a third layer is added to the graph, representing the monoamine connections (wireless network). Furthermore, the assumptions about electrical and chemical synapses made in Hizanidis et al. [54] related to intra-and inter-connectivity are dropped, and the three-layer neuronal network is created using the actual connectivity data of the three synapse types (see Methods 5.1). In this section, we present two approaches to finding communities in this biological multilayer network.

Communities Based on Topology
The communities are first evaluated using a Multilayer-Louvain community detection algorithm (see Methods 5.4), which yields 8 communities instead of 6 as in the aggregated case discussed in the previous section (cf.  wireless network, presenting few intra-community links, is distributed almost randomly across the network, the chemical layer for instance strongly reflects the underlying community structure in the adjacency matrix. Since chemical connections make up most of the edges in the network, the algorithm optimized the community partition mainly based on the chemical layer. Comparing this partition to the one in Figure 2, one can observe that the clear separation between edge types into intra-and inter-community edges is not achieved when using the biological connectome without assumptions. The unclear separation of edge types in the partition shows its effects when analyzing the dynamics of the system using the Hindmarsh-Rose equations (see Methods 5.2). Figure 6 shows the parameter scans for the global level of synchronization γ , the chimera-like index χ γ and the metastability index λ γ for two different wireless coupling strengths. For the level of synchronization γ 1 to γ 8 within the individual communities (see Figures S1, S2). First of all, the global level of synchronization of the system is highly reduced when using the real connectome, which can be explained by the previously mentioned edge distribution. Since the electrical layer synchronized the nodes within communities in the designed partition, the intercommunity coupling could easily be tuned using the chemical coupling strength. Therefore, a clear chimera-like region could be observed in the parameter scans in Figure 3. In the case of the real connectome, it is not possible to tune intra-and intercommunity coupling separately, since all edge types are distributed across and within communities.
Even though the parameters used to identify chimera-like states (γ , χ γ , and λ γ ) are significantly lower than in the model by design, the timeseries of the neuron membrane potential p in Figure 7 suggest that the system adopts three different synchronization patterns, depending on the dynamical coupling strengths g wl , g el and g ch . The corresponding coupling strengths and synchronization parameters are noted in Table 1. Figures 7A-C show the evolution of p for g wl = 0.0, and Figures 7D-F show p for g wl = 0.2. One observes that for high electrical coupling strengths (Figures 7A,D), the system is synchronized, meaning that most of the nodes follow the same periods of bursting (green) and quiescent (blue) states. However, one can see singular nodes that fall out of the synchronized pattern (especially in communities 1 to 3), which correspond to the nodes that are not electrically coupled. Figures 7C,F depict the system in the desynchronized regime, where one cannot distinguish any clear pattern of synchronization. The desynchronization becomes higher when the wireless coupling strength is increased. In Figure 7B, the system seems to exhibit chimera-like behavior at first glance. However, the chimera-like index for this coupling parameter set is very small (see Table 1). Even though community 3 seems less synchronized, one can read from the synchronization parameters (cf. Figure S1) that in fact all communities present very low synchronization. This behavior does not change significantly when adding wireless coupling (cf. Figure 7E).
Using the Multilayer-Louvain community detection approach to partition the network, one observes a system expressing different synchronization patterns depending on the interplay of the coupling strengths of the three layers. However, even though these patterns are visible in the evolution of the p-variable (Figure 7), the topological community partition cannot reflect these patterns in the parameter scans ( Figure 6). Therefore, the question about the meaning of community detection can be raised. The following section proposes an alternative way of detecting communities in multilayer networks, which is not based on the topology or link distribution, but on the correlation of nodes with respect to their dynamic variable p.

Communities Based on Dynamics
In the previous section it was shown that the Multilayer-Louvain partition already leads to significant differences in the synchronization behavior of the distinct communities. Yet, as this barely becomes apparent in the chimera-like index, we investigated an alternative approach of finding communities, based on dynamical correlations between the time series of the p-variable. The heuristic algorithm leading to the correlationbased partition is described in detail in Methods 5.5. It is worth mentioning that this algorithm only changes the partition of the network, it does not change the nodes' dynamics when compared to section 3.1, as the underlying graph itself remains untouched. However, the synchronization measures (i.e. γ m , χ γ , and λ γ , see Methods 5.3) do intrinsically depend on the particular partition. Figure 8 shows that the correlation-based partition strongly differs from the partition found with the Multilayer-Louvain algorithm. The most striking difference is the stronger heterogeneity among inter-community links in the correlationbased partition. The outgoing and incoming links between communities are relatively evenly distributed in the topologybased Multilayer-Louvain partition. In the correlation-based partition however, the two largest communities (3 and 6) are much stronger connected than the rest of the network. Furthermore, it can be seen that the two smallest communities (5 and 8) present no electrical connection to any other community in the network. In general, the electrical sub-network shows only very few inter-community links compared to the strong intra-community coupling (self-loops) of the two largest communities (3 and 6).
Regarding the dynamical properties, the dynamical correlation-based partition leads to qualitatively similar results as the Multilayer-Louvain partition. Compare Figures 6  and 9. In particular, Figure 9B clearly shows that the highest values of the chimera-like index χ γ are still obtained for high electrical couplings and small chemical couplings. Increasing the wireless coupling as in Figure 9E reduces the value of χ γ , similar to what has been observed previously in the dynamical analysis of the Multilayer-Louvain partition. However, the value of the highest chimera-like index (χ γ ≈ 0.14, obtained at g el = 1.96 and g ch = 0.04) is significantly higher for the correlation-based partition. Moreover, it is also higher than the corresponding metastability index (λ γ ≈ 0.07). Therefore, we may indeed call the state "chimera-like" [38,54,[75][76][77][78].
The reason for the high chimera-like index can be observed in the space time plots of the p-variables in Figure 10. As was mentioned before, the dynamics of single oscillators do not depend on the partition. In order to compare the results from the correlation-based partition with the Multilayer-Louvain partition, we decided to show the time series of p using the same coupling constants as in Figure 7. The different ordering of the nodes, though, leads to a significantly higher homogeneity between nodes within one community, which is especially apparent for the largest communities (3 and 6) in Figures 10A,D. As a consequence, the respective levels of synchronization (γ 3 ≈ 0.40 and γ 6 ≈ 0.67) are higher than for the other communities, which raises the chimera-like index. The reason for this high level of synchronization is the strong intra-community coupling of the two large communities in the electrical layer, since the electrical coupling strength is the very high in parameter set (g el = 1.80).
For time series with small electrical coupling (see Figures 10C,F), the communities are not separated as clearly, particularly when wireless coupling is applied to the system, which can be observed in Figure 10F. Especially the largest community (3) is much less synchronized, which can be explained by a stronger influence of the chemical and wireless layers, where a high number of intra-community links exist. This demonstrates that the dynamical correlation-based algorithm seems to preferentially sort the network according to the electrical sub-network, as this layer seems to be most important for the overall synchronization of the network. The adjacency matrices in Figures 8A-C support this observation: Only the electrical layer shows a large number of intra-community links, whereas the other layers present no clear visible structure. This is again in contrast to the Multilayer-Louvain partition, where the chemical layer also shows a pronounced community structure (see Figure 5B). In the case of intermediate electrical, chemical coupling and no wireless coupling (see Figure 10B), the distinct communities can still be identified, yet the level of synchronization in the large communities does not suffice to reach a high chimera-like index. Adding wireless coupling as shown in Figure 10E does not lead to significant changes in the values.
For a full review of the different coupling strengths of to the system that lead to the time series in Figure 10, as well as the consequent values of the synchronization parameters γ , χ γ , and λ γ , please refer to Table 2.

DISCUSSION
Interesting synchronization patterns were found using different modeling approaches in the multilayer network of C. elegans.  Table 1. Red lines separate different communities.
They were quantified based on the pairwise Euclidean distance between the dynamical variables p, q and n of the underlying Hindmarsh-Rose system (see Methods 5.3).
Following the approach of Hizanidis et al. [54], we first assume purely electrical connections within the communities, and chemical synaptic intercommunity coupling (section 2). This results in a designed coupling scenario, where chimera-like states are clearly visible due to a strong separation of connection types. While the electrical sub-network synchronizes the communities within themselves, the chemical sub-network only allows for connections across communities.
Moving toward a more biologically-inspired modeling (section 3), these synchronization states are more difficult to observe. Since the edge types are not clearly separated anymore and it is therefore impossible to tune intra-and intercommunity coupling separately, the nodes within one community cannot synchronize as easily. This is especially the case when partitioning the network with the Multilayer-Louvain algorithm: the synchronization patterns are visible in the time series (see Figure 7), but the synchronization indices are very low (see Figure 6).
We also discussed an alternative way to identify correlated clusters in the network, namely to sort nodes in communities according to the Pearson correlation matrix of the p-variable (see Methods 5.5). In this case, the community structure is dominated by two large communities with a high amount of electrical self-loops (see Figure 8) that present a strong synchronization (see Figure 10). There are two small communities that share no electrical links to the rest of the network and therefore scarcely synchronize with nodes from the other communities. This promotes the emergence of chimera-like states. Further insights on the dynamical correlation-based partition can be obtained by investigating the neuronal functions of the nodes. In Table 3, one can see that the highly synchronized communities (3 and 6) contain 75% of the motor neurons of the system. The synchronization of motor neurons under varied coupling strengths is in harmony with experimental findings, for example in rats [79].
In this context, a question could be raised regarding the multilayer nature of the studied network. Since the three layers do not share the same number of neurons (only 253 of the 279 neurons are connected by electrical gap junctions), a certain  Figure 3, except for g wl . Refer to Figures S3, S4 for the parameter scans including γ i . group of nodes are prone to remain desynchronized for certain combinations of g el , g ch and g wl . However, the strong biological interplay between synapse types is crucial to the understanding of the neuronal network as an entity [80,81]. Therefore, the connectome should be modeled as a multilayer network.
Keep in mind that the studied three-layer network contains information only about the electrical, chemical and monoamine connections. Another layer could be added for the neuropeptide wireless network, which was not included since many neuropeptide receptors, as well as ligands for many neuropeptide receptors are unknown. Also, the distance over which neuropeptide signaling can occur is uncharacterized for many of them.
Furthermore, concerning the synchronization metric based on Euclidean distances (see section 5.3), the threshold parameter which defines the limit between synchronized and desynchronized nodes has been set to δ = 0.01 as in reference [73]. This value could be adapted to better suit the system and the three-dimensional distances.
This work presents an approach for analyzing the complex biological network of C. elegans using metrics of synchronization based on Euclidean distances and a new method of finding clustered nodes by correlating their dynamical variables. The underlying framework can be extended for multiple complex network applications.

Datasets
The gap junction and chemical synapse networks of a hermaphrodite C. elegans have been obtained in [66], and are available through WormAtlas [82]. The associated adjacency matrices are computed for the electrical and chemical layers, where we omitted the neuromuscular junctions since we are only interested in neuronal interaction. Note that this dataset does not include the 20 pharyngeal neurons. Hence, we work with the somatic giant component of the neuronal network.
The electrical sub-network consists of 253 neurons and 890 synapses or gap junctions from 517 unique neuron pairs (including 3 self-connections). A total of 352 out of 517 neuron pairs have only one synapse between them, while the other 165 pairs show multiple parallel connections, with a maximum value of 23. This means that the respective symmetric (undirected) adjacency matrix has weights varying from 1 to 23 for connected neurons.
The chemical sub-network contains 253 source and 268 target neurons, the union of both sets is composed of 279 neurons, which is the total number of nodes of the modeled C. elegans network. There are 6,294 synapses from 2,575 unique sourcetarget neuron pairs. A total of 1,362 out of the 2,575 neuron pairs have only one synapse, while the other 1,211 pairs have multiple  Table 2. Red lines separate different communities. The time series are identical to Figure 7, but for a different ordering of nodes.
synaptic connections with a maximum value of 37. Therefore, the associated asymmetric (directed) adjacency matrix has weights varying from 1 to 37 for connected neurons.
The wireless sub-network of this study is restricted to the monoamine network in Bentley et al. [69] and is available in Bentley et al. [83]. Again, the pharyngeal neurons are excluded. This network by itself can be thought of as a directed quadripartite network composed of a source neuron, a neurotransmitter, a receptor and a target neuron. The considered wireless network is composed of 16 source neurons, 4 monoamine neurotransmitters, 16 associated protein receptors and 215 target neurons. As a first approach to the implementation of the wireless connectivity within our dynamical model, we are only interested in an effective connectivity between a source and  Figure 10.  a target neuron. For this purpose, we reduce the quadripartite nature by assigning binary weights to the adjacency matrix: 1, if there is any path between the neurons and 0, if there is no possible connection from a source to a target neuron through any neurotransmitter and matching receptor. The final directed adjacency matrix includes 2,282 edges. The functional classification of the neurons in three categories (sensory, motor and interneurons) has also been obtained from [66] and manually created based on the dataset in WormAtlas [84]. We excluded the male data, since we only consider the hermaphrodite information. If the description of a particular neuron includes both interneuron and motor characteristics, we choose to describe it as motor neuron, as well as we define amphid neurons to be sensory, following the approach in Varshney et al. [66].
For details relevant to our study on the individual neurons and their characteristics, please refer to the Table S1.

Hindmarsh-Rose Dynamics
We consider a network of neurons locally characterized by Hindmarsh-Rose dynamics [85], a model intended to describe the transition between the stable resting state and the stable limit cycle of neurons. The model was intended for two types of coupling (electrical and chemical), we propose extending it by a third coupling term to account for the wireless connections in the studied network, as described by the following equations: where i = 1, . . . , N is the neuron index, p i is the membrane potential of the i-th neuron, q i is associated with the fast current, either Na + or K + , and n i with the slow current, for example Ca 2+ . The parameters of Equation (1) are chosen such that a = 1, b = 3, c = 1, d = 5, s = 4, p 0 = −1.6, and I ext = 3.25, for which the system exhibits a multi-scale chaotic behavior characterized as spike bursting [86]. The parameter r modulates the slow dynamics of the system and is set to 0.005 so that each neuron lies in the chaotic regime in the absence of coupling [54]. For these parameters, the Hindmarsh-Rose model enables the spike-bursting behavior of the membrane potential observed in experiments made with single neurons in vitro. We choose random initial conditions in the time series shown and find that changing the initial conditions does not change the longterm synchronization behavior. In our time series analysis, we always remove the transients and average over a long time period. The chaotic behavior of the Hindmarsh-Rose model has been studied in earlier publications with slightly different parameters, which led to the investigation of a plethora of chaotic phenomena using spike-counting techniques [87] and detailed bifurcation analysis [88]. The connectivity structure of the electrical synapses is described in terms of the Laplacian matrix L, whose elements are defined as L ij = A el ij − δ ij k el i , where k el i is the degree of node i in the electrical layer and δ ij = 1 if i = j, δ ij = 0 otherwise. A el is the symmetric adjacency matrix whose elements are A el ij = 0 if there are electrical synapses connecting the neurons i and j, and A el ij = 0 otherwise. The strength of the electrical coupling is given by the parameter g el and its functionality is governed by the linear function H(p) = p.
The connectivity structure of the chemical synapses is described by the adjacency matrix A ch , whose elements are A ch ij = 0 if there are chemical synapses between neurons i and j, and A ch ij = 0 otherwise. The nonlinear chemical coupling is described by the sigmoidal function S(p) = {1 + exp[−λ syn (p − θ syn )]} −1 , which acts as a continuous mechanism for the activation and deactivation of the chemical synapses. The associated coupling strength is noted as g ch . For the chosen set of parameters, |p i | < 2 and thus (p i − V syn ) is always negative. Therefore, the chemical coupling is excitatory if V syn = 2. The other parameters are θ syn = −0.25 and λ syn = 10, following references [89,90].
The wireless connectivity structure is described by the adjacency matrix A wl . It is also considered nonlinear; however, much slower than the chemical synaptic coupling [69]. Intuitively, the exponential function S(p) in the denominator can be decreased by replacing λ syn withλ syn ≪ λ syn . We choseλ syn = 1, V syn = 2 and θ syn = −0.25, as for the chemical coupling. Furthermore, the wireless coupling is considered as an additional disrupting signal to the synchronization of the network. It is therefore treated like excitatory chemical synapses.

Level of Synchronization
We adapt an approach based on Kemeth et al. [73] in order to compute a level of synchronization of the studied network. Instead of considering the local curvature, which is optimized for a ring network, we calculate the pairwise Euclidean distances between the variables {x 1 , x 2 , . . . , x N } at every time step t, with x i = (p i , q i , n i ). For all t, we obtain a set of all possible distances between a set of N nodes: Two nodes i and j are now defined to be synchronized if ||x i (t) − x j (t)|| ≤ δ and desynchronized if ||x i (t) − x j (t)|| > δ, where δ = 0.01 · D max is a threshold value. The value D max is the maximum possible Euclidean distance between a pair of nodes: where x max = (p max , q max , n max ) and x min = (p min , q min , n min ) are vectors containing the maximum and minimum values of the dynamical variables for all time steps t and nodes N. Therefore, two nodes are defined to be synchronized at time t if their Euclidean distance does not exceed 1% of the maximum possible distance, which is well defined for every space-time series.
Based on the set of Euclidean distances, we can measure the amount of spatially coherent nodes at each time step t. For this purpose, we consider a different set of distances, containing only those that are smaller than the threshold value δ: The fraction between the number of distances within the range of the threshold value and the possible number of distances then results in the amount of synchronized node pairs. Note that the number of node pairs grows at a rate of N 2 . It is therefore necessary to take the square root of this value, in order to make it comparable across network sizes. We call the resulting value "level of synchronization:" If γ (t) is only computed for a certain community m, it is called γ m (t), representing the level of synchronization of community m at time t. For a network consisting of M communities, it is now possible to compute the chimera-like index: at time t as proposed in Shanahan [91] and also used in Hizanidis et al. [54], where γ m (t) M denotes the average level of synchronization at time t over all communities m. Thus, the only difference to Shanahan [91] is the application of the Euclidean-distance-based level of synchronization γ (t) instead of the Kuramoto order parameter. The temporal mean then defines the time-averaged chimera-like index of the network: Similarly we can compute the metastability index: (8) of community m, where γ m (t) T denotes the temporal mean of γ m (t) over all time steps. The average over all communities yields the metastability index of the whole network: The subscript γ is utilized to emphasize that these parameters differ from the parameters in Shanahan [91] and Hizanidis et al. [54] in the way the underlying level of synchronization is computed. In order to compare community partitions with different numbers of communities, it is important to know that the ranges of χ γ and λ γ depend on the number of communities M of the studied network. The chimera-like index of a chimera state, where half of the communities is completely synchronized and the other half is desynchronized, becomes: since the deviation from the mean is 0.5 for every community. The same considerations lead to a maximum metastability index of:λ which is approximately 0.25, since the total number of time steps T is large. Hence, we obtain the chimera-like and metastability indices normalized to unity: and,

Multilayer-Louvain Community Detection
The communities discussed in section 3 are computed based on a multiplex Louvain community algorithm [92]. In single-layered graphs, the key to finding communities usually lies in optimizing the modularity function Q [93]: where A ij is the graph's N×N adjacency matrix, m l = 1 2 N i,j=1 A ij is the total link weight in the network, and k i = N j=1 A ij is the weight incident to node i. The weight of a link between i and j corresponds to the number of connections between these two nodes (see Methods 5.1). δ g i , g j = 1 if nodes i and j are in the same community, and δ g i , g j = 0 otherwise. Therefore, the term A ij − k i k j 2m l quantifies how strongly the two nodes will be coupled in the studied network, compared to how strongly they would be coupled in a random network. In the algorithm, the function Q is computed for every pair of nodes iteratively until it reaches a maximum value.
For multilayer networks, the modularity as defined in Equation (14) is not well suited as it does not differentiate if nodes are connected by different layers. In order to extend the modularity to multilayer applications, consider a network with S layers. We define the degree of node j within the same layer as ij denoting the adjacency matrix in layer σ . The generalized modularity function Q multilayer for a multilayer network with S layers is defined as [92]: where m l,σ = 1 2 N j=1 k (σ ) j is the total link weight within layer σ , and µ = S σ =1 m l,σ is used for normalization similar to m l in Equation (14). Note that in the case of the considered C. elegans network, there are no inter-layer connections, since every connection type (electrical, chemical, and wireless) represents an independent sub-network. In this study, we consider 3 layers with σ ∈ {el, ch, wl} that are shown by black, red and green link color in Figures 11A,B for the multilayer connectome and correlation-based matrix, respectively. However, Equation (15) can be extended to be used for the study of multiplex networks, in which inter-layer connections exist [92].

Dynamical Correlation Community Detection
We present a heuristic approach to finding meaningful communities based on the dynamics of the system. While previous approaches aimed to find a community structure based on the topology, we propose an algorithm which partitions the network based on the nodes' correlations of the p-variable. Figure 12 shows a schematic description of the algorithm. In order to gain insight on the synchronization of the time series for each pair of nodes, we compute the Pearson-correlation matrix from the p-time series of all nodes In the hereby created matrix, every entry represents the correlation value of two time series of the two respective nodes (cf. the matrix in Figure 12, top left corner).
In order to find a community partition in the correlation matrix, we employ the stochastic block model approach from the graph_tool framework [94]. Since this framework does not intrinsically support weights in the network, we created a graph with a discrete number of edges between two distinct nodes that depends on the link weight in the correlation matrix; higher link weight therefore corresponds to a larger number of edges.
For every parameter set (i.e. every combination of the three different coupling strengths g el , g ch and g wl ) we obtain one correlation matrix, on which we apply the graph_tool algorithm. Note that, for some parameter sets, the algorithm will not find a "reasonable" partition. In particular, some partitions may consist of very small communities with only very few nodes. This is problematic in terms of the level of synchronization γ m , since one node in a small community m plays a bigger role in the synchronization of this community; this implies stronger fluctuations of γ m . Therefore, very small communities FIGURE 12 | Schematic description of the dynamical correlation-based community detection algorithm. In the sorted matrices red lines represent the borders between two communities.
(especially communities consisting of only two or three nodes) can have a significantly stronger influence on the chimera-like and the metastability index than large communities. This is why we only consider community partitions with at least six nodes per community.
Another constraint applied to the partitions is a lower boundary for the level of synchronization in at least one community. The constraint is needed, because nodes do not synchronize as easily in the system based on the connectivity data (see section 3). However, a highly synchronized community is crucial to finding chimera-like states. The threshold value used to filter out partitions containing low-synchronized communities was set to γ thr = 0.30. This is a reasonable compromise between reaching a high level of synchronization in at least one community and still keeping a relatively high number of partitions.
The algorithm finds 582 partitions that satisfy the constraints out of an initial set of 50 · 15 · 7 = 5250 possible partitions (g el ∈ [0.04, 0.08, ..., 2.00], g ch ∈ [0.02, 0.04, ..., 0.30] and g wl ∈ [0.00, 0.05, ..., 0.30]). Subsequently, we iterate over all pairs of nodes (i, j) and count how often they are found in the same community for the 582 partitions. In other words, if a pair of nodes always finds itself in the same community, the counter will be 582, while a pair that is always found in distinct communities will receive a counter of 0. This then leads to a 2dimensional histogram as can be seen in Figure 12 in the bottom right corner.
As a final step, this histogram is sorted using the graph_tool algorithm in order to find a merged community partition that contains information over a large set of parameter values g el , g ch , and g wl . The resulting sorted histogram is shown in Figure 12 in the bottom left corner and the network is visualized in Figure 11B. There is one community consisting of nodes that almost always find themselves in the same community. This implies that the time series of the corresponding nodes have a high correlation value for almost every parameter set. All the results from section 3.2 were created using this community partition.
Please note that the proposed algorithm is only one possible way of finding partitions based on a system's dynamical behavior. In fact, it invites to further explore the interplay between topology and dynamics.

AUTHOR CONTRIBUTIONS
AP, LM, and JR had the lead analyzing the network data and performing the numerical simulations. All authors designed the study, analyzed the results, and wrote the manuscript.