Resilience of Nematode Connectomes Based on Network Dimension-reduced Method

The whole map of nematode connectomes provides important structural data for exploring the behavioral mechanism of nematodes, but to further reveal the functional importance and resilience pattern of nematode neurons, it is necessary to effectively couple the regulatory relationship between neurons and their topology. Here, with a typical signal excitation function we propose a model to capture the interacting relationship between the neurons, because a differential equation depicts the activity of a neuron, n neurons mean we need high-D differential equations to capture the neural network. With mean-field theory, we decouple this N-dimension question into a one-dimension problem mathematically. In our framework, we emphatically analyze the characteristics, similarities and differences of the structure and dynamical behaviors of the neuronal system for Caenorhabditis elegans and Pristionchus pacificus. The comparing results of simulating method and theoretical approach show that the most important homologous neurons between C.elegans and P.pacificus are I2 and NSM, which may lead to their different behavior characteristics of predation and prey. At the same time, we expect that the x eff index can be used to reveal the importance of neurons for the functional evolution and degeneration of neural networks from a dynamic perspective. In the hermaphroditic and male C.elegans, we test the control level of the intermediate neuron groups over the output neuron groups and the single neuron. These results suggest that our theoretical approach can be used to reveal the effects of bio-connectivity groups, potentially enabling us to explore the interaction relationship of neural networks in humans and animals.


INTRODUCTION
Due to special physiological structure and easy modeling, nematode has become the primary model to reveal the neurons structure and functional mechanism for humans and animals [1,2]. After four decades of exploration, nematodes are the first organisms, for which the wiring diagram of their entire nervous system has been mapped at the cellular level [3,4]. Even Witvliet et al reconstructed full brain of eight isogenic Caenorhabditis elegans individuals across postnatal stages to investigate how it changes with age [5]. This contributes to a more detailed and accurate understanding and exploration of human and animal behavior at the neuronal level, as well as an analysis of the functional importance of individual neurons at the systemic level. Therefore, we coupled the structural information of the neuronal map of the nematode and the functional relationship between the neurons, predicted the functional involvement of specific neurons or synapses in defined behavioral responses [6,7].
Resilience is system's ability to retain its basic functionality while errors, failures and environmental changes occur, which universally presents in most dynamical systems [8]. Thanks to the rapid development of network theory, the description of the nonlinear dynamics that governs the interactions between the neurons have been extracted. According to external disturbance for neuron systems, deleting interneurons or synapses may disconnect the networks, or make a great impact on the individual behavior [9]. The methods to explore the nematode's resilience mainly can be divided into two categories. 1) Biological experiments: such as the gene labeling and the neuron ablation. They can track some specific synapses or their function in detail. The labeling approach can capture the directionality of synaptic connections, and its quantitative analysis of synapse patterns display excellent concordance with electronic micrograph reconstructions [10,11]. However, this approach takes a lot of energy and is not universal. 2) Theoretical experiments: using different mathematical methods and engineering mechanisms, which can abstractly study the function of outstanding nodes in nematodes [12][13][14][15][16]. For example, according to the theory of symmetry group, they found that the symmetry of the neural network has directly biological significance, and its correctness can be strictly proved by using the mathematical form of symmetry group. This form makes it possible to understand the importance of the structure-function relationship [17]. In addition, Yan also applied network control principles to the connectomes, which reveals both neurons with known importance and neurons which was previously unknown [18]. Within a network control principle, they add input signal to control the output of neurons by a linear framework. The results remain reliable with a small amount of disturbance to the reference connectors, but large disturbances are likely to cause distortion.
Here we construct a framework to explore the resilience of nematode connectomes: we use the signal excitation model to describe and solve quantitatively structure-function relationship of the nematode, which combines the biological behavior with dynamic behavior. The application of signal excitation mode is key to revealing how different behaviors, which contributes to understand the mapping from network structure to function [19]. Since the network contains a large number of neurons, each activity of synapse node or neuron node can be seen as a solution of the high-D differential equations. However, calculating the stable state of the neural system could be difficult or almost impossible, especially when the system is large-scale. With a network dimension-reduction method, we can derive an effective one-dimensional dynamic model which captures the system function of the neural networks. Finally, we use a weighted average activity x eff to measure the neural performance of the whole nematodes. Meanwhile, we introduce a weighted average connectivity β eff to express the structural strength of the nematode neuron system.
With the application of our framework into the neuron networks of P.pacificus and C.elegans, the differences between the predator and prey behaviors of these two nematodes can be clearly quantified at the neural network level, and the key neurons leading to the differences in nematode behaviors can be easily identified as well. Similarly, for hermaphroditic and male C.elegans, we apply our approach to discuss the core control of the intermediate neuron groups to output neuron groups. Core control means that the interneurons have great influence over the output neurons or all the neurons after neuronal perturbation.
Actually, our approach is a general mathematical framework, which couples the complex structure of the system with the nonlinear activation function, and enriches the research methods for exploring the functional behaviors of nematode. We can reveal the underlying principle of the neural network in this way. In the future, our theoretical methods can be applied to more specific linkage groups, to do targeted quantitative research [20].

Neural Dynamics Model
There are nonlinear regulation relationships between each synapse in neural networks. Taking hyperbolic tangent excitation function as an example in this paper, we give a neural network decoupling method with a general network structure. For the stimulation relationships between the trillions of neuron synapses in human brain or animals, it would enable the individuals show a much more complex nonlinear phenomenon of affine transformation. The activation function can be used to describe the coupling mechanism among synapses as Where, x i describes the activity of the neuron i, I is the basal activity of the neurons, R is the inverse of the death rate, α is the firing threshold, and J is the maximal interaction strength between pair of neurons. The adjacency matrix A depicts the topology of the neuron networks. The coefficient n governs the steepness of the sigmoid function, analogously to the Hill coefficient n in gene regulatory networks [18]. Actually, the Eq. 1 is made up of two parts: I − xi R describes the growth rate of the neuron itself, J 2 N j 1 A ij (1 + tanh(n(x j − α))) describes the positive excitation effect of its neighbor neurons on neuron i.

Network Dimension-Reduced Method
When the system of Eq. 1 tends to a relatively stable state, which means that the activities of the N individuals in the network are constant values, so the N-dimensional nonlinear rate equations all equal to 0. However, calculating the stable state of the system requires to solve the N-dimensional nonlinear rate equations numerically or analytically, which could be very difficult, or almost impossible especially when the system is large-scale. We present a new dimensionality reduction to explore the interactive relationship between neuronal connectomes for Frontiers in Physics | www.frontiersin.org October 2021 | Volume 9 | Article 731941 nematodes with mean-field approximation. This approximation is exact in the limit where the node activities are uniformly distributed.

Theoretical Framework for Network Resilience
Once we get the value of x i for the neuron i from the solutions of Eq. 1, we can quantify the performance and the connectivity strength of the neural networks with a weighted average metrics on the network topology following the way in [21]. Defining an operator for the degree sequence of A as Here, y i is a variable of neuron i, and y is the weighted average of all the neurons. s out j is the outdegree of neuron j. In contrast, s in i represents the indegree of neuron i. This operator is used to calculate the weighted average y over all nearest neighbor nodes of network A.
In order to simplify the N-dimension of Eq. 1 and approximate it to 1D formalism, we use a mean field hypothesis. Here we define For the states of all nodes in the system, Eq. 3 can be rewritten by the Hadamard product as Where we used the Hadamard product • , namely In order to reduce the dimensionality of this set of equations, we apply Γ(x) 1 T Ax In the process of dimension reduction, we introduce two variables to describe the state of the system, where x eff measures the neural performance of the whole nematodes, β eff expresses the structural strength of the nematode neuron system, Lastly, we have successfully decouple this N-dimension question into a one-dimension problem mathematically. Bring Eqs 7, 8 into Eq. 6, one gets Subsequently, when the system comes to a stable state, it is obvious that dx eff dt 0. Hence, we can further get that Equation 10 describes the theoretical relationship between the topology of the neuron networks and its dynamical performance. Actually, by Eq. 10 we have decoupled the complex system of Eq. 1 into a 1D problem.

Simulation Method
From Eq. 10, we can obtain the theoretical curves of x eff for β eff . To verify the correctness of our dimension-reduced method for the system resilience of x eff , we use the ODE45 function to solve the neural dynamics of Eq. 1, which can ensure we get the numerical results of x i . Then, with the steady state vector x and Eq. 2, which can separate the system from the topology parameter β eff and the behavior parameter x eff . The steps to simulate the neural dynamics are following 1. use ODE45 function to solve Eq. 1. We can obtain the stable vector x. In our study we set I n 2, R α J 1. The setting of this parameter we refer to [16]. Assuming the network has N neurons, the solution of the steady state vector is It is easy to know that in the steady state vector, the value of column i represents the steady state activity x i of neuron i.
2. compute the x eff and β eff from Eq. 9. For the vector x obtained in step 1 and the network indegree sequence s in , we can get the system topology parameter β eff and the behavior parameter x eff with the operator Γ(y).
3. compare the simulation results with the theoretical solutions. From the theoretical curve of Eq. 10, we test whether the simulation values of β eff and x eff fall on the theoretical curves. In our experiments, the simulation points agree well with the theoretical curve, which shows the feasibility of our approach.

Datasets for Neural Networks
In this study, we used two data sets: pharyngeal nervous systems of C.elegans and P.pacifica [22], all neuronal connections of hermaphrodite and male C.elegans [3]. The groundbreaking work of Cook and Bumbarger et al. provides a wealth data of neural networks, which can help explore the wiring problem from the perspective of structure. We want to couple the information of structure and dynamical function, so as to probe into the complex wiring problem of nematode from the network function and system resilience.
The connection matrix for pharyngeal nervous systems of C.elegans and P.pacifica was derived from the data set used by Bumbarger in his earlier work [22]. The matrix contains neuron connections and weights. The data and its experiment are used for homology comparison of two nematode in our study.
In fact, because of the complex wiring of the nematodes, their neurons have a precise functional grouping: input neuron, intermediate neuron, output neuron. So, for the dataset of Connectome adjacency matrices of hermaphrodite and male C.elegans are gap junction, which were provided by Steven J. Cook [3]. Notably, they presented the complete wiring diagram of an animal nervous system for the first time in 2019, including adult C.elegans of both sexes. It is an important milestone in the field of neuron science. The data corresponding to one connectome has the following information: the node, the weight, neuron names in each group. e.g., InterNeurons contains all the neurons that belong to the "InterNeurons" group, neuron names in each subgroup, e.g., InterNeurons_1 contains all the neurons that belong to the "InterNeurons_1" subgroup of the "InterNeurons" group.

RESULTS
3.1 Resilience of C.elegans and P.pacificus P.pacificus ( Figure 1A) and C.elegans ( Figure 1B) have highly similar synaptic structures, but quite different functions in the pharynx [23], we mainly explore the different connective structure and function in pharyngeal neurons.
We measured the functional differences of the connective structure for C.elegans and P.pacificus with x i . As shown in Figure 1C, by exploring the impact of the 14 pharyngeal homologous neurons, we found that the values of M1 for P.pacificus and C.elegans show the greatest difference, which means that M1 is the most functionally different synapse between the two nematodes. M1 is a motor neuron, and their structure centralities of M1 in C.elegans and P.pacificus are quite similar [22]. Our framework can help identify the potential critical neurons, which may be inconsistent with the traditional structural comparison.
Then, with the framework we explored the differences of the predation behavior for the two nematodes through the system resilience index x eff . As shown in Figure 1D, first of all, the activities of C.elegans and P.pacificus were in good agreement with our predicted theoretical curves. Secondly, we were awared that the x eff of the P.pacificus simulation point is much higher than C.elegans. That is to say, the total activity of the P.pacificus is much higher than C.elegans. Nematode feeding strategies differ greatly [24,25]. P.pacificus have a necromenic association with scarab beetles [22]. Their dauer larvae rest on the insect and resume development after the beetle's death to feed on microbes on the decaying carcass. P.pacificus can be easily cultured on bacteria [26], however, which is also predatory on other nematodes. In fact, the difference of predation strategy leads to the difference of neuron function. For example, differential but coordinated regulation of pm1 and pm3 in P.pacificus represents motor output that is specific to predatory feeding and does not exist in C.elegans [27,28]. Actually, the key characteristic of P.pacificus is that it takes predation behavior actively for a long time [29,30], while the C.elegans, which is preyed on, just eats the food around it. Thus, from the neural level, we can conclude that the system resilience index x eff could be used to reveal the essential differences of the predation behavior for the two nematodes that C.elegans are prey, and P.pacificus is more likely to be predators [31]. Thirdly, we tested 14 pharyngeal neurons, which showed that NSM and I2 are the neuron with the greatest reduction compared with the original simulation point in Figure 1D. It meant that NSM and I2 are the two most important homologous neurons. I2 is an interneuron, which play the role of coordinating corpus and tooth contractions during predation. I2 are more highly connected in P.pacificus and may function as network hubs. Although there are no more detailed studies showing the important function of NSM neurons, the high impact of which on the system resilience in our results is sufficient for further attention of biological experiments. Compared with Bumbarger's earlier biological experiments [22], the difference is that they concluded that I1 and I2 are the most important candidate neurons from the view of the system structure, but we found NSM and I2 are the two most important homologous neurons while considering the synaptic structural and functional characteristics simultaneously.
In addition, a phenomenon occurs in our experiment from Figure 1D, after neuronal perturbation (node deleting) of a few homozygous neurons, the x eff of P.pacificus decreased, but x eff of C.elegans increased. M4 is one of these special neurons. With its necromenic beetle association, P.pacificus has an intermediate position between the microbivorous C.elegans and true parasites [32,33]. In other words, C.elegans are prey and P.pacificus is more likely to be predators. As one of the important feeding neurons in the pharynx, M4 is closely related to predation strategy. So we made a bold guess that M4 may degenerate in P.pacificus, but evolve in C.elegans. More details in the Discussion. As shown in Figure 2A, in hermaphrodite C.elegans, no matter removing which group of interneurons, the line at the bottom of the Figure 2A is always stable, which means that the x i of the BodyWallMuscles is steady. When removing the OtherEndOrgan intermediate neurons, what stands out in our results is the x i of whole neurons (green line) decreases mostly. However, in Figure 2C, for the sexOtherorgan intermediate neurons in male C.elegans, all line corresponding to the sexOtherorgan drop sharply, which both have great effect on the neural network and the output neurons. As a result, we conclude that although these intermediate neurons have no direct control over the output neurons, they still play an important role in the overall functional mechanism of C.elegans. Hence, even though there is no decrease in the activity of the output neurons after the ablation experiment, the average activity of all neurons decreases significantly. In other words, the control mechanism over the output neurons and the whole neurons are not congruent.

Quantitative Analysis of Control Model in Hermaphrodite and Male C.elegans
Ultimately, we try to identify the intermediate neurons group with the greatest control level. We measure the neuron disturbance by degree in the intermediate group with core control, and calculate the effect of neuron control by the change of the mean value of x i . In Figure 3, the neurons in Y-coordinate are arranged in degrees. As can be seen from the Figure 3A, in the case of male C.elegans' sexOtherorgan intermediate group of neurons, the value decreases mostly after the removal of CEPshVR neurons. So CEPshVR has the maximum control over the SBmuscle output neurons group. As shown in Figure 3C, the value in male C.elegans' sexOtherorgan intermediate neurons decreases sharply with the disturbance of Int neurons. Consequently, the Int neurons have the maximum control over the Smuscle output neurons group. Similarly, as shown in Figure 4A, in the case of hermaphrodite C.elegans' OtherEndOrgan intermediate neurons, the HSNR neurons have the maximum control over the sexspecificcells-muscle output neurons group. Briefly, for most of the interneuron groups, the Frontiers in Physics | www.frontiersin.org October 2021 | Volume 9 | Article 731941 higher the degree, the bigger circle in Figures 3, 4, the greater the controlling impact. To illustrate more clearly, we mark in Figures  3B,D, 4B, The gray area in the middle displays the specific neurons of the interneuron group, where the neurons of the core control level are marked in yellow.

DISCUSSION
In our framework, the coupling method of network structure and dynamical interactions enables us to establish a dynamical model of neuron networks to explore the system performance. To quantitatively solve the high-dimension rate equations, we provide a way in which multi-dimensional features between neurons can be reduced into one-dimension equation. Our approach can help the realization of the mapping from structure to function, and the quantitative measurement from function to control for neuron systems. Excepting this, we noticed that x eff of the P.pacificus simulation point is much higher than C.elegans. With its necromenic beetle association, P.pacificus has an intermediate position between the microbivorous C.elegans and true parasites [34,35]. The omnivorous feeder P.pacificus should have the most complex metabolic pathways for nutrition and protection against defense and prey, comparing with the microbivorous C.elegans [36,37]. So we infer that the essence of this phenomenon is the distinction of predation behavior caused the difference of synaptic connection, which changes the structure of nematode network.
Furthermore, we hope that our research could be used to illustrate evolution and degeneration of neuron: M4 most likely degenerates in C.elegans, evolve in P.pacificus. Combining the differences in the two nematode predation strategies, and the functions of M4, our evidence is as follows: M4 are known to be one of the major excitatory neurons in C.elegans, which is required for posterior isthmus peristalsis [38][39][40], it meant that M4 plays a crucial role in the act of eating or swallowing in C.elegans. However, the ablation experiments made by Edelman and Garry have shown that MC are the only neurons required for rapid eating in P.pacificus. But the ablation of other cholinergic pharyngeal neurons, such as M2S and M4, only slightly reduced feeding rates [41,42]. In addition, Trojanowski's previous results demonstrate that this robust and evolutionarily adaptable network is highly degenerate at both the neural and genetic levels, the same behavior can be stimulated by multiple neurons and different types of receptors [41]. Avery and Horvitz even found out, in the C.elegans gene ced-3 Mutant, probably MSpaaaaap (the sister of M4), can sometimes take over M4's function. Using a laser microbeam to kill cells, they found that one of the extra cells in ted-3 worms, tentatively identified as MSpaaaaap, could become a functioning M4 neuron about half of the time, although it rarely or never fully replaced M4 [38]. Given the above, we can not be absolutely sure that M4 most likely degenerates in C.elegans, evolve in P.pacificus, but we hope our theoretical framework can help explore the function of neurons in the evolution and degeneration process.  To illustrate the power of our framework, we analyzed the control model of hermaphroditic and male C.elegans. Based on the data of Cook et al [3], our experiment abstracted the neuron connection group from the real network to the theory level by a nonlinear dynamical model, and obtained a measurable control index of the intermediate neurons to the output neurons. The control between different functional groups need to take into account both structural and functional characteristics. It means that our framework may help design more targeted and accurate biological ablation experiments. To understand vividly, we mapped the control model of the three-layer neurons, which is accurate to the individual neurons [43,44].
In sum, our theoretical research perfects the exploration of the functional structure of nematodes, and breaks the limitation of analyzing the function and control ability from the structure. Comparing with the linear control model to predict motor neuron controllability, Yan et al. use network control principles, we emphasize the growth rate of the node itself and the influence of its neighbors on the node [18,45] Importantly, this is a more efficient way to quickly target groups of neurons, even individual neurons or synapses that have core controlling power. In this view, our theoretical approach contributes to exploring the importance and the evolution mechanism of the individual neurons, and proves how to calculate the control level with the dynamics. These observations are not limited to nematode, but are also applicable to connectome model organism with complete data and clear background. However, the real interactions of nematode connection network are much more complicated than our dynamical model [46,47]. We use this model for theoretical exploration, there may be some difference between our parameter setting and the real system. We need to refer to more biological experiments for further parameter setting and to use more complex dynamical equations or generative models to capture the interactions in our next work.

AUTHOR CONTRIBUTIONS
DD and WX performed the analysis. WX validated the analysis and drafted the manuscript. DD and SS reviewed the manuscript. WX designed this study. All authors have read and approved the content of the manuscript.