Autapses promote synchronization in neuronal networks

Neurological disorders such as epileptic seizures are believed to be caused by neuronal synchrony. However, to ascertain the causal role of neuronal synchronization in such diseases through the traditional approach of electrophysiological data analysis remains a controversial, challenging, and outstanding problem. We offer an alternative principle to assess the physiological role of neuronal synchrony based on identifying structural anomalies in the underlying network and studying their impacts on the collective dynamics. In particular, we focus on autapses - time delayed self-feedback links that exist on a small fraction of neurons in the network, and investigate their impacts on network synchronization through a detailed stability analysis. Our main finding is that the proper placement of a small number of autapses in the network can promote synchronization significantly, providing the computational and theoretical bases for hypothesizing a high degree of synchrony in real neuronal networks with autapses. Our result that autapses, the shortest possible links in any network, can effectively modulate the collective dynamics provides also a viable strategy for optimal control of complex network dynamics at minimal cost.

An autapse is referred to as a special synapse that connects the axon and dendrites of the same neuron. Autapses were first discovered in 1972 in the pyramidal cell of neocortex cerebri (Golgi preparations of rabbit occipital cortex) 1 . Because of their aberrant and incestuous structure, autapses were originally regarded as "a wiring error" in the development of the brain and therefore were conceived as playing no actual role in the biological functions of the underlying neural network. Twenty-five years later, autapses were found commonly in the neural circuits of the cat's visual cortex 2 , suggesting that autapses may actually have certain biological usage 3,4 . Sign of ubiquity of autapses emerged subsequently, when they were discovered in multiple areas in the brain such as neocortex, hippocampus, cerebellum, substantia nigra, and striatum [5][6][7][8] . There was experimental evidence that autapses are sparse and are present only on a small fraction of the neurons in the underlying brain neuronal network 5,6 . With respect to the functional role of autapses, significant improvement in the spike-time precision of the neocortical interneurons was found in the presence of inhibitory autapses 6 , and excitatory autapses were found to be essential to persistent activities of certain neurons of Aplysa 8 . In theoretical and computational studies, issues that have been addressed concerning the roles of autapses include: persistence in recurrent neural networks 9 , generation of oscillatory behavior of a single neuron 10 , switching among distinct dynamical states (e.g., quiescent, periodic or chaotic) 11 , induction of wave patterns in a regular network of neurons 12 , signal detection in stochastic neurons 13 , emergence of coherence resonance in single neurons and in neural networks 14 , and promotion of rhythmic propagation in neuronal networks 15 .
A fundamental issue in complex networks with significant implications to biology and physiology is synchronization 16 , a ubiquitous phenomenon in natural systems 17 . A widely known example is epilepsy, a disorder characterized by seizures, which affects over 50 million people worldwide. The classical understanding is that seizures are associated with hypersynchrony 18 . This principle, however, has been challenged [19][20][21][22] , and the dynamical origin of epilepsy in relation to synchrony has been continuously debated 23,24 . As a matter of fact, the lingering and unanswered question in this field of medicine is whether epileptic seizures are associated with an increased or a decreased level of neuronal synchrony. To certain extent, this question may be addressed through the approach of data analysis: by analyzing EEG (electroencephalogram) or ECoG (electrocorticogram) data recorded from epileptic patients and calculating measures of synchrony, one hopes to be able to assess the possible causal role of synchronization in triggering a seizure. Indeed, there were previous efforts in developing seizure detection and characterization frameworks based on partial synchrony among multichannel brain data, with the finding that, depending on the type of seizures, there can be either an enhancement or a reduction in the degree of Figure 1. The impact of a single autapse on synchronization in a toy neuronal network. (a1) Without any autapse, the network has four nodes and four edges, where each node is a Hindmarsh-Rose neuron. (b1-d1) Network structure when a single autapse (represented by the red curve with an arrow) is present at node 1, 2, and 4, respectively. (a2-d2) For the network structures in (a1-d1), respectively, the dynamical behaviors of the network in terms of synchronization. Shown in each panel is a plot of the x variable from each node versus the averaged value of this variable over all the nodes during the time evolution. When there is global synchronization, all the variables are equal to their average value at any instant of time, tracing out a straight line segment along the diagonal. Any deviation from the diagonal signifies lack of synchronization. The uniform coupling parameter is ε = 1 and the time delay associated with the autapse is τ = 4.
Scientific REPORts | (2018) 8:580 | DOI:10.1038/s41598-017-19028-9 where x is the membrane potential, y and z represent the transport rates of the fast and slow ion channels corresponding to the spiking and bursting variable, respectively, I is the external stimulating current that determines the nature of the isolated neuronal dynamics, e.g., periodic or chaotic. To appreciate the ability of autapses to promote synchronization for complicated dynamics, we choose the parameters of the HR neuron so that it exhibits a chaotic attractor (see Methods). A single autapse can be present at any node, and there are three distinct cases, as shown in Fig. 1(b1-d1), respectively (due to the symmetry in the network structure, an autapse on node 3 is equivalent to one on node 2). To determine if there is global synchronization in the network, we plot one dynamical variable, e.g. x i (t), from each node (for i = 1, …, 4), versus the network averaged variable defined as x ave (t) = ∑ i x i (t)/4 for a relatively large time interval. When the network is completely synchronized, we have x i (t) = x ave (t) for i = 1, …, 4, so the dynamical trajectory will trace out a straight line segment along the diagonal in the (x i , x ave ) plane. Any deviation from the diagonal indicates lack of synchronization in the system. From Fig. 1(a2-d2), we see that, without any autapse (a2), the network is not synchronized. When an autapse is present at node 1 or 4, which has the largest and the smallest degree, respectively, there is global synchrony in the network, as shown in Fig. 1(b2) and (d2), respectively. Adding an autapse at node 2 (or 3) will not result in global synchronization (c2). These results indicate that the presence of a single autapse at a proper location in the network can promote synchronization.
Intuitively, from the network synchronization theory [33][34][35][36][37][38][39][40][41][42][43][44][45][46][47] , the introduction of an autapse is likely to discourage global synchronization because the time delay associated with the autapse effectively makes the dynamics of the neuron hosting the self-loop non-identical to other neurons in the network at which no autapses are present. Previous simulation of an isolated neuron showed 11 that, with the addition of an autapse, the neuronal dynamics could be changed drastically, e.g., from chaotic to periodic. If the dynamics of the individual nodes are characteristically different, the synchronization manifold existing for identical nodal dynamics may be destroyed 48 . The phenomenon of autapse promoted synchronization, as exemplified in Fig. 1(b2) and (d2), thus seems to be counterintuitive. However, an examination of the dynamical behavior of neuron 1 (or neuron 4) reveals that, even in the presence of an autapse, its trajectory is identical to that without the autapse. This observation suggests that, in spite of the autapse, the network still possesses a global synchronization manifold. In Methods, we provide a heuristic argument for the existence of a synchronization manifold in the presence of an autapse. (It is worth noting that, depending on the time-delay parameter, the dynamics within the synchronization manifold may differ from those of the individual, isolated neurons. For the time-delay parameter used in Fig. 1(b), there are periodic dynamics in the synchronization manifold).
To demonstrate that autapses can promote synchronization in a more quantitative manner, we exploit the concept of master stability function (MSF) 49 . Given a networked system possesses the synchronization manifold, the MSF method proposed in 1998 by Pecora and Carroll 49 has become the golden standard in the study of synchronization in complex dynamical systems. To appreciate the physical significance of the MSF method, we note that, as the network approaches a global synchronization state, the variables of the nodes approach each other asymptotically due to the mutual entrainment caused by the coupling. The subspace in which the synchronous solution lies is the synchronization manifold whose dimension is typically much smaller than that of the full phase space. Specifically, let x s (t) be the synchronous state of the network and {δx i (t)} be the infinitesimal perturbations added on node i at time t. Whether the perturbed system is restorable to the synchronous state can be inferred from the MSF 49 .
s s with δy and σ = ελ being the perturbation in the eigenmode space and the generic coupling strength, respectively. The quantity ε represents the uniform coupling strength, and {λ} are the eigenvalue spectrum of the Laplacian matrix defining the network structure (see Methods). Denoting Λ as the largest Lyapunov exponent calculated from the above equation, we have that the variation of Λ with σ gives the shape of the MSF curve. When the coupling configuration and the synchronous dynamics (i.e., the dynamics of an isolated node) are given, the MSF curve is independent of the particular network structure. If the dynamics in the synchronization manifold is stable with respect to perturbations in the transverse subspace, synchronization can be physically realized. In this case, the value of Λ is negative for all the transverse perturbation modes σ i = ελ i for i = 2, …, N (N is the number of coupled oscillators and the mode associated with λ 1 = 1 describes the motion along the synchronization manifold. See Methods for details). For the general nodal dynamics and coupling configuration, the value of Λ is negative only in certain interval 36,49,50 of the generic parameter. As such, for a network of coupled nonlinear oscillators to be synchronizable, the necessary condition is that all the generalized coupling strengths of the transverse modes {σ i } fall entirely the interval of negative Λ. (A mathematical formulation of the MSF is given in Methods). We calculate the MSF curve for the small neuronal network in Fig. 1 without or with an autapse at each node, to assess the impact of autapse on synchronization. For a fixed value of ε, the argument of σ (i.e., the generalized coupling parameter) is directly proportional to the eigenvalue λ of the network coupling matrix. For convenience, we choose ε = 1 and plot Λ versus λ (see Methods). Figure 2(a) shows such a plot for τ = 4 (the same parameter values as in Fig. 1). We have Λ < 0 for λ ∈ (σ 1 , σ 2 ) with σ 1 ≈ −0.7 and σ 2 ≈ 1. The structure of the small neuronal network gives λ 1 = 1 = σ 2 and λ 1 > λ 2 ≥ λ 3 ≥ λ 4 . Stable synchronization of the network is thus determined solely by the the smallest eigenvalue λ N : synchronization is (is not) achievable for λ N > σ 1 (λ N < σ 1 ). Without any autapse Fig. 1(a1), we have λ 4 = −0.73. Since λ 4 < σ 1 , the network is not synchronizable -in agreement with the numerical result in Fig. 1(a2). When an autapse is present at neuron 1 Fig. 1(b1), we have λ 4 = −0.5 > σ 1 , so the network is synchronizable, as indicated by Fig. 1(b2). Similarly, for an autapse at neuron 4 Fig. 1(d1), we have λ 4 = −0.5 > σ 1 so there is synchronization Fig. 1(d2). However, when the single autapse is attached to neuron 2 Fig. 1(c1) we have λ 4 = −0.72 < σ 1 , so the network is non-synchronizable, as indicated by Fig. 1(c2).
In general, when the network structure is given so that all eigenvalues {λ i } (i = 1, …, N) are known, whether stable synchronization can arise is determined by the MSF curve, where the range of the negative values of the MSF, namely the stable region, determines whether all the transverse eigenmodes are stable. For network systems without time delay (the general model studied in literature), the existence and the size of the stable region depend only on the local nodal dynamics and the coupling function, but not on the coupling strength 49,50 . For our neuronal network with autapses, however, the stable region depends also on the coupling strength parameter and the time delay. To uncover the effect of these two parameters on the stable region of the MSF curve, we decrease the coupling strength to ε = 0.8 and calculate the corresponding MSF, as shown in Fig. 2(b). Comparing Fig. 2(a) and (b), we see that not only the range but also the shape of the stable synchronization region are modified: the stable region now consists of three separated subregions: (−0.82, −0.75), (−0.72, 0.36), and (0.61, 0.67). To obtain a full picture of the distribution of the stable regions in the three-dimensional parameter space (ε, λ, τ), we show in Fig. 2(c-f) the value of Λ in four different parameter planes. In general, network synchrony is sensitive to the network structure (characterized by the eigenvalue spectrum {λ i }), the coupling parameter (ε), and the time delay parameter τ associated with the autapse. The sensitive dependence of stable synchronization region on variations in ε and τ is a typical feature of ragged synchronization 51 . The remarkable phenomenon is that, regardless of the changes in the parameters, insofar as there is an autapse, stable synchronization regions [blue regions in Fig. 2(c-f) persist. That is, for an unsynchronized neuronal networked system, the presence of a single autapse is able to induce stable synchronization.  Fig. 1, the stability of the MSF, as characterized by Λ, for different combinations of the system parameters. The synchronization state is stable when Λ < 0 for all the transverse modes, and will be unstable when Λ < 0 for any of the transverse modes. All results are obtained by solving Eq. (9) numerically (see Methods). (a) For ε = 1 and τ = 4, Λ versus λ. We see that Λ is negative in the region λ ∈ (−0.7, 1), indicating that stable synchronization will be achieved when all the transverse modes are located in this parameter interval. (b) For ε = 0.8 and τ = 4, Λ versus λ, where now Λ is negative in three subregions: λ ∈ (−0.82, −0.75), λ ∈ (−0.72, 0.36), and λ ∈ (0.61, 0.67). (c) For fixed ε = 1, contour plots of Λ in the parameter plane (λ, τ), where the blue color indicates parameter regions of stable synchronization (Λ < 0) and the red color specifies the regions where synchronization is unstable or cannot be realized physically (Λ > 0). (d-f) Contour plot of Λ in the parameter planes (λ, ε), (λ, τ), and (λ, ε) for fixed τ = 4, ε = 0.8, and τ = 6, respectively. The general phenomenon is that, through variations of the parameters (λ, τ, and ε), the region of stable synchronization can be readily modified. Autapse centralities and synchronization optimization. Suppose a small set of autapses are distributed into a neuronal network, to which set of neurons should the autapses be attached so that the network achieves the maximum synchronizability? To address this question for the general case of multiple autapses is difficult. We thus first consider the relatively simple case of a single autapse and develop a theoretical criterion to determine the optimal location for the autapse. Equivalently, it is necessary to find the autapse centrality, a quantity that measures the impact of a specific autapse on the network synchronizability. From the MSF formalism, we have that the autapse induces a shift in the eigenvalues of the network coupling matrix. Specifically, let G 0 and G be the network coupling matrices without and with the autapse, respectively. Letting λ j 0 and λ j be the jth eigenvalue of the respective matrices, we write the eigenvalue shift as λ . In the typical case where the MSF curve possesses only one bounded stable region, the network can reach a global synchronous state only if the eigenmodes in the transverse subspace associated with the largest and the smallest nontrivial eigenvalues (λ 2 and λ N , respectively) are both stable. The impact of the autapse on the network synchronizability can then be characterized by the shifts in the two nontrivial eigenvalues: , with i denoting the neuron that receives the autapse. A large value of λ Δ i 2 or λ Δ N i signifies a more significant impact. We thus propose the quantities λ Δ i 2 and λ Δ N i as two autapse centrality measures. To characterize the autapse centralities quantitatively, we develop a perturbation-based approach by treating the autapse as a small alteration to the network structure. With the perturbation, the normalized network matrix can be written as G = G (0) + G (1) , where G (1) is the perturbation matrix induced by the autapse. Assuming that the autapse is attached to neuron i which has degree k i , we have the elements of G (1) as if there is a link between neurons i and j, and For a densely connected network, we have  k 1 i and so  g g ij ii . Under the approximation , we have that the matrix G (1) has only one non-zero element: g ii (1) , leading to the following simplified version of Eq. (1): The smallest eigenvalue λ N can be treated in a similar way. We have To verify Eqs (2) and (3) numerically, we present in Fig. 3 the two autapse centrality measures for complex networks of distinct topologies. In particular, Fig. 3(a-f) are for a scale-free, random, and small-world network, respectively. The three networks have the same size N = 100 and average degree 〈k〉 = 6. Introducing a single autapse onto each neuron alternatively, we calculate the two centrality measures: the variations of the two extreme eigenvalues, Δλ 2 and Δλ N . Figure 3(a,b) show the autapse centralities versus the nodal index and degree, respectively. From Fig. 3(a), we see that the centralities exhibit large variations among the neurons, while Fig. 3(b) shows that the hub neurons are less sensitive to the autapse than the non-hub neurons in spite of the large fluctuations in the values of Δλ 2 and especially of Δλ N . Figure 1 Fig. 3(e). The results for the small-world network are similar to those for the scale-free network Fig. 3(f). Overall, there is a reasonable agreement between the two centrality measures calculated numerically and predicted theoretically.
We now address the issue of multiple autapses. By introducing autapses to a judiciously chosen set of neurons (one autapse to each neuron), we seek to optimize network synchronization. We again exploit Eqs (2) and (3), which suggest that the optimal set of neurons should be selected according to the largest possible values of the autapse centralities calculated based on a single autapse. Suppose the number of available autapses is m. We arrange the centrality values for individual nodes in a descending order and choose the m neurons from the top of the list. To test this optimization strategy, we compare its performance with those of two alternative strategies: degree-based and random, where for the former, the set of neurons (each receiving one autapse) are selected according to the nodal degree in the descending order, while for the latter, the set is chosen randomly. In simulations, we use the scale-free network in Fig. 3(a), with the same nodal dynamics and coupling function as in Fig. 2(a). As the stable region of the MSF curve is bounded Fig. 2(a), the relevant autapse centrality is Δλ N : the set of neurons to the receive the autapses can be chosen according to the descending order of the values of Δλ N . We calculate the network-averaged synchronization error of the variable x: number of autapses m is increased through a critical point (denoted as m c ), the synchronization error decreases to zero, as shown in Fig. 4. For our autapse centrality based optimization strategy, we have m c ≈ 45, while the degree based and random strategies give m c ≈ 52 and m c ≈ 60, respectively. Among the three strategies, ours yields the minimum number of autapses required for the network to be fully synchronized.
Generality and robustness of autapse-centrality based synchronization strategy. We address a number of issues to further substantiate the main result that autapses could promote global synchronization in neuronal networks, including: (1) the performance of autapse-centrality based synchronization strategy for other types of complex networks; (2) the impacts of parameter mismatches among neurons on global synchronization; and (3) the influences of the system parameters on the minimum number of autapses required for global synchronization.
So far we have demonstrated that the autapse-centrality based strategy outperforms the conventional degree based and random strategies in promoting global synchronization in scale-free network, as shown in Fig. 4. Our theoretical analysis implies that this result should hold for the general network model, regardless of the network topologies. An example is given in Fig. 5, which displays the synchronization performance of the three strategies for the random and small-world networks. We see that, as for the case of scale-free networks in Fig. 4, the autapse-centrality based strategy always requires the least number autapses in achieving global network synchronization.
In real systems, parameter mismatches among the neuronal dynamics and the network links are inevitable. To gain insights into whether autapses could promote synchronization in networks consisting of non-identical neurons, we consider the toy network model in Fig. 1 and introduce small mismatch in the parameters (I, ε, and τ).
To be concrete, we vary the coupling strength w aut and the time delay τ aut of the autapse and investigate how network synchronization is deteriorated as the values of w aut and τ aut are deviated from those of the synaptic connections. Setting w = 1 and τ = 4 for the synaptic connections, we calculate the synchronization error 〈δx〉 as a function of w aut . The numerical results are shown in Fig. 6(a). We observe 〈δx〉 = 0 for w aut > 0.3. To see how mismatched coupling strength affects synchronization, we plot in Fig. 6(b) and (c) the neuron trajectories for w aut = 0.1 and w aut = 0.5, respectively, where there is still global synchronization for the latter but for the former, synchronization is lost. This indicates that, even when the mismatch of w aut is about 50% Fig. 6(c), synchronization can still be achieved. In fact, given w aut > 0.3, global synchronization will be maintained regardless of the mismatch amplitude Fig. 6(a). The underlying reason for the robustness of global synchronization to the mismatch of w aut is that, due to the normalized coupling scheme (see Methods), the dynamics of the synchronous manifold is independent of the autapse strength. Figure 6(d) demonstrates the robustness of network synchronization in the presence of mismatch in time delay, where the behavior of 〈δx〉 versus Δτ = τ aut − τ aut is shown for fixed w aut = 1. We see that the network is well synchronized for a small amount of mismatch in the time delay. For example, we have 〈δx〉 < 0.1 for Δτ ∈ (−0.08, 0.1). For a relatively large amount of mismatch, synchronization is lost, as shown in Fig. 6(e) for Δτ = 0.2. In general, for Δτ ≠ 0, the local dynamics of the autapsed neuron are different from those of the other neurons, making global synchronization difficult. However, as shown in Fig. 6(f), if the mismatch amount Δτ is not too large, all neurons in the network can still be well synchronized.
In general, the number of autapses required to synchronize a complex network depends on the network parameters such as the connecting density, the degree distribution, the delay parameter, and the coupling strength among the neurons. For the scale-free network in Fig. 4 of 100 nodes and average degree 〈k〉 = 6 and by the parameters ε = 1 and τ = 4, m c = 45 autapses are required for global synchronization. This critical number of autapses, however, could be changed significantly by varying the network parameters. To shown an example, we fix the network structure and the parameter τ used in Fig. 4 but changing the value of ε to 0.95, and recalculate the number of required autapses for different synchronization strategies. The numerical results are shown in Fig. 7. We see that, for the autapse-centrality based, degree based, and random strategies, the critical numbers of autapses required for global synchronization are reduced to, respectively, to m c = 19, 31, and 35.

Discussion
In this paper, we study neuronal autapses, a type of structural anomalies that give certain neurons a time-delayed, self-feedback loop. Autapses, when they were first discovered 1 , were thought of playing no specific biological or physiological role. Only later were their potential impacts in biological functions found or suggested [3][4][5][6][7][8] . In view of the fundamental importance of synchrony in biology, physiology, and biomedical sciences, we focus our study on the possible role of autapses in modulating neuronal synchrony through computation and theory. Our main finding is that, for a complex neuronal network, even the existence of autapses on a small fraction of the neurons can promote synchrony significantly. In particular, using both a small toy network and larger networks of distinct complex topologies, we introduce autapse to a fraction of neurons and analyze the change in the network synchronizability. Based on a systematic analysis of the impact of a single autapse at different nodes in the network, we develop centrality measures to analyze the situation of multiple autapses in terms of their optimal placement in the network to realize robust synchronization of the entire network. Our work provides the computational and theoretical foundation for hypothesizing the positive role of autapses in promoting synchronization in neuronal networks.
It is an accepted notion that synchrony in neuronal networks is strongly correlated with certain neurological diseases such as epileptic seizures 18 , but whether such a disease can be attributed to an increasing level of synchrony has been a controversial issue [19][20][21][22] . At the present it remains difficult to ascertain, through the traditional approach of EEG or ECoG data analysis, whether a general, well defined causal relation between an elevation of synchrony and the occurrence of seizure exists. The main idea underlying our work is to focus on "unusual" features or structural anomalies in neurons and to study their role in modulating neuronal synchrony. Suppose for certain type of seizure, the responsible neuronal network in a brain region can be identified. Whether there exist structural anomalies in some neurons in this region can then lead to insights into the interplay between synchrony and the particular type of seizure. Generally, the issue of synchronization goes much beyond the scope of epileptic seizures with broad implications to biology and physiology 16 .
We discuss a few technical issues. Firstly, to make the model theoretically tractable, we assume that the isolated neuronal dynamics are identical and the links (including the autapses) have the same coupling strength and time delay. However, in real systems, parameter mismatches among the neuronal dynamics and network links can be expected. We have gained preliminary insights into whether autapses can promote synchronization in networks of non-identical neurons, by introducing small mismatches in parameters (I, ε, and τ) in the toy network model in Fig. 1. Our computations reveal, for example, when the time delay parameter has a 2% mismatch among the network links, global synchronization can still be achieved. Secondly, in our study we adopt the scheme of normalized couplings to enable computation and analysis of the MSF and the proposal of autapse centralities. If the couplings are not normalized, such a theoretical analysis would not be possible and it appears at the present that one must rely on numerical simulations to uncover and quantify the role of autapses in synchronization. Thirdly, we assume linear feedback couplings in the network. While this type of coupling scheme does have physiological support (e.g., in terms of electrical synapses underlying the neuronal connections), it remains unknown whether other types of couplings, e.g., those with chemical synapses, can lead to effective modulation of network synchronization by autapses.
Finally, we point out the significance of our work in general network science and engineering with respect to the problems of synchronization and control. A focus in this area of research has been on the role of long-range connections, such as those in small-world networks 53 , in promoting network synchronization 36 . From a practical standpoint, long range connections may be costly 54 . For example, for neural networks in the human brain 55 , while there are white fiber tracts among distant cortical areas which are essential for the brain functions, the network connectivity is still dominated by short-distance, local connections, due to the minimization principle of the axonal length and energy consumption. A similar issue arises in other realistic systems 54 such as traffic networks, power grids, communication network, and the Internet. Our finding that autapses, the shortest possible links in any network, can play a positive role in promoting synchronization, is encouraging from the perspective of optimal control of complex network dynamics at a minimal cost.

Methods
Model description. Our model of neuronal network with autapses reads where i, j = 1, 2, …, N are nodal indices, N is the network size, x i the state vector of the i th neuron, ε is the coupling parameter, F i (x) is the velocity field governing the dynamics of the i th neuron when it is isolated, H(x) is the coupling function, τ ij denotes the time delay of signal propagating from neuron j to i, and η i (t) represents the independent, identically distributed Gaussian white noise added to the membrane potential of each neuron, with 〈η i (t)〉 = 0 and 〈η i (t)η j (t)〉 = (2 × 10 −10 )δ ij . The network structure in the presence of autaptic connections is characterized by the coupling matrix C = {c ij }. For non-diagonal elements, we have c ij = c ji = 1 if there is a link between neurons i and j, otherwise c ij = 0. For the diagonal elements, we set c ii = δ il , where δ il is the Kronecker delta function and V = {l} (with l = 1, …, m) denotes the set of neurons with autaptic connections. Specifically, we have c ii = 1 if neuron i possesses an autapse, and c ii = 0 otherwise. The quantity k i = ∑ j ≠ i c ij is the number of links connected to neuron i (i.e., the degree). To capture the essential network dynamics subject to autapses while making the model theoretically tractable, we assume identical time delays: τ ij = τ. Biologically, this approximation is reasonable for synapses within the same cortical minicolumn 55 . The scheme of diffusive (linear feedback) coupling assumed in Eq. (4) models the realistic electrical synaptic interactions (e.g., of the gap-junction type) among neurons 16,56 .
In simulations, we set the parameters of the HR oscillator as (a, b, c, d, r, s, x R , I) = (1, 3, 1, 5, 6 × 10 −3 , −1.6, 3.2), for which the isolated HR neuron exhibits chaotic bursting dynamics 28 . The coupling function is chosen to be H(x) = [x, 0, 0] T , i.e., the neurons in the network are coupled through their membrane potentials. We use the Bogacki-Shampine algorithm 57 to simulate Eq. (4), which is a time-delayed, stochastic system of coupled differential equations. The integration time step is δt = 1 × 10 −3 . Persistence of a global synchronization manifold in the presence of an autapse. We argue that, given that the autapses have an identical time delay with respect to the synaptic connections, the neuronal network possesses a synchronization manifold. Consider the network in Fig. 1(b1) where the 1st and 2nd terms on the right side represent, respectively, the local dynamics of the autapsed neuron and the coupling signals it received from the neighboring neurons. Assume that the network is globally synchronized, and denote x s as the dynamical state within the synchronous manifold, we have x i (t − τ) = x s (t − τ) (for i = 1, …, 4). Equation (5) can then be rewritten as For neurons without autapse (j = 2, 3, 4), their dynamics are governed by which is identical to Eq. (5). That is, in the synchronous state all neurons in the network, regardless of whether they are autapsed or not, follow the same dynamical equation. Further, because of the normalized coupling scheme, the equations are independent of the network structure. A global synchronization manifold thus persists in the presence of autapses, regardless of their locations in the network.