Introduction

The ion current flowing through the gap junction between neurons changes the membrane potential of those neurons, which is called signal transmission by electrical synapses.(Carew and Kandel, 1976) Electrical synapses have characteristics such as simple structure, fast signal transmission, bidirectional communication (in most cases), and a wide operation range for voltage including the subthreshold region.(Bennett, 1997) Based on this, it has been reported that electrical synapses play various roles, such as excitation or inhibition of neurons, synchronization or desynchronization of neural activity, promotion or attenuation of rhythms, and improvement of signal-to-noise ratio, etc.(Connors, 2017) Here, rhythm refers to the periodic active pattern of neurons, and in particular, the rhythm in the subthreshold region helps neurons ignite the action potential simultaneously according to the cycle of the rhythm, and a typical example of this can be seen is the excitatory cells in inferior olivary nucleus.(Bazzigaluppi et al., 2012; Benardo and Foster, 1986; Giocomo et al., 2007; Khosrovani et al., 2007; Long et al., 2002) Spontaneous and robust subthreshold membrane potential oscillations were observed in the inferior olivary cells.(Benardo and Foster, 1986; Chorev et al., 2007; Devor and Yarom, 2002; Lampl and Yarom, 1997; Llinas and Yarom, 1981) The inferior olivary cells are intra-connected only by electrical synapses, and it was reported that the presence of electrical synapses is largely involved in the synchronization of subthreshold oscillations.(De Zeeuw et al., 2003; Leznik and Llinás, 2005; Leznik et al., 2002; Llinás, 2014; Long et al., 2002; Turecek et al., 2016) Although electrical synapses and subthreshold oscillations are functionally closely related, theoretical and computational exploration of this has been lacking. The occurrence of subthreshold oscillations and their role in neural networks have been explored through pioneering neuron models,(Roach et al., 2018; Tchumatchenko and Clopath, 2014; V-Ghaffari et al., 2016) but the establishment of an entire connectome between numerous neurons belonging to regions such as the inferior olivary nucleus remains a task to be solved in the future.

In the field of connectome, information on the indexing of all individual neurons, muscles, and organ cells present in the nematode Caenorhabditis elegans (C. elegans) and the entire connectivity between those cells by chemical and electrical synapses was recently updated.(Cook et al., 2019) This available connectome data of C. elegans identifies connectivity weights proportional to the structural size of each synapse, with 1,433 electrical synapses present among a total of 469 cells, including neurons, muscles, pharynx, and other organ cells. Some examples of synchronized activity of a small number of neurons through electrical synapses in C. elegans have been reported. During the gentle nose touch response, nose mechanoreceptor neurons connected to electrical synapses simultaneously increased in activity, thereby serving as a coincidence detector.(Chatzigeorgiou and Schafer, 2011) And a pair of GABAergic motor neurons controlling the bowel movement program repeated every 50 seconds used synchronized activity through electrical synapses at their axon terminal.(Choi et al., 2021) However, direct observation or indirect evidence of spontaneous subthreshold oscillations has NOT been confirmed in C. elegans. Nevertheless, a theoretical and computational model consideration of the signal transmission of subthreshold oscillations on the entire electrical synapse network among cells with various cell types could be a great challenge.

In general, the wave properties of membrane potential are well-known, such as constructive summation of the two excitatory postsynaptic potentials or deconstructive summation of the excitatory and inhibitory postsynaptic potentials. Likewise, it is experimentally confirmed that subthreshold oscillations have wave properties as fluctuations observed in membrane potential measurements and experience interference when propagated in the neural tissue space.(Chiang and Durand, 2023; Gupta et al., 2016) Therefore, interference should be considered important when subthreshold membrane potential waves propagate through the neural network and cross each other at one spatial point.

Here, we attempted to describe a situation in which wave signals, mimicking the subthreshold membrane potential waves, are propagated and interfered on a network that mimics the electrical synapse network of C. elegans (the cells become nodes, the electrical synapses become edges, and the subthreshold waves can flow in both directions on the edges). In the field of physics, the Anderson localization phenomenon has been well known, and it is a general wave phenomenon due to the wave interference between multiple scattering paths. Depending on the strength of the scattering of waves, the result of wave interference could be constructive or destructive, namely, waves could transmit or localize through the media, correspondingly. The theoretical and numerical method for describing the so-called Anderson localization problem usually employed the tight-binding Anderson Hamiltonian. A relevant model was developed to calculate the wave’s quantum percolation using the tight-binding Anderson Hamiltonian on a network lattice consisting of nodes and edges, where the interference by all possible propagation paths was considered when the quantum mechanical probability wave was propagated by experiencing the percolation disorder on the network lattice.(Anderson, 1958; Chang et al., 1995; Meir et al., 1989; Shapir et al., 1982; Thomas and Nakanishi, 2016) In this study, we benefited from this model, and by treating the probability wave as a general wave signal and replacing the network lattice with a synapse network model we are interested in. We could calculate the interference by all possible paths which were considered when the wave signals propagate in the synapse network.

By applying the theoretical concept and computational methodology introduced above and well established in the field of physics, we calculated the signal transmission coefficient between arbitrary two cells of the model connectome as a function of the wavenumber (the quantity of inverse of wavelength) of the wave signal. Based on the signal transmission coefficients for all cell-pairs estimated under various wavenumber conditions, we unraveled 1) the existence of a threshold wavenumber above which the signal transmission disappears, 2) regions of cell-pairs with wavenumber-dependent transmission map, 3) long-distance and regular patterned signal transmission, and 4) major hub cell-pair common across multiple wavenumbers. When comparing with the results for a situation that didn’t consider the wave interference, we revealed the role of the interference when subthreshold waves propagate on the C. elegans’ electrical synapse network model. Although many aspects of the actual neural network and electrical synapses were simplified or omitted in our model study, our results provided the insight of the fundamental role of the electrical synapse network with subthreshold oscillations.

Results

Cell-pairs in the model connectome of C. elegans could belong to regions where the wave signal could be transmitted (or not transmitted) below (or above) a threshold wavenumber

In our circuit model mimicking the anatomical gap junction network of C. elegans (details described in “Construction of our circuit model” part of Methods) (Fig. 1), we calculated all the transmission coefficient Tij of the wave signal for each of the 109,746 cell-pairs 〈i, j〉 between the total 469 cells (details described in “Calculation for the transmission coefficient” part of Methods). Here, the transmission coefficient was a function of energy E of the incoming wave, and the E value was a function of the wavenumber k of the wave signal. And according to the wavenumber, E value was set up to exist in the range [−10, 10], and we calculated all the Tij(E) of the entire cell-pairs at E values of 801 points listed by equal intervals over the entire range (ΔE = 0.025).

Anatomical gap junction network and our circuit model.

(A) In the sample of the network diagram on the left, there are white-filled nodes stand for cells, and the edges of a solid line stand for the anatomical gap junction between them. The graph on the right is the distribution of 1,433 gap junction weights obtained from the connectome data of C. elegans. (B) In the sample of the network diagram on the left, equally spaced two virtual nodes were added between cells connected by a gap junction, and the virtual node is presented as a black-filled circle. The solid line stands for the connection between the cell and the virtual node that gave a weight of 1, and the dotted line stands for the connection between the two virtual nodes that gave our processed weight between 0 and 1. Our processed weights are calculated from the anatomical gap junction weights, and the graph on the right shows its distribution.

First, the average transmission coefficient of all cell-pairs 〈i, j〉 for each E value, 〈Tij(E)〉All ij cellpairs, was obtained and drawn as a function of E (Fig. 2A). The 〈T(E)〉 graph showed an incomplete but quite symmetrical form based on E = 0. And on the real Y-axis scale, high 〈T(E)〉 values were broadly identified in these three ranges of the E values: [−1.575, −1.275], [−0.350, 0.350], and [1.275, 1.575], Also on the log Y-axis scale, several relatively high 〈T(E)〉 values were locally identified in the form of peaks in the remaining ranges of the E values. Meanwhile, on the log Y-axis scale, the signal transmission of all cell-pairs almost disappeared in the E value range at both ends. In the field of solid-state physics, a threshold energy value in which electrical conductivity disappears in the energy band is called a "mobility edge". We were able to define a threshold wavenumber, namely a signal mobility edge, even in the transmission of the wave signal of our circuit model. For example, considering that 〈T(E)〉 = 10−7 is a value that can be obtained when only one cell-pair out of all cell-pairs shows 1% transmission and all others are completely unable to transmit the wave, this 10−7 value can be regarded as an extreme disappearance of signal transmission throughout the entire circuit. Thus, by taking this reference value, E = ±5.650 in our circuit model became a signal mobility edge. If the reference value is set more conservatively lower than now, the absolute value of E corresponding to the signal mobility edges will be larger than now.

Wavenumber-dependent transmission coefficient.

Average transmission coefficient for all cell-pairs as a function of E values (E value is a function of the wavenumber of the subthreshold wave). (A) The upper panel shows the average transmission coefficient graph on the Y-axis of the real scale and the lower panel shows the same graph on the Y-axis of the log scale. The reference value (10-7) of the signal mobility edge is indicated in the lower panel as a dot-dash line. (B) The transmission coefficient of all cell-pairs at E = 0.225 is represented as a heatmap, which is the highest average transmission coefficient in our result. (C) The transmission coefficient of all cell-pairs at E = 5.650 is represented as a heatmap, which is one of the signal mobility edges. (B, C) In heatmap, the X and Y axes are arranged according to the index of 469 cells, and the index follows that of the original connectome data. Instead of displaying all cell names, only seven cell types were indicated.

Next, the transmission coefficient of all individual cell-pairs Tij(E = 0.225) was represented as a heatmap at the E value (or the corresponding wavenumber) where the 〈T(E = 0.225)〉 value was the highest (Fig. 2B). Here, we confirmed that not all cell-pairs transmit the wave signals with even transmission coefficients. Most of the strong transmissions were found in the cell-pairs of intra-body-wall muscles, and a small number of cells belonging to sensory neurons, motor neurons, and other organ cells also showed strong transmission in cell-pairs between them or to body-wall muscles. The rest of the cells except these cells were almost isolated separately without exchanging the wave signals with any cells even at this E value, which was the highest transmitted state on average.

In addition, the transmission coefficient of all individual cell-pairs Tij(E = 5.650) was represented as a heatmap at the E value (or the corresponding wavenumber) which we defined as the signal mobility edge (Fig. 2C). As expected, all cells were isolated separately without exchanging the wave signals with each other, and the entire circuit stop transmitting wave signals of this wavenumber.

Transmission map of subthreshold waves on the electrical synapse network model depended on values of wavenumber

Between the two E values discussed above (E = 0.225 and E = 5.650), the average transmission coefficient of all cell-pairs, 〈T(E)〉, did not simply increase or decrease, but fluctuated along the E value (Fig. 2). Therefore, it was expected that the composition of the cell-pairs, which showed strong transmission, also changes according to the E value. To efficiently investigate the changes in the composition of the cell-pairs, instead of inspecting the Tij(E) heatmaps in all E values, we selected 31 specific E values and tried to analyze the composition of the cell-pairs only in these E values. We chose the specific E values on the following three criteria: (1) the 〈T(E)〉 value was the local maximum (the shape of a peak) in the 〈T(E)〉 graph, (2) the components of the Tij(E) showed greater than 0.5 at least one cell-pair, (3) The E value was positive number (Fig. S1). The reasons why we used these conditions are as follows (1) we assumed that the peak occurred in the 〈T(E)〉 graph due to the emergence of new cell-pairs with strong transmission, (2) we set the reference value for strong transmission to be 0.5 or higher, and (3) the E values where peak occurred in the 〈T(E)〉 graph showed left-right symmetry based on E = 0. (Additional explanations are in “Auxiliary paragraph 1” of Supplementary Materials.)

The transmission coefficients of all individual cell-pairs were represented as a heatmap at each E value (or the corresponding wavenumber) for the 31 specific E values we selected (Left panels in Fig. 3 and S2-9). To closely observe the composition of cell-pairs with strong transmission, only cells (Table S1-3) belonging to the cell-pairs with Tij(E) > 0.5 were selectively exhibited in these heatmaps, not for all 469 cells. To confirm the spatial status of the cell-pairs with strong transmission on the electrical synapse network model, a network diagram was prepared first. Since drawing all 469 cells as a network diagram was complicated and poorly visible, we presented nodes for only 173 cells that belonged to the cell-pairs with Tij(E) > 0.5 in the positive E values. Electrical synapses between cells were represented as edges, and the thickness of the edge proportional to our processed weights (wij) was used. (Additional explanations are in “Auxiliary paragraph 2” of Supplementary Materials.) On the prepared network diagram, we marked the cell-pairs with strong transmission that satisfies Tij > 0.5 (Right panels in Fig. 3 and S2-9).

Wavenumber-dependent transmission map.

The cell-pairs with strong transmission (Tij(E) > 0.5) at (A) E = 0.225, (B) E = 0.800, (C) E = 1.450, and (D) E = 3.875 (or the corresponding wavenumbers) are shown, respectively. The left panels represent the transmission coefficient between cells belonging to the cell-pairs with strong transmission at each E value as a heatmap. Here, the X and Y axes are the cell index, and the cell name corresponding to the cell index can be found in Table S1-3. The right panels display the cell-pairs of strong transmission at each E value as red bidirectional arrows of the same thickness on the network diagram. In the network diagram, the cell name is written inside the circle that stands for each node (for the body-wall muscles, the cell name is abbreviated as follows; ex) dBWML1 → dL1), and the color of the node stands for the cell type (Red: Pharynx cells, Orange: Sensory neurons, Yellow: Inter neurons, Green: Motor neurons, Light blue: Body-wall muscles, Purple: Other end organs, Magenta: Sex-specific cells). Edges indicated by black solid lines stand for the electrical synapses among the cells, and the thickness of the edges is proportional to our processed weight. In the network diagram, only 173 cells that belonged to the cell-pair with strong transmission at least once in positive E values are represented, and the remaining cells and the electrical synapses by them are omitted from display. The virtual nodes of our circuit are also omitted from the display.

As expected above, the region of the cell-pairs with strong transmission changed depending on the value of E (or wavenumber), which we named "Wavenumber-Dependent Transmission Map (WDTM) of subthreshold waves on the electrical synapse network model". However, it should be also noted that our WDTM was a result derived from a network where only electrical synapses were considered, not a general neural network where chemical and electrical synapses coexist.

We looked at the cell type of the cell-pairs that make up the 31 WDTMs (Fig. 3 and S2-9). In the two E value ranges of [0, 0.350] and [1.275, 1.575], the high 〈T(E)〉 value was broadly identified (Fig. S1), and many cell-pairs of intra-body-wall muscles were observed in the corresponding WDTMs of these E values (Fig. 3A, 3C, S2, S3A-B, and S5). Since the body-wall muscles consist of four muscle strands: dorsal/ventral-left/right strands, the cell-pair of intra-body-wall muscles can be divided into two cases: a cell-pair of intra-strand and a cell-pair of inter-strand. Following this classification, the cell-pairs of both intra-strand and inter-strand were found in the WDTMs of the E value range of [0, 0.350] (Fig. 3A, S2, and S3A-B), whereas only the cell-pairs of intra-strand were found in the WDTMs in the E value range of [1.275, 1.575] (Fig. 3C and S5). On the other hand, in the remaining WDTMs, a small number of cell-pairs within cells with less than ten members were mostly found (Fig. 3B, 3D, S3C-D, S4, and S6-9), and the WDTMs of the two largest cases were shown in Fig. 3B and 3D. The composition of cell-pairs with various combinations of cell types (sensory/inter/motor neurons, body-wall/sex-specific muscles, pharynx, and other organ cells) was shown in these remaining WDTMs.

Cell-pairs in body-wall muscles on long-distance with regular patterned transmission were found in the WDTMs

When looking at each WDTM on the network diagram, it shows various patterns spatially. In particular, in the case of the body-wall muscles, the cell-pairs of intra-strand existed, ranging from short-distance to the longer distance than the half of the full length of the strand (Fig. 3A and 3C). Here, the "Distance" was meant by the length of the shortest path along the network between two cells of a cell-pair. If two cells are directly connected to an edge, the distance between the two cells is one. In the WDTMs of the E value range of [1.275, 1.575], cell-pairs with regular patterned transmission of the body-wall muscles were observed through the network diagram (Fig. S5). Only the cell-pairs of intra-strand were found in these WDTMs, and the configuration of the distance of these cell-pairs was different for each WDTM (multiple of 3 at E = 1.350, multiple of 2 at E = 1.450, multiple of 7 at E = 1.500, and multiple of 4 at E = 1.575).

On the other hand, when the cell-pairs of inter-strand of the body-wall muscles appeared in the WDTMs as shown in Fig. 3A-B, the head mesodermal cell (hmc) seemed to play an important role in building a bridge between them. This would be recognized from the network diagram itself of black edges without WDTM, and hmc was strongly connected near the neck of the body-wall muscles (seventh and eighth muscles) for all four strands. In the WDTM of Fig. 3A, there were many cell-pairs of inter-strand between dorsal-right and ventral-left strands, and there were also many cell-pairs between the ventral-left strand and hmc, which directly revealed the role of hmc as a bridge between the two body-wall muscle strands. In the WDTM of Fig. 3B, it was shown that the role of hmc was indirectly revealed as the cell-pairs of inter-strand, which existed only among the neck of four body-wall muscle strands.

Additionally, we reaffirmed the important role of hmc while confirming the distance characteristics of the cell-pairs with strong transmission. Under the conditions used when preparing the network diagram (the transmission coefficient Tij(E) of a cell-pair was 0.5 or higher in the range of positive E values at least once), a total of 146 cell-pairs with strong transmission were identified excluding the cell-pairs of intra-body-wall muscles, and their distance was investigated (Table S4). (Additional explanations are in “Auxiliary paragraph 3” of Supplementary Materials.) Most of the cell-pairs in Table S4 showed the kinds of short-distances 1 to 3, which consisted of various cell types. On the other hand, the cell-pairs with long-distances in Table S4 ranging from 4 to 17 were mainly pairs between hmc and body-wall muscles. The strong transmission between hmc and the tail muscle, which corresponds to the longest distance, passes through the intra-strand connections from the neck to the tail of body-wall muscles.

A few cell-pairs of strong transmission with long-distance in Table S4 turned out to be not those between hmc and the body-wall muscles. In there, one inter neuron (AVBR) and two motor neurons (AS11, VD04) formed a cell-pair with the body-wall muscles, respectively, and four sensory neurons (CEPDL, IL1R, IL1DR, OLQDL) formed cell-pairs between them.

Major hub cell-pairs with the strong transmission of signals for many wavenumbers exist

We looked at cell-pairs frequently appearing in all 31 WDTMs in order to identify the major hub cell-pairs that is commonly used for various wavenumbers. First, for each cell-pair 〈i, j〉, we calculated the average appearance rate in all 31 WDTMs, 〈Θ[Tij(E) − 0.5]〉E∈{31 selected E}, and represented this as a heatmap (Fig. 4). Where, the Θ[x] was a step function, which value was 1 if x ≥ 0 or 0 if x < 0. Thus, the average appearance rate was 1 in the case of a cell-pair that appeared in all 31 WDTMs, and 1/31 ≈ 0.03 in the case of a cell-pair that appeared only once. In this study, we declared a cell-pair as a major hub cell-pair if the average appearance rate of a cell-pair was 3/31 ≈ 0.1 or higher. The major hub cell-pairs were the cell-pairs of intra-strand of the body-wall muscles (Fig. 4). As seen in the network diagram, the edges of each strand of the body-wall muscles from the neck (seventh or eighth muscle) to the tail (twenty-third or twenty-fourth muscle) were simply connected in a line, unlike the messy edges among most neurons in the network (Fig. 3). In general, the complexity of connections among most neurons increase with the number of transmissible paths of wave propagation, and the sum of phase changes coming from various path-length differences causes the deconstructive interference. Therefore, the simple (or regular) connections in each strand of the body-wall muscles from the neck to the tail reduced the number of transmissible paths of wave signals, which suppressed deconstructive interference and thus resulted in the strong transmission for various wavenumbers.

Average appearance rate as cell-pair with strong transmission in the WDTMs.

The average appearance rate of all cell-pairs is represented as a heatmap. Here, the X and Y axes are arranged according to the index of 469 cells, and the index follows that of the original connectome data. Instead of displaying all cell names, only seven cell types were indicated. The highest rate is 7/31 ≈ 0.23, and four cell-pair groups showing rates of 3/31 ≈ 0.1 or higher are identified and numbered in the heatmap. (1) PVDL-PVDR pair, (2) hmc-vBWML6, hmc-vBWML7, and hmc-vBWML8 pairs, (3) CANL-CANR, CANL-exc_cell, and CANR-exc_cell pairs, (4) many cell-pairs of intra-strand of the body-wall muscles.

Also, the cell-pairs between hmc and the neck of the ventral-left strand of the body-wall muscles were also identified as major hub cell-pair (Fig. 4). As seen in the network diagram, since the edges between hmc and the neck muscles were very strong compared to other edges (Fig. 3), the relative contribution of wave signals coming from other cells could be almost ignored, and thus resulted in the strong transmission for various wavenumbers.

The cell-pairs among two sensory neurons (PVDL, PVDR) and three other organ cells (CANL, CANR, exc_cell) were also identified as major hub cell-pairs (Fig. 4). The cell-pair between CANL and CANR was the most major hub cell-pair, with the highest average appearance rate of 7/31 ≈ 0.23. As seen in the network diagram, the three other organ cells were very strongly connected by the edges to each other, and the two sensory neurons were very strongly connected by the edges to both ALA inter neuron and hmc (Fig. 3). Although there were cases where they were strongly connected by the edges to each other, such as between PVDL and CANL or between PVDR and CANR, they were not major hub cell-pair. Therefore, the strong edge shown on the network diagram was not a sufficient condition to implicate a major hub cell-pair.

Discussion

The fundamental aspect of this study is that the interference phenomena of the subthreshold oscillations propagating on a neuronal circuit was described by a theoretical and computational model study of the reference system C. elegans for which the structural anatomical connectome was completely known. We investigated the wavenumber-dependent transmission of the sinusoidal wave signals between all cell-pairs on the connectome model of C. elegans by considering the interference effect caused by multi-paths of signaling.

The results related to the wavenumber-dependent transmission map expanded the way we used to view signal transmission on electrical synapses. The difference in signal transmission with the presence or absence of interference effect on electrical synapses could be illustrated by considering the effective conductance between two nodes in the electric circuit with or without considering the interference effect (Fig. 5). For a situation without considering the interference, we calculated the effective resistance and effective conductance for all cell-pairs by imitating a cell as a node and the electrical synapse as an electrical resistor connecting nodes as shown in Fig. 5A (details in “Calculation for effective conductance” part of Methods). Several inter neurons and other organ cells showed high effective conductance, while pharynx, body-wall muscles from the neck to the tail, and sex-specific muscles showed very low effective conductance (Fig. 5B). Since inter neurons have many complex connections among them, many electrical current paths exist between the two inter neurons, which act as parallel connections of multiple electrical resistors. This results in lower effective resistance (or higher effective conductance). However, for a situation with considering the interference, many transmissible forward and backward paths of the signal propagation rather caused the deconstructive interference, preventing the transmission of the signal. For each cell-pair 〈i, j〉, we calculated the average transmission coefficient for all E values and represented this as a heatmap (Fig. 5C-D). Which clearly showed that most of the neurons could not exchange wave signals with each other except for very few neurons. Wave signals were effectively delivered only within a small number of very strongly connected cell groups (PVDL/R, CANL/R, exc_cell, hmc) or within a large number of very regularly connected cell groups (body-wall muscle strands from the neck to the tail).

Differences depending on whether the electrical synapse network model considers interference.

(A) In the electrical synapse network model without considering interference, when the electrical synapse is treated as an electrical resistor, the effective conductance (the reciprocal of the effective resistance) experienced by the external direct current power applied to an arbitrary cell-pair was calculated, and this value was represented as (B) a heatmap for all cell-pairs. (C) In our circuit considering interference of wave signals, the average transmission coefficient for all E values (or all wavenumbers) of an arbitrary cell-pair was calculated, and this value was represented as (D) a heatmap for all cell-pairs. (B, D) In heatmap, the X and Y axes are arranged according to the index of 469 cells, and the index follows that of the original connectome data. Instead of displaying all cell names, only seven cell types are indicated.

In this study, we used the structural anatomical information from C. elegans’ electrical synapse network, but no synchronized rhythmic activity induced by subthreshold membrane potential oscillations has been reported for C. elegans so far. However, we aimed at investigating the signaling characteristics caused by the interference and to discover a feasible phenomenon related to the propagation of the subthreshold oscillations on a neuronal circuit of living nervous system. As to the synchronized rhythmic activation observed like in the mammalian inferior olive nucleus, we think that there might be an electrical synapse network in there that connect cells very strongly or very regularly. The plausible possibility according to our model study is that the constructive interference of subthreshold membrane potential waves with a specific wavenumber may generate the synchronized rhythmic activation. We hope that the results in our study would serve as the worthwhile framework and knowledge for designing the future experimental studies not only for inferior olive nucleus, C. elegans but also other living systems.

Methods

Construction of our circuit model

We used the contents of the "hermaphrodite gap jn symmetric" tab in the "SI 5 Connectome adjacency matrices, corrected July 2020.xlsx" file of the C. elegans connectome dataset as data for the anatomical gap junction network.(Cook et al., 2019) The 469 cell names and the gap junction weight (gij) of all cell-pairs connected by gap junction were provided in the tab, and the weights had a value distribution from a minimum of 1 to a maximum of 401 (Fig. 1A). The weights were proportional to the spatial size of the corresponding gap junction, which we hypothesized that the larger the weight, the stronger the connection between the cells and thus the better transmit the wave signal to the connected cell.

Our circuit model described C. elegans’ anatomical gap junction network as follows (Fig. 1B). First, each cell was treated as a node in our circuit. Each gap junction connected to a cell-pair was treated as an edge connecting two nodes in our circuit, and we inserted two virtual nodes (VNs) between those two nodes to be connected sequentially through VNs (Cell-VN-VN-Cell). Here, all edges inside our circuit were set to have the identical length a. The edge of Cell-to-VN (or VN-to-Cell) representing like axon or dendrite of neuron gave a weight of 1, and the edge of VN-to-VN representing a gap junction gave our processed weight (wij = min{1, gij⁄40}). Our processed weight (wij) was proportional to the anatomical gap junction weight (gij) and normalized as a value between 0 and 1 but was fixed to 1 if gij > 40. In our circuit model, the weight of the edge meant the transmission coefficient of the signal transmitted through this edge, so if the weight was 1 corresponded to complete transmission, and 0 corresponded to complete reflection (or unavailable edge), respectively. And within our circuit, the characteristic of the signal being spontaneously attenuated, and disappearing was excluded. Therefore, when a signal incoming from the outside to a node in our circuit was reflected to the outside at the same node and simultaneously transmitted to the outside at another node, the sum of the transmission and reflection coefficients to the outside always remained at 1 without loss.

Calculation for the transmission coefficient

We benefited from the method of calculating the transmission coefficient according to the energy E of the incoming wave when the electron probability wave incident on one point of the network lattice was detected on another point. This methodology, called the quantum percolation model, considered the propagation and interference of waves on the network lattice using the tight-binding Anderson Hamiltonian.(Anderson, 1958; Chang et al., 1995; Meir et al., 1989; Shapir et al., 1982; Thomas and Nakanishi, 2016) In our methodology, we replaced the classical wave signals of complex number expression representing subthreshold membrane potential waves instead of the electron probability wave and replaced the network lattice with our circuit describing the electrical synapse network model built above. Then, the transmission coefficient calculated by the quantum percolation model told how well the subthreshold waves were transmitted from one cell i to another cell j while undergoing interference by multiple propagation paths in the electrical synapse network model. Here, the electrical synapses act as wave scatters that interrupt the propagation of the subthreshold waves from cell i to cell j.

In the quantum percolation model, the tight-binding Anderson Hamiltonian was defined as follows,

Where the |m⟩ (also expressed as the ψm) is as the electron probability wave on the m-th node, and the ⟨m| means its complex conjugate . We used ψm as the classical wave signal observed on the m-th node, which of complex number expression made it easy to express the phase shift of the wave and calculate the interference. The amplitude of this classical wave signal on the m-th node was defined as the In the Hamiltonian, the εm is the on-site energy, and we regarded it as the resting membrane potential of the m-th cell and set it to 0, the identical value as the relative reference value of the resting membrane potential in all cells. The 〈mn〉 represents a pair of nearest-neighbor sites in the Hamiltonian, we replaced it with a node-pair connected by the edges in our circuit. The Vmn as a hopping matrix represents the entire structure of the network lattice in the Hamiltonian, which we implemented the structure of our circuit by setting it to Vmn = 1 for a node-pair between cell and VN, and Vmn = wmn for that between two VNs. Since the total number of the cells was 469, and we inserted two VNs for every 1,433 electrical synapses, the newly constructed tight-binding Anderson Hamiltonian for our circuit was represented as a matrix of 3,335 by 3,335.

By solving the time-independent Schrödinger equation for the Hamiltonian of our circuit, we sought to obtain the reflection and transmission coefficients of the cell-pair between the IN-cell and the OUT-cell from the ψIN of the IN-cell where the wave signal from the outside was incident and reflect, and the ψOUT of the OUT-cell where the wave signal from our circuit transmitted to the outside. We defined ψIN and ψOUT as follows using complex numbers r and t with amplitudes and phases,

Where the ψIN was the sum of 1 and r corresponding to the incoming wave and the reflective wave, respectively, and the ψOUT was t corresponding to the transmitted wave. The reason why the amplitude and the phase of the incoming wave on the IN-cell were 1 and 0, respectively was because we set that as a reference amplitude and a reference phase in our circuit. The amplitude of the wave reflected from the IN-cell to the outside became the reflection coefficient with R = rr, and the amplitude of the wave transmitted from the OUT-cell to the outside became the transmission coefficient with T = tt, and T + R = 1 had to be always satisfied, as mentioned above.

This methodology also assumed both an external injector and an external detector. The external injector injected the generated wave into the network lattice and received the reflective wave from the network lattice. The external detector received the transmitted wave from the network lattice. The waves on the external injector ϕInjector and the external detector ϕDetector were extended from the definitions of ψIN and ψOUT above, respectively, and determined as follows,

The external injector and detector were set to be connected to the IN-cell and OUT-cell by a perfect conducting wire with a length of a, respectively. The wave outgoing to the external injector (detector) from the IN-cell (OUT-cell) was determined by multiplying the phase difference of eika to the wave on the IN-cell (OUT-cell), whereas the wave incoming to the IN-cell from the external injector was determined by multiplying the phase difference of eika to the wave on the IN-cell. Insert these terms into the time-independent Schrödinger equation for the Hamiltonian of our circuit, the following equation was derived describing the entire system consisting of our circuit and the external injector and detector,

Where the E value on the right side representing the energy of the entire system was equivalent to the energy E of the incoming wave from the external injector.

In the quantum percolation model, the energy E of the incoming wave was defined as a tight-binding energy function as follows,

Where γ meant the binding energy between the nearest-neighbor atoms in the tight-binding model, but we had to set it to an arbitrary value. In this study, it was set to γ = 5, and the decision process was described in the below section. Applying this energy function, the above equation for the entire system was expressed as the matrices as follows,

Where m and n meant the 3,333 remaining nodes except two nodes for IN-cell and OUT-cell. For a wavenumber k in the range of [0, π/a] (or an E value in the range of [−2γ, 2γ]), all components of both the left-hand square-matrix and the right-hand column-vector were fully expressed as complex numbers, and then the r and t values in the left-hand column-vector were obtained by matrix multiplication of the inverse matrix of the left-hand square-matrix and the right-hand column-vector. From squared absolute the r and t values, the reflection and transmission coefficients of the cell-pair between IN-cell and OUT-cell at the E value (or the corresponding wavenumber k) were deduced, respectively.

Searching for optimal gamma

To determine the appropriate γ of the above energy function, we preliminarily calculated the transmission coefficient for all cell-pairs between the 469 cells under all four conditions of γ, and compared the average transmission coefficient graph as a function of the E value 〈Tij(E)〉All ij cellpairs (Fig. S10). Where four values of 1, 2.5, 5, and 10 were attempted for γ, and equally spaced a hundred values were used for the wavenumber k in the range of [0, π/a]. According to the above energy function, the minimum/maximum range of the E value differed as [−2γ, 2γ], and as γ was larger, the points of the 〈T(E)〉 graph were placed sparsely on the X-axis of the E value. In each of the 〈T(E)〉 graphs for these four conditions of γ, similarly high 〈T(E)〉 values occurred in the vicinity of E = 0, ±1.4. Based on this preliminarily estimation, we confirmed that the occurrence of high 〈T(E)〉 values appeared as a function of the E value regardless of γ. Therefore, we considered γ as just determining the minimum/maximum range of the E value.

The appropriate γ we wanted in this study was large enough to provide a wide E value range to check the signal mobility edge in the 〈T(E)〉 graph, and the appropriate γ had to be small so that the points in the 〈T(E)〉 graph did not become too sparse. It was the condition of γ = 5 that satisfied this demand out of the four conditions of γ. Because the condition of γ = 2.5 was not sufficient to confirm the signal mobility edge, and the condition of γ = 10 had too sparse points along the X-axis.

Error report on our calculation for the transmission coefficient

In principle, the sum of transmission coefficient Tij(E) and reflection coefficient Rij(E) of any cell-pair at any E value always had to satisfy 1 because our circuit model did not consider spontaneous attenuation of the wave signal. However, in the calculating program we wrote and used for this study, the error occurred with the sum of the two coefficients exceeding 1 (by the tolerance was 10-5) for a total of 5,213 cell-pairs in three specific E values (45 cell-pairs at E = −1.000, 3,836 cell-pairs at E = −0.375, and 1,332 cell-pairs at E = 1.000). We performed the calculations on a total of 109,746 cell-pairs for a total of 801 E values, so these errors correspond to 0.006% of the total calculation results. To eliminate the effect of these errors in this study, we post-corrected the two coefficients of the 5,213 cell-pairs in the three specific E values in which the error occurred to Tij(E) = 0, Rij(E) = 1.

Calculation for effective conductance

An electrical resistor network model was constructed that mimics C. elegans’ electrical synapse network, with the 469 cells replaced with electrical conducting nodes and the 1,433 electrical synapses replaced with electrical resistors connecting these nodes. The effective resistance and its reciprocal, the effective conductance, when a direct current power source outside the electrical resistor network model contacted a cell by its positive pole, and another cell by its negative pole, respectively, were calculated.(Klein and Randic, 1993) The effective conductance shows how well electric current without interference flows between the two cells contacted to the external power source.

We expressed the electric voltage at node i as vi, the electric current flowing from node i to node j as μij, the resistance of the electrical resistor between node i and node j as χij, the conductance of that as . The anatomical gap junction weights (gij) were used as the conductance value of the electrical resistors in the network, ωij = gij, with arbitrary units.

By Ohm’s law, each electrical resistor in the network satisfies the following equation,

The sum of the electric currents entering each node by Kirchhoff’s current law satisfied the equation below,

Where (+) and (−) meant a node in which the positive and negative poles of the external power source were in contact, respectively, and the capital letter I was the total electric current supplied by the external power source. When Ohm’s law above was substituted on the left-hand of this equation, it was expressed as viji ωij − ∑ji ωijvj, and then the equation was able to express in matrix and vector as follows,

Where Lω was the Laplacian matrix of the electrical resistors’ conductance matrix (ω), defined as was a column-vector with an electric voltage at all nodes as a component, and was a unit column-vector of the identical dimension as , with only the component of the (+) node ((−) node) having a value of 1 and all other components having a value of 0. By applying the pseudo-inverse matrix of the Laplacian matrix (Lω+) on both sides of the equation, we obtained a particular solution as .

The effective resistance applied between the positive and negative poles of the external power source was defined by Ohm’s law of the electric voltage difference between the (+) and (−) nodes and the total electric current, as follows,

Substituting here the particular solution of obtained above was as follows,

As a result, we were able to obtain the effective resistance for any one cell-pair by calculating the pseudo-inverse matrix of the Laplacian matrix (Lω+) from the electrical resistors’ conductance matrix (ω). The effective conductance was defined as the reciprocal of the effective resistance, Geff = (Reff)−1.

Data and materials availability

All data are available in the manuscript or the supplementary materials.

Additional information

Acknowledgements

We appreciate two Nobel Laureates, Prof. Erwin Neher and Prof. Kurt Wuthrich, for their insightful discussions on the impact of this work. We have benefited a lot from their suggestions in the course of this study.

Funding

This work was supported in part by both pCoE program of DGIST [grant number: 23-CoE-BT-01] and a grant from iPT, Korea.

Author contributions

IC and SK conceived the idea for this work. TC, IC, and SK set up the model, conducted the calculation, and analyzed the data and results. IC and SK wrote the manuscript. All authors reviewed and agreed on the contents of the manuscript.

Declaration of interests

The authors declare no competing interests.

Supplementary Materials

Auxiliary paragraph 1

The 31 specific E values in this study were not fixed and were of a property that can change as a researcher adjusts the resolution of the E value. Perhaps the higher the resolution of the E value, the more specific E values will be found. And it was arbitrary that our reference value for strong transmission was 0.5, and this reference value can be raised or lowered depending on the researcher’s intention, and the lower the reference value, the greater the number of cell-pairs with strong transmission will be found.

Auxiliary paragraph 2

Because of the omitted the rest of cells and their electrical synapses in the network diagram, some nodes in the diagram look like remote islands, but it should be remembered that the 173 cells were connected as one large cluster through both the drawn and the omitted electrical synapses.

Auxiliary paragraph 3

Because the number of the cell-pairs of intra-body-wall muscles was overwhelmingly large and it was relatively easy to measure the distance compared to those of other cell-pairs, the investigation for the distance excluded them (for example, the distance of two body-wall muscles in the same strand easily calculate from the given number in the name of the muscles).

Selected 〈T(E)〉 peaks.

The 31 selected peaks are shown as dots on the graph of the average transmission coefficient for all cell-pairs. To make the dots more visible, only the E value range of [−1, 5] was drawn out of the total [−10, 10]. The upper panel shows the dots and the graph on the Y-axis of the real scale and the lower panel shows the same things on the Y-axis of the log scale.

Wavenumber-dependent transmission map.

The cell-pairs with strong transmission (Tij(E) > 0.5) at (A) E = 0.000, (B) E = 0.100, (C) E = 0.175, and (D) E = 0.225 (or the corresponding wavenumbers) are shown, respectively. The left panels represent the transmission coefficient between cells belonging to the cell-pairs with strong transmission at each E value as a heatmap. Here, the X and Y axes are the cell index, and the cell name corresponding to the cell index can be found in Table S1-3. The right panels display the cell-pairs of strong transmission at each E value as red bidirectional arrows of the same thickness on the network diagram. In the network diagram, the cell name is written inside the circle that stands for each node (for the body-wall muscles, the cell name is abbreviated as follows; ex) dBWML1 → dL1), and the color of the node stands for the cell type (Red: Pharynx cells, Orange: Sensory neurons, Yellow: Inter neurons, Green: Motor neurons, Light blue: Body-wall muscles, Purple: Other end organs, Magenta: Sex-specific cells). Edges indicated by black solid lines stand for the electrical synapses among the cells, and the thickness of the edges is proportional to our processed weight. In the network diagram, only 173 cells that belonged to the cell-pair with strong transmission at least once in positive E values are represented, and the remaining cells and the electrical synapses by them are omitted from display. The virtual nodes of our circuit are also omitted from the display.

Wavenumber-dependent transmission map.

The cell-pairs with strong transmission (Tij(E) > 0.5) at (A) E = 0.275, (B) E = 0.325, (C) E = 0.475, and (D) E = 0.575 (or the corresponding wavenumbers) are shown, respectively. The heatmaps of the left panels and the network diagram of the right panels are represented in the same way as in Fig. S2.

Wavenumber-dependent transmission map.

The cell-pairs with strong transmission (Tij(E) > 0.5) at (A) E = 0.750, (B) E = 0.800, (C) E = 0.850, and (D) E = 0.925 (or the corresponding wavenumbers) are shown, respectively. The heatmaps of the left panels and the network diagram of the right panels are represented in the same way as in Fig. S2.

Wavenumber-dependent transmission map.

The cell-pairs with strong transmission (Tij(E) > 0.5) at (A) E = 1.350, (B) E = 1.450., (C) E = 1.500, and (D) E = 1.575 (or the corresponding wavenumbers) are shown, respectively. The heatmaps of the left panels and the network diagram of the right panels are represented in the same way as in Fig. S2.

Wavenumber-dependent transmission map.

The cell-pairs with strong transmission (Tij(E) > 0.5) at (A) E = 1.700, (B) E = 1.750, (C) E = 1.800, and (D) E = 1.875 (or the corresponding wavenumbers) are shown, respectively. The heatmaps of the left panels and the network diagram of the right panels are represented in the same way as in Fig. S2.

Wavenumber-dependent transmission map.

The cell-pairs with strong transmission (Tij(E) > 0.5) at (A) E = 2.000, (B) E = 2.225, (C) E = 2.450, and (D) E = 2.650 (or the corresponding wavenumbers) are shown, respectively. The heatmaps of the left panels and the network diagram of the right panels are represented in the same way as in Fig. S2.

Wavenumber-dependent transmission map.

The cell-pairs with strong transmission (Tij(E) > 0.5) at (A) E = 3.000, (B) E = 3.375, (C) E = 3.475, and (D) E = 3.525 (or the corresponding wavenumbers) are shown, respectively. The heatmaps of the left panels and the network diagram of the right panels are represented in the same way as in Fig. S2.

Wavenumber-dependent transmission map.

The cell-pairs with strong transmission (Tij(E) > 0.5) at (A) E = 3.625, (B) E = 3.750, and (C) E = 3.875 (or the corresponding wavenumbers) are shown, respectively. The heatmaps of the left panels and the network diagram of the right panels are represented in the same way as in Fig. S2.

Searching for optimal gamma.

The average transmission coefficient for all cell-pairs as a function of E value (E value is a function of both γ and wavenumber k) was calculated under the four γ conditions (1, 2.5, 5, 10), respectively. Here, the wavenumbers used the same set as a hundred values equally spaced in the range of [0, π/a]. The upper panel shows the average transmission coefficient graphs as the Y-axis of the real scale and the lower panel shows the same graphs as the Y-axis of the log scale. The reference value (10-7) of the signal mobility edge is shown in the lower panel as a dot-dash line.

Cell names corresponding to the indices in the heatmaps for wavenumber-dependent transmission map.

Each two columns shows the index and the corresponding cell name of the cells that belonged to the cell-pairs with strong transmission at the E value (or the corresponding wavenumber) indicated in the first row. The color painted on the name box stands for the cell type (Red: Pharynx cells, Orange: Sensory neurons, Yellow: Inter neurons, Green: Motor neurons, Light blue: Body-wall muscles, Purple: Other end organs, Magenta: Sex-specific cells).

Cell names corresponding to the indices in the heatmaps for wavenumber-dependent transmission map.

It is represented in the same way as Table S1.

Cell names corresponding to the indices in the heatmaps for wavenumber-dependent transmission map.

It is represented in the same way as Table S1.

List of cell-pairs with strong transmission except intra body-wall muscles pairs.

The first column represents the shortest distance along the network connecting the two cells of the cell-pair. The second and third columns represent the cell names of the cell-pair, and the color painted on the name box stands for the cell type (Red: Pharynx cells, Orange: Sensory neurons, Yellow: Inter neurons, Green: Motor neurons, Light blue: Body-wall muscles, Purple: Other end organs, Magenta: Sex-specific cells). The fourth, fifth, and sixth columns represent the transmission coefficient, the E value, and the wavenumber when the cell-pair shows the strongest transmission in the positive E value range, respectively.