Gap junctions desynchronize a neural circuit to stabilize insect flight

Insect asynchronous flight is one of the most prevalent forms of animal locomotion used by more than 600,000 species. Despite profound insights into the motor patterns1, biomechanics2,3 and aerodynamics underlying asynchronous flight4,5, the architecture and function of the central-pattern-generating (CPG) neural network remain unclear. Here, on the basis of an experiment–theory approach including electrophysiology, optophysiology, Drosophila genetics and mathematical modelling, we identify a miniaturized circuit solution with unexpected properties. The CPG network consists of motoneurons interconnected by electrical synapses that, in contrast to doctrine, produce network activity splayed out in time instead of synchronized across neurons. Experimental and mathematical evidence support a generic mechanism for network desynchronization that relies on weak electrical synapses and specific excitability dynamics of the coupled neurons. In small networks, electrical synapses can synchronize or desynchronize network activity, depending on the neuron-intrinsic dynamics and ion channel composition. In the asynchronous flight CPG, this mechanism translates unpatterned premotor input into stereotyped neuronal firing with fixed sequences of cell activation that ensure stable wingbeat power and, as we show, is conserved across multiple species. Our findings prove a wider functional versatility of electrical synapses in the dynamic control of neural circuits and highlight the relevance of detecting electrical synapses in connectomics.

Insect asynchronous flight is one of the most prevalent forms of animal locomotion used by more than 600,000 species. Despite profound insights into the motor patterns 1 , biomechanics 2,3 and aerodynamics underlying asynchronous flight 4,5 , the architecture and function of the central-patterngenerating (CPG) neural network remain unclear. Here, on the basis of an experiment-theory approach including electrophysiology, optophysiology, Drosophila genetics and mathematical modelling, we identify a miniaturized circuit solution with unexpected properties. The CPG network consists of motoneurons interconnected by electrical synapses that, in contrast to doctrine, produce network activity splayed out in time instead of synchronized across neurons. Experimental and mathematical evidence support a generic mechanism for network desynchronization that relies on weak electrical synapses and specific excitability dynamics of the coupled neurons. In small networks, electrical synapses can synchronize or desynchronize network activity, depending on the neuron-intrinsic dynamics and ion channel composition. In the asynchronous flight CPG, this mechanism translates unpatterned premotor input into stereotyped neuronal firing with fixed sequences of cell activation that ensure stable wingbeat power and, as we show, is conserved across multiple species. Our findings prove a wider functional versatility of electrical synapses in the dynamic control of neural circuits and highlight the relevance of detecting electrical synapses in connectomics.
With over a million known species, insects comprise the largest group of animals on Earth 6 . Their considerable evolutionary success has been attributed to small body size and the ability to fly. These two features provide access to unutilized niches and rapid translocation, but aerodynamic constraints in small flyers require high wingbeat frequencies, and space constraints demand miniaturization of the central nervous controllers for flight 7 . In 75% of all flying insect species, highly specialized, indirect, asynchronous flight muscles form an oscillatory system that generates wingbeat frequencies of 100-1,000 Hz by reciprocal stretch activation of antagonistic wing muscles to ensure forward propulsion at low Reynolds numbers 1,8 . The flight motoneurons (MNs) that innervate asynchronous flight muscles fire at much lower frequencies, therefore not activating the muscles on a cycle-to-cycle basis. Nonetheless, power output is regulated by a CPG network in the central nervous system that controls MN firing frequencies to adjust the myoplasmic calcium levels that, in turn, regulate wingbeat frequency and amplitude 1 . Although asynchronous flight has emerged independently 7-10 times during evolution 8 , neither the principles of CPG architecture for generating MN output from the miniaturized central nervous system of asynchronous flyers nor the functional consequences thereof have been identified.

Asynchronous flight motor patterns
To quantify asynchronous flight patterns and decipher CPG architecture, we used the firing output of the five identified MNs (MN1-5) innervating the dorsal longitudinal wing depressor muscle (DLM) of the genetic model system [9][10][11] Drosophila melanogaster as well as other insect species (to test for generality). The DLM provides the force for wing downstroke, consists of six muscle fibres, each of which is innervated by one of five identified MNs 9-11 (MN1-5; Fig. 1a). MN1-4 each target one of the four most ventral DLM fibres ipsilateral to their somata, whereas MN5 innervates DLM fibres 5 and 6 on the side contralateral to the MN5 soma (Fig. 1a). This neuromuscular architecture is conserved across insect species examined (locust 12 , moth 13 , blowfly 14 ).
In vivo recordings of MN1-5 from their DLM target muscle fibres using simultaneous laser-based wingbeat detection during tethered flight show that each MN fires only every approximately 20th to 40th wingbeat (Fig. 1a-c). Although the firing frequencies of MN1-5 can vary between animals and are adjusted on demand 1,15,16 , within a given animal and power demand, all five MNs always fire at the same frequencies (Fig. 1b) and with similar variance in the interspike Nature | Vol 618 | 1 June 2023 | 119 interval (ISI; coefficient of variation is 0.22 ± 0.07 (mean ± s.d.) for each MN). MN1-5 firing frequencies directly control muscular tension and stretch activatability by adjusting myoplasmic calcium levels and, therefore, wingbeat frequency and stroke amplitude 1,17 . Accordingly, in single animals, alterations in power demands go along with changes in MN1-5 firing frequency that are proportional to wingbeat frequency changes ( Fig. 1c (grey dots)). A linear relationship between average MN firing and average wingbeat frequency within the normal flight working range was further confirmed by recordings in 100 animals, although interanimal variability yields a lower correlation between MN firing rate and wingbeat frequency (Fig. 1c (red dots)). This correlation is increased when analysing the changes in the MN firing rate in relation to changes in the wingbeat frequency (Extended Data Fig. 4e-g). Thus, the central nervous system controls asynchronous flight muscle power output neither by the recruitment of different motor units, nor on the scale of single wingbeats, but the frequency of MN1-5 population firing is the key regulator of wing power production 1 . We show that MN excitability is well suited to dynamically regulate wing power output. First, MN1-5 respond to constant input with slow tonic firing (Fig. 1d). Second, a nearly linear input-output relationship as observed in f-I curves (3-30 Hz; Fig. 1e) covers and even exceeds the working range observed during tethered flight (approximately 3-12 Hz; Fig. 1c). The excitability of MN1-5 is therefore tuned to linearly translate synaptic input into tonic firing to regulate wingbeat in the working range of flight.
Despite equal firing frequencies at any given wingbeat frequency 1,15 (Fig. 1a,b), the five MNs in this system do not fire in synchrony. Instead, as previously suggested 18 , their spikes are splayed-out in time with firing phases across neurons dispersed approximately equidistantly, resulting in stereotyped preferred sequences of MN1-5 spiking (Fig. 1g), which we name splay states. A network splay state means that each individual cell in a network of N neurons fires regularly, yet with a constant, non-zero phase relationship with respect to the firing of each other neuron. Firing phase differences between neurons approximately correspond to 1/Nth of a neuron's ISI or a multiple (<N) thereof. Importantly, the sequence of MN1-5 firing can change intermittently, but the CPG robustly slides back into one or two of the most preferred splay states. Notably, the same splay states are preferred across animals (Fig. 1f), indicating that they result from hard-wired CPG circuitry. Splayed-out firing results in characteristic phase relationships between each pair of MN1-5. These phase relationships are conserved across individuals (Extended Data Fig. 1a-c). The phase relationships of all MN1-5 pairs show a gap around phase zero, indicating out-of-phase firing. Given that multiple splay states exist during flight (Fig. 1f,g), firing of MN pairs is not precisely phase locked (Extended Data Fig. 1a-c). Yet, restricting phase histograms to Data are mean ± s.d. c, MN firing frequency and wingbeat (WB) frequency (the red bars indicate working ranges) are linearly related within an animal (grey dots; correlation coefficient, r 2 = 0.63; P < 0.0001, two-tailed t-test) and with larger variance also across animals (red dots; n = 100; correlation coefficient, r 2 = 0.31; P < 0.0001, two-tailed t-test). d, The firing responses of MNs (top traces) to current injections of different amplitudes (bottom traces). e, The mean MN response frequency ( f ) and injected current amplitude (I) are approximately linearly related for 2-30 Hz (n = 15 animals), therefore exceeding the normal MN firing frequencies observed during flight (inset; ~3-12 Hz, data from c). Data are mean ± s.e.m. (main plot) and median ± range (inset). f, During flight, MN1-4 spikes are dispersed in time (in splay state) with characteristic sequences. Each animal switches between different splay states during flight, but the same splay states are preferred across individuals (n = 8). The box plots show the median (centre line), quartiles (box limits) and range (error bars). g, Timing of MN1-4 spikes in four subsequent splay states (1423 (red); 1243 (turquoise); 1234 (green); 1324 (dark blue)) during flight. Article one splay state narrows the phase relations (Extended Data Fig. 1d). Moreover, synchronous spikes between MN pairs occasionally occur during switches between different splay states, but not within a stable splay state (Extended Data Fig. 1d).
Our characterization of phase relationships between pairwise combinations of the MN1-5 enables us to test for across-species conservation of the motor patterns (see Supplementary Videos 1-3 for the flight of three species). This reveals a notable similarity between D. melanogaster and the gold fly Lucilia sp. (Extended Data Fig. 2a,b). Similar phase relationships have also been suggested for Calliphora erythrocephala 14 , Eucaliphora lilaea and Musca domestica. Moreover, using the MN4-MN5 pair (Extended Data Fig. 2c), our measurements indicate conservation of CPG architecture between two Drosophila species (D. melanogaster and Drosophila hydei), other dipteran genera (Calliphora and Musca) and, to a certain degree, even between Diptera (for example, flies) and Hymenoptera (for example, honey bee). Given that asynchronous flight evolved multiple times independently 8 , splayed-out firing has probably provided selective benefits over millions of years. To gain mechanistic and functional insights into the CPG network that controls asynchronous flight, we next analysed the network principles underlying splay state motor pattern generation and the resulting functional benefits for flight performance.

Network splay states produced by a minimal CPG
From invertebrates 19 to mammals 20 , the timing and pattern of MN activation during locomotion typically rely on networks of premotor interneurons. By contrast, we provide evidence that splayed-out MN1-5 firing is generated by interactions between the MNs themselves. First, unpatterned optogenetic activation of excitatory, cholinergic input to MN1-5 increases MN firing frequencies without changing their phase relationships, preserving a characteristic gap around phase 0 (Extended Data Fig. 3). We confirmed both physiologically 21 and anatomically that MNs receive excitatory cholinergic input to their dendrites (Extended Data Fig. 4a,b), and that all MNs share common synaptic input (Extended Data Fig. 4c,d) as previously proposed 18 . Tonic common excitatory synaptic drive to MN1-5, as induced by optogenetic stimulation, is therefore transformed into out-of-phase MN1-5 firing. Second, unpatterned, tonic, optogenetic stimulation of the five MNs (see Extended Data Fig. 5a,b for selective expression of transgenes in MN1-5) during flight increases their firing frequency and, therefore, wingbeat power, but it does not change phase relationships between MNs, as exemplified by the same characteristic gap around phase 0 for the MN4-MN5 pair (Extended Data Fig. 5). The strength of this dip around phase 0 in MN phase histograms is strongly correlated with full network splayness (definitions are provided in the Methods) in both simulated and experimental data (Extended Data Fig. 6). Similar phase histograms of MN pairs before and during optogenetic stimulation therefore indicate unaltered splayness. The generation of continuous splayed-out firing after selective, unpatterned, optogenetic stimulation of all five MNs indicates that interactions between MNs shape the firing patterns. An alternative possibility would be synaptic output from MN1-5 to feedback interneurons, which then generate the pattern. However, neither trans-synaptic mapping nor genetic markers for central output synapses from MN1-5 provide evidence for postsynaptic partners in the ventral nerve cord (VNC; Extended Data Fig. 7). Thus, it seems probable that the timing and pattern of MN activation does not require patterned activity of interneurons and, instead, the MN1-5 ensemble constitutes a minimal CPG.
One possible mechanism to transform common excitatory input to an ensemble of MNs into dispersed firing is lateral inhibition among MN1-5 by chemical synapses 15 . We reject this possibility by combining targeted genetic manipulation of MNs with in vivo recordings during flight. Knockdown (KD) of receptors for inhibitory transmitter (GABA A receptors (GABA-ARs) and glutamate-gated chloride channel (GluCl)) increases MN firing frequencies, but the phase relationships remain similar (Extended Data Fig. 6e,f). Lateral inhibition through chemical synapses is therefore not required for pattern generation.
Another possibility to create a neural network exclusively from MN1-5 is to connect them with electrical synapses. This has previously been suggested 22,23,24 , but experimental evidence has been lacking, and electrical coupling seems difficult to reconcile with desynchronized firing. We have tested this by genetic manipulation of innexins 25 , the invertebrate counterparts of connexins 26 , which comprise the pore-forming proteins of electrical synapses. ShakingB (ShakB) is the innexin expressed in the Drosophila escape circuit, including in MN1-5. KD of shakB using RNA interference (RNAi) in MN1-5 reduces electrical coupling below the detection threshold in paired patch-clamp recordings (Fig. 2g).
Genetic manipulation of shakB in the five MNs disrupts the splay state without affecting the firing frequencies (Fig. 2a,b). As exemplified for the MN4-MN5 pair, in control animals (n = 7), firing of MN5 is inhibited before and after the occurrence of MN4 spikes (Figs. 2a,b (top)), therefore resulting in a characteristic gap around phase 0 (Extended Data Figs. 3b and 5e) that is indicative for network splayness (Extended Data Fig. 6a-d). By contrast, with shakB KD in MN1-5, the gap around phase 0 is absent and firing is more random with a slight preference for in-phase firing ( Fig. 2a,b (middle)). Thus, electrical synapses are required for firing desynchronization in this small CPG. Recordings of all five MNs during flight confirm that shakB KD in MNs impairs their normal phase relationships and, therefore, the splay state (compare Extended Data Fig. 8 with Extended Data Fig. 1). This contradicts the common notion that electrical synapses function to synchronize network activity 27 .
Key to this role of gap junctions is weak electrical coupling, because increasing electrical coupling by overexpression of ShakB in MN1-5 causes firing synchronization (Fig. 2a,b (bottom)). Thus, the in vivo data indicate that weak electrical synapses cause firing desynchronization, but stronger coupling synchronizes firing. Weak electrical coupling between MN1-5 is also supported by anatomical experiments. Although we have previously not observed diffusion of small dye tracer molecules through gap junctions between MN1-5 16,21 , in Drosophila neurons, dye uptake during iontophoresis and subsequent diffusion through shakB-encoded electrical synapses are increased by KD of the fragile-X mental retardation protein (FMRP) (encoded by dfmr1) 28 . In this genetic background dye coupling of MNs is reliably observed (Fig. 2c).
In summary, a minimal CPG of five electrically coupled MNs is sufficient to pattern splayed-out firing across power demands and does not rely on additional interneurons or chemical synapses. First, lateral inhibition by chemical synapses is not required (Extended Data Fig. 6e,f); second, we found no evidence for chemical output synapses from MN1-5 in the VNC (Extended Data Fig. 7); third, unpatterned optogenetic activation of either presynaptic cholinergic neurons or the MNs during tethered flight increases MN1-5 firing rates and wingbeat frequency (Extended Data Figs. 3 and 5), but the phase relationships between MNs remain unaltered (compare Fig. 2a with Extended Data Figs. 3b and 5e). Finally, weak electrical synapses between MNs are required to generate normal phase relationships (Fig. 2a,b). However, this raises the question of how electrical synapses cause firing desynchronization.

Electrical coupling within the CPG
A first step in deciphering the mechanism is to characterize the electrical synapses. Dual in situ patch-clamp recordings of MN pairs confirm electrical coupling and demonstrate that the electrical synapses are weak, bidirectional and non-rectifying (Fig. 2d). Non-rectifying, because depolarizing and hyperpolarizing current injections (20 ms duration) into one MN cause gap junctional potentials in the other MN. The relationship between presynaptic and postsynaptic charge is linear (Fig. 2e,f). Bidirectional, because the direction of charge transfer can be reversed (Fig. 2d). Compared with common CCs (postsynaptic charge/presynaptic charge 29 ) estimated in different types of neurons 29 (0.02-0.2), electrical synapses between the five MNs are weak, but coupling is twice as strong for the MN1-MN2 and MN3-MN4 pairs (coupling coefficients (CC) = 0.023 ± 0.003) compared with all of the other possible combinations of MN1-4 pairs (CC = 0.01 ± 0.0027; Fig. 2f).
The possibility that electrical synapses may desynchronize neural networks under specific conditions has been suggested by theoretical work 30-33 . Experimentally, a transient desynchronization of electrically coupled cerebellar Golgi cells has been described in one study, but for the specific condition of sparse input 34 . There, the transient desynchronization is attributed to a more effective gap junctional transmission of the afterhyperpolarization (AHP) of the action potential compared with its depolarizing overshoot 34 . Similarly, we found that the AHP of MN spikes is transmitted more effectively (CC, 0.042 ± 0.04) through ShakB-mediated electrical synapses than the brief spike overshoot (CC, 0.01 ± 0.004; Fig. 2h). CCs can differ for different components of the action potential because the duration of the presynaptic signal and the time constant of the postsynaptic membrane shape the junction potential 29 . Paired in situ current-clamp recordings of MNs that were induced to fire tonically by somatic current injection indicate that firing of one MN can depress firing of the other during and shortly after the AHP (Fig. 2i). However, the preconditions for desynchronization, as described for electrically coupled cerebellar Golgi cells 34 , are probably not fulfilled in the insect asynchronous flight CPG. First, network desynchronization by electrical synapses during flight is not transient but permanent. Second, network desynchronization manifests not only under sparse synaptic input regimes, but through the full range of synaptic input that occurs during flight. Thus, although we cannot rule out a contribution of the AHP, a general mechanism for small network desynchronization without a pronounced AHP (which is probably not sufficient; see above) or inhibitory synapses (not necessary for MN desynchronization; Extended Data Fig. 6e,f) is required. Indeed, theoretical considerations of network connectivity as well as cell-intrinsic excitability provide insights into how the observed splay states are generated.

Coupling strength and excitability class
According to the theory of coupled phase oscillators, splayed-out network states are obtained if pairs of neurons have a preference to fire out of phase, that is, they are phase-repulsive 24 . However, for networks of more than two neurons, a strict antiphase state that maximizes the difference between the firing phases of two neurons cannot be achieved across all neuronal pairs; these systems are therefore called frustrated 35,36 . Typically, they can organize into a splay state that equidistantly maximizes phase distances among neurons to minimize frustration. The phase difference between individual pairs of neurons is therefore not maximal (that is, smaller than half the ISI), but firing phases of neuron pairs are still separated across the population of cells, which reflects the splay state that we observe in vivo (Fig. 1a,f,g and Extended Data Figs. 1 and 6a-d).
Whether pairs of neurons tend to fire out of phase depends on the specific combination of synaptic connectivity and cellular excitability (in other words, the dynamics of action potential generation). Cellular excitability classes characterize qualitatively different types of spiking dynamics and are associated with distinct mathematical bifurcation types at spike onset (that is, threshold). Given the weak electrical coupling among MN1-5, in theory, repulsive phase coupling 24  Article to generate a gap junctional network of MN1-5 (Fig. 3a). Subjecting the single neuron model to a mathematical bifurcation analysis within the physiological parameter range identified a cellular excitability class that our analysis predicts to favour out-of-phase firing (through a homoclinic spike-onset near the saddle-node loop (SNL) bifurcation; Methods and Extended Data Fig. 9a). Indeed, five identically, weakly electrically coupled MNs with cellular dynamics near the SNL point robustly exhibit a desynchronized splay state (Fig. 3b). Note that the minimal conductance-based model does not contain a pronounced AHP ( Fig. 3b (spike shape)), showing that the presynaptic AHP is not mandatory for network desynchronization.
By contrast, systematic variation in electrical coupling strength reveals that weak electrical coupling is required because models with non-zero, small CCs (<0.05), as obtained in vivo (Fig. 2f), yield splayed-out firing (Fig. 3c), whereas, for CCs between 0.05 and 0.21, network in-phase synchronization increases and, for CCs of >0.21, the network state is synchronized (Fig. 3c). To test these model predictions experimentally, we manipulated gap-junction strength genetically and quantified the synchrony of the MN4-MN5 pair from in vivo recordings during flight (Fig. 3d). In control animals with weak gap junctions, the synchronization index (Methods) is low, similar to model simulations with weak gap junctions. RNAi KD of electrical synapses increases

The mechanistic core of the splay state
We next examined why cellular excitability dynamics close to the SNL point favour splay states. For regularly firing cells with all-or-none action potentials, three main dynamical excitability classes-shaped by cell-intrinsic properties including the ion channel compositionexist 40 . Mathematically, these correspond to three distinct spike-onset bifurcations: the subcritical Hopf, the saddle-node-on-an-invariant cycle (SNIC) and the saddle homoclinic orbit (HOM) bifurcations ( Fig. 3e (top row)), which in turn qualitatively determine the phase dependence of a neuron's sensitivity to inputs (also termed the phaseresponse curve (PRC)) in a bifurcation-specific and, therefore, excitability-class-specific manner ( Fig. 3e (middle row)). For the network state, the coupling function, which combines the influence of the PRC and electrical coupling, is decisive (definitions are provided in the Methods), therefore assigning a direct impact on network dynamics to both gap junctions and cellular voltage dynamics 38 . Specifically, an asymmetry of the PRC shape favours stable phase relationships of spiking in the network; these are predicted by the fixpoints of the coupling function's odd part (that is, fixpoints that can be read off the asymmetric components of the function; Methods). Our analysis shows that, for weak reciprocal electrical coupling, out-of-phase firing of identical neurons can be obtained only in combination with a PRC that is monotonically decreasing around phase 0.5 (a graphical explanation is provided in Extended Data Fig. 9b,c). Such a PRC is found for HOM excitability, including dynamics close to the border of SNIC firing, the latter corresponding to the SNL point introduced above (Fig. 3e (middle row)). Here, weak electrical coupling combined with such cellular dynamics results in a stable fixpoint at out-of-phase firing (Fig. 3e (bottom row)). By contrast, for the other two spike-onset bifurcations, synchronous in-phase firing is expected (compare also the stable fixpoints for the phases between neurons; Fig. 3e (bottom row)). Precisely, we predict that MN1-5 should dynamically reside close to the SNL point. Here, the PRC shape enables a network splay state, yet the slope of the neurons' firing-rate versus current-input (f-I) curves is not as steep as deeper in the HOM regime, therefore enabling smooth control of MN firing frequency and wingbeat power, as observed in vivo ( Fig. 1c-e). To test this model prediction, the excitability class of MN1-5 has to be manipulated in vivo. Shab potassium channels constitute around 50% of the delayed rectifier current in these MNs (Extended Data Fig. 10). Bifurcation analysis of the single-cell MN model reveals an SNL bifurcation (Extended Data Fig. 9a). Moreover, increases in total Shab channel conductance, g shab , induce several transitions in excitability class 37 , from HOM (at low g shab ), through the SNL point into a SNIC (medium g shab ) and a Hopf bifurcation (high g shab , Fig. 3e and Extended Data Fig. 9a). At both excitability dynamics resulting from higher g shab (SNIC and Hopf), the weakly coupled MN1-5 network exhibits a synchronized state, whereas SNL dynamics favour splay states (as described above). We tested this model prediction experimentally by genetic manipulation in vivo. Targeted overexpression of Shab in MN1-5 causes a near doubling of Shab delayed rectifier current (Fig. 3f) without compensatory changes in the amplitudes of non-Shab delayed rectifier current (Extended Data Fig. 10). Recordings of the MN3-MN4 pair were used to test whether this caused firing synchronization in vivo. In controls, both MNs fire out of phase ( Fig. 3g (top) and Extended Data Fig. 1). By contrast, with Shab overexpression, the MN3-MN4 pair shows markedly increased synchronization (Fig. 3g (bottom)). Variability in synchronization strength between animals (Fig. 3h) is probably caused by different overexpression strengths and variable levels of other delayed rectifier channels. We also cannot exclude some compensatory regulation of yet other currents. However, pooling the data from all recordings (n = 10) revealed a significant increase in firing synchrony after overexpression of Shab (Fig. 3h), as predicted by an increase in g shab in the network model (Fig. 3h). By contrast, our model predicts that a decrease in g shab shifts MNs deeper into the HOM regime but does not cause a transition in excitability class ( Fig. 3e and Extended Data Fig. 9a). Thus, with reduced g shab in MN1-5, the network model predicts desynchronization (Fig. 3e). We tested this prediction experimentally. Compared with the control (n = 7 animals) targeted RNAi KD of Shab in MN1-5 (n = 11 animals) reduced the Shab current by around 70% (Extended Data Fig. 10). Paired in vivo recordings of MN4 and MN5 during flight confirm a low synchronization index as predicted by the model (median synchronization indices are 0.559 for control flies (n = 11 animals) and 0.598 for Shab KD (n = 11) and were not statistically different; P = 0.393, U = 74, two sided Mann-Whitney U-test).
The theoretical principles underlying this mechanism are independent of model details (instead they depend on the dynamical excitability class). Weak electrical coupling between neurons with SNL excitability therefore provides a generic mechanism that suffices to produce the observed desynchronized (splay) states in small networks. Importantly, this mechanism holds across different firing rates and does not require a pronounced AHP, but the AHP may further stabilize desynchronization. Although previous studies have provided mechanisms by which networks with electrically coupled neurons can be configured into desynchronized, phase-locked firing 41 -such as inhibitory synapses overriding the synchronization by electrical synapses 42 or interacting with the spiking synchrony 43,44 , silencing of electrical synapses in a network with mixed chemical and electrical synapses 42,45 , or strongly asymmetric electrical synapses and specific input regimes 46desynchronization by the electrical synapse alone, without the need of additional chemical synapses, other network motives or feedback from the network, is conceptually new.

Sequence preference requires heterogenous coupling
Our in vivo recordings show that MN1-5 firing is not only desynchronized, but organized into preferred sequences of cellular activation (within the splay state; Fig. 1f,g). This raises of the question of how the preferred sequences are generated. Our modelling (5 min simulations in network models with SNL neurons and noise) demonstrates that homogenous electrical coupling does not show sequence preference ( Fig. 3j (left)). By contrast, heterogenous electrical coupling between MN1-5 ( Fig. 3i) as observed in situ (Fig. 2f) produces the same preferred splay states with similar sequence statistics (Fig. 3j (right)) as observed in vivo across animals (Fig. 1f). We conclude that heterogenous weak electrical coupling among neurons with HOM excitability (near SNL) sufficiently explains the preferred splay states in this small network.

Splay state serves stable wing power
MN1-5 splayed-out firing in multiple species suggests that it provides a functional benefit. We tested this hypothesis by comparing the wingbeat frequency between normal splay state and synchronous MN firing, the latter induced by increasing electrical coupling through genetic manipulation (Figs. 2a and 3c). In vivo, synchronous MN firing causes 8× higher wingbeat frequency fluctuations over time compared with firing in splay state (Fig. 4a). The amplitudes of such fluctuations are close to changes in wingbeat frequency in response to optomotor input 1,16 and are therefore functionally relevant for flight altitude control. In the steady-state, muscular stretch sensitivity and wingbeat frequency are directly proportional to the myoplasmic calcium levels 1 , but the Article dynamics of myoplasmic calcium in DLM fibres is unclear. To link MN firing patterns to wingbeat frequency, we measured the kinetics of muscle fibre membrane potential and myoplasmic calcium concentration changes after MN1-5 firing (Fig. 4b). An MN spike causes a rapid depolarization of the DLM fibre membrane (Fig. 4b (top)), which is followed by a myoplasmic calcium signal with a rise time constant of around 6 ms and a decay time constant of about 80 ms (Fig. 4b (middle)). If muscular stretch sensitivity was controlled by the myoplasmic calcium levels with a ms temporal resolution, changes in wingbeat frequency during simultaneous MN firing should follow the same time course as changes in myoplasmic calcium across all six DLM fibres. This is precisely what we found (Fig. 4b (bottom)). Synchronous MN firing causes peak wingbeat frequencies after 7.3 ± 1.3 ms (myoplasmic calcium rise time constant τ = 6.2 ± 1.7 ms) that decline with a time constant τ = 83 ± 27 ms (calcium decay time constant is 82 ± 6.7 ms).
To determine the kinetics of myoplasmic calcium also at different MN firing frequencies, we stimulated a single MN electrically at different frequencies (0.5, 1, 2, 5, 10 Hz), spanning the normal range during flight, and imaged the resulting calcium signals in the respective DLM target fibre. Discrete myoplasmic calcium signals were observed for the entire frequency range tested ( Fig. 4c and Supplementary Video 4), although summation as caused by an incomplete return to resting calcium levels between MN spikes starts at 2 Hz and increases at higher frequencies.
However, for firing frequency ranges of MN1-5 that are observed during normal flight (Fig. 1c), in each muscle fibre, myoplasmic calcium fluctuates around each corresponding MN spike. Given that we did not observe any differences in the temporal coupling of MN firing and myoplasmic calcium across different stimulation frequencies or different DLM fibres, these data can be extrapolated to all DLM fibres for splay state (Fig. 4d (top)) as compared to synchronous MN firing (Fig. 4d (bottom)). Given that each DLM fibre is innervated by one MN only (Fig. 1a), MN firing in splay state causes splayed-out calcium signals across the different muscle fibres (Fig. 4d (top)), each peaking around 6 ms after the respective MN spike and declining with a time constant of around 80 ms (Fig. 4b (middle)). The bottom trace shows the resulting average myoplasmic calcium across fibres (Fig. 4d (top)). By contrast, during synchronous firing of MN1-5, the calcium signals in all DLM fibres are time locked, which results in much larger calcium fluctuations across DLM fibres (Fig. 4d (bottom)). This extrapolation of the average myoplasmic calcium levels across the DLM fibres predicts eightfold smaller fluctuations in wingbeat frequency over time in normal splay state than for synchronous MN firing. We confirmed this prediction by averaging changes in wingbeat frequency around MN spikes for splayed-out ( Fig. 4e (blue line)) as opposed to synchronous firing ( Fig. 4e (green line)). Thus, the splay state serves to minimize wing power fluctuations.
Taken together, splay-state asynchronous flight motor patterns are conserved across individuals and species, are produced by a minimal CPG of weakly electrically coupled MNs, and serve constant wingbeat power output at a given power demand. This provides a comprehensive view of the asynchronous flight CPG network structure and the resulting functional consequences for one of the most abundant forms of locomotion on earth. Moreover, we provide a theoretical background for firing d esynchronization/synchronization in small gap-junctional networks in the context of their specific cellular excitability. The underlying mechanism is generic, predicting desynchronizing functions of electrical synapses beyond the motor system of insects. Electrical synapses can therefore be used for operations such as sign-reversal and ensemble firing desynchronization, a functional versatility and impact on neural circuit dynamics that was previously attributed to chemical synapses alone.

Online content
Any methods, additional references, Nature Portfolio reporting summaries, source data, extended data, supplementary information, acknowledgements, peer review information; details of author contributions and competing interests; and statements of data and code availability are available at https://doi.org/10.1038/s41586-023-06099-0.   Publisher's note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

Animals
D. melanogaster were reared in 68 ml transparent food vials (Kisker Biotech, 25 mm × 95 mm) filled with 10 ml of standard cornmeal/ glucose/yeast/agar diet at 25 °C and 60% humidity under a 12 h-12 h light-dark cycle. For experiments, adult 2-5-day-old male flies were used. We used Canton-Special (Canton-S) as the wild-type control. To selectively target the DLM MNs (MN1-5) with UAS-driven transgenes we used a split GAL4 driver line that combines GMR23H06 (BDSC, 49050 discontinued) and GMR30A07 (BDSC, 49512) from the Rubin Collection 47 . The GMR23H06 enhancer expresses the activating domain (AD) and the enhancer GMR30A07 expresses the DNA-binding domain (DBD). Overlap of both expression patterns-and, therefore, functional GAL4 expression-is restricted to the MN1-5 (Extended Data Fig. 5a,b). The following transgenes were expressed in MN1-5 to manipulate gap junction expression and function. shakB RNAi (BDSC, 57706) targets all shakB isoforms and was used to KD gap junctions in MN1-5. For gain of function through overexpression of gap junction protein the shakB(N+16) isoform was used (gift from P. Phelan). KD of inhibitory synaptic transmission through GluCls or RDL GABA-ARs was conducted by targeting the expression of UAS-RNAi KK for GluClα (VDRC, 105754) or UAS-Rdl-RNAi (VDRC, 41103) to MN1-5 with the DLM MN GAL4 driver GMR23H06. To test for electrical synapse function by genetic manipulation, UAS-RNAi-KD of shakB was compared with UAS-ShakB overexpression and empty UAS-control with 7 biological replicates in each group. To test for effects of GABA-ARs and GluClRs on flight patterns, UAS-RNAi-kd for each receptor type was expressed in MN1-5 and the effect on flight patterns was compared with the respective genetic controls, with 8 biological replicates for GluClRs and 10 replicate animals for GABA-ARs. For all genetic manipulations, the order of experiments was fully randomized and analysis was conducted without knowledge of the genotype. Other insect species were obtained from pet shops as feed insects and Apis mellifera (at 4 to 6 weeks of age) from the local apiculture of the Johannes Gutenberg University Mainz.
For in vivo optogenetic manipulation of the activity of MN1-5 during flight, we expressed UAS-XXL channelrhodopsin (XXL-ChR) (BDSC, 58374) 48 either under the control of our DLM-MN-specific split GAL4 driver (split GMR23H06+GMR30A07) to enable light activation of MN1-5 (7 replicate animals) or, alternatively, in presynaptic cholinergic neurons under the control of Cha-GAL4 (11 replicate animals). Direct light activation of MN1-5 during flight was used to test whether the CPG correctly shapes MN1-5 firing into preferred sequences even when firing rates are artificially increased. Light activation of presynaptic cholinergic neurons during flight was used to test whether the CPG translates artificially presented, unpatterned synaptic input into patterned output from MN1-5. For in vivo calcium imaging in DLM fibres (6 replicate animals) the genetically encoded calcium indicator UAS-GCaMP8f 49 was expressed under the control of Act88-GAL4 50 .
We used targeted expression of FMRP RNAi in MN1-5 (w;UAS-FMRP RNAi / 23H06-ADZ UAS-CD4-td-GFP;30A07-DBD/+) to visualize weak electrical coupling by dye coupling. A reduction in FMRP is reported to increase both the uptake of small tracers into Drosophila neurons during iontophoresis and, consequently, also increase dye uptake into neurons that are electrically coupled through shakB-encoded gap junctions 28 .
No ethical oversight was required for work with insect species. The work was according to the guidelines of research with transgenic invertebrate animals of the state of Rhineland-Palatinate in Germany.

Generation of UAS-Shab flies
Shab cDNA was provided by F. Sigworth 51 in pBluescript KSM. Shab was then PCR amplified with the following primers: TTCAGGCGGCCGCGGCTCGAGaacttaaaaaaaaaaatACAAATGGTCGG GCAATTGCAAGG (forward primer; the XhoI restriction site is underlined, the Shab start is bold and underlined) and CCTTCACAAAGATCCTCTAGA gactcactatagggcgaattgg (reverse primer; the XbaI restriction site is underlined, the Shab terminal sequence is bold and underlined). The primers included overhangs matching the pJFRC81 attB-site-directed expression vector (36432, Addgene). pJFRC81 was cut with XhoI and XbaI restriction enzymes according to the provider's protocol (New England Biolabs, R0146S, R145S) to remove the GFP sequence from the vector. Using NEBuilder according to provider's instructions (New England Biolabs, E5520S), the amplified Shab cDNA was inserted into the cut pJFRC81 vector, thereby replacing the GFP. The PCR-amplified Shab cDNA was not purified before assembly with NEBuilder. The assembled construct was then transformed into Escherichia coli Dh5α bacteria and streaked out onto ampicillin (50 μg ml −1 ) containing LB media agar plates. Overnight colonies were picked and grown in ampicillin containing LB media (50 μg ml −1 ). DNA was isolated using the Maxiprep kit (Qiagen, 12162). Using specific primers, we determined that the cDNA used was Shab splice variant K (Supplementary Table 5). The correct sequence was confirmed by sequencing. For generation of transgenic flies, the Shab-pJFRC81 attB vector was injected into fly embryos at the attP2 landing site by BestGene. For experiments, homozygous male w;;10xUAS-Shab-pJFRC81attP2 flies were crossed to female w;GMR23H06-ADZ;GMR30A07-DBD (DLM-split GAL4 driver) flies and male F 1 progeny were tested.

Electromyography/extracellular recording of MNs
The in vivo activity of MN1-5 during tethered flight was monitored by extracellular voltage recordings from their respective muscle fibre. One MN action potential causes one large, transient postsynaptic depolarization ( Fig. 4b (top trace)) in its target DLM fibre (for innervation, see Fig. 1a). As each DLM fibre is innervated by one MN only, and each MN spike causes the same distinct evoked postsynaptic response, the EMG spikes recorded from the DLM fibres are a 1:1 reflection of the activity patterns of MN1-5.
To prepare flies for in vivo EMG recordings during tethered flight, animals were briefly cold-anaesthetized (20 s in an empty 68 ml plastic vial on ice), transferred dorsal side up onto a cold metal plate (~3 °C) and glued (clear glass adhesive (Duro; Pacer Technology)) with head and thorax to a triangle-shaped tungsten hook (0.1 mm diameter). Curing of the glass adhesive was induced by exposure to ultraviolet light (Mega Physik Dental Cromalux-E Halogen Curing Light Unit) for 45 s. The flies were kept individually for 10 min to recover from the cold anaesthesia. The rested flies were then mounted in the set-up to a clamp attached to a micromanipulator. Tarsal contact to a small polystyrene bead prevented unwanted engagement into flight. Next, a light barrier, consisting of a red laser, an aperture to reduce beam width to approximately 1 mm and a LED light sensor, was positioned so that the light beam was broken by each wingbeat. After positioning of this laser-based wingbeat detector, tungsten electrodes for extracellular MN recordings were inserted. A reference electrode was inserted into the last abdominal segment. To record the activity of MN5 and MN4 from the dorsal-most DLM fibres, we inserted one sharpened tungsten wire into the dorsal thorax in front of the anterior dorsocentral bristle, making sure not to cross the midline of the thorax. This resulted in extracellular recordings of two units that could easily be distinguished by amplitude and shape. Confirmation of unit identity was achieved by adjusting the tungsten insertion depths and monitoring the resulting amplitude changes of the recorded units. As MN5 innervates the dorsal-most DLM fibres, it is the first unit to appear after shallow dorsal electrode insertion. Moving deeper lets a second unit appear that belongs to MN4, which innervates the fibre just beneath. This MN4 unit increases in amplitude with increased insertion depth in DLM4, while the amplitude of MN5 decreases. This procedure allowed unambiguous allocation of the two units to MN5 and MN4. To record up to 5 DLM units (MN1-5) simultaneously, additional tungsten electrodes were inserted. A second electrode was inserted four small bristles in front of the first electrode, which recorded MN4 and MN3. A third electrode was inserted anterior on the same line as electrodes 1 and 2 at the position where small bristles begin to appear on the thorax. Again, unit identification was conducted by altering electrode insertion depth at the beginning of each experiment. Starting from dorsal and slowly moving ventrally, the first unit to appear is MN5, the second one is MN4, the third one is MN3 and so on. Electrode depths were adjusted until the subsequent unit that appeared while moving deeper had approximately half of the amplitude of the first unit, so that both were easily distinguishable. In some recordings, it was also possible to pick up three clearly separable units simultaneously, but spike sorting was more time consuming.
All extracellular recordings were amplified at 1,000× using the AM-Systems 1700 extracellular amplifier. Electrophysiological and wingbeat recordings as well as a frame trigger signal from high-speed video (see below) were digitized with the Axon Digidata 1550B (Molecular Devices), acquired with AxoScope (v.10.7) and stored on a PC. Spike sorting was conducted offline with the spike sorting function in Spike2 (v.7.2). Data were further analysed using Spike2, and custom Python routines were created using Jupyter notebook to count pattern probabilities and create phase histograms. For additional functions, the Python libraries NumPy, pickle, SciPy, Matplotlib and seaborn were imported. All of the Python scripts are available on request. Additional data analysis and statistics were conducted using Microsoft Excel Professional Plus 2019 (v.2110) and GraphPad Prism (v.9.2.0).
All of the animals that were successfully recorded for at least 10 min of flight were included in the analysis of phase relationships. For spike processing, the DLM recordings (.abf format) were imported into Spike2 and sorted with the wavemark function, therefore creating templates for each unit. Template-matching spikes were visualized separately and synchronous spikes (the sum of the amplitude and waveform of two distinct MN units) were split into the two according units. The timepoints of the peaks of all marked spikes were then imported into single-event channels. The first and last 3 s of each flight bout were discarded to restrict the analyses to stable in-flight motor patterns. At the beginning of flight bouts, the MN spikes occur at higher frequency and with the highest irregularity, whereas, at the end of a flight bout, the firing frequency decreases. Phase histograms were created from the event channels. Each ISI of an event channel with the timepoints of the units of a given MN was compared to when the units of a second MN occur within the ISIs and counted into bins. All ISIs were normalized to 1 (phase 0 to phase 1), divided into 100 bins, the events of the second MN unit were counted into the respective bins and divided by the total number of events (relative occurrence; Extended Data Figs. 1, 2, 6 and 8). To create phase histograms with the phase 0 in the middle of the diagram (Fig. 2b and Extended Data Figs. 3b, 5e and 6e,f), the last 50 bins were copied in front of the first 50 bins and the x axis was changed to range from phase −0.5 to phase 0.5.
All insect species larger than Drosophila were cold anaesthetized for up to 20 min in the fridge at 5 °C and fixated in a 3D printed device to hold the insect in place. The .stl file is available on request. Instead of a rectangular tungsten hook, the wire is bent in a half circle and glued onto the outer circumference of the posterior dorsal region of the thorax to increase the adhesive surface.

High-speed video
In addition to electrophysiological recordings of MN activity and laser-based wing detection, selected times of tethered flight were simultaneously recorded as high-speed video using the Photron FAST-CAM Mini UX100 camera and Photron FASTCAM Viewer software (PFV, v.3.6.9.0) at 5,000 fps. As illumination, two infrared lights (Sygonix IR illuminator with 48 LEDs) were used. Videos of Drosophila (Supplementary Video 1), goldfly (Supplementary Video 2) and honey bee (Supplementary Video 3) during tethered flight are provided.

Electrolytic sharpening of tungsten electrodes
The tips of tungsten rods (diameter of 100 μm for Drosophila, 125 μm for all larger insects) were sharpened electrolytically by repeated dipping in a NaNO 2 KOH solution (10.3 M NaNO 2 and 6.05 M KOH in double-distilled H 2 O). Tungsten rods of 1 cm length were crimped onto metal rods. For this, the wire was placed onto the rod and covered with a 0.5 mm ferrule. A crimping tool compressed it to a pluggable tungsten wire electrode. Now the electrode was placed into an alligator clip, connected to a stimulator (Grass SD9 square pulse stimulator) and repeatedly dipped into the sharpening solution, which was connected to the other pole of the stimulator. A monophasic current with 100 Hz, 40 V and 1 ms duration was applied. The tip was frequently moved in and out of the solution to form a thin tip. Finally, electrodes were rinsed with distilled water.

Double patch-clamp recordings of MN pairs
After removing the legs and wings, 2-3-day-old adult Drosophila males were pinned into a Sylgard-coated lid of a 35 mm Petri dish and fixed ventral side up with minute pins in the head and through the tip of the abdomen. After submerging in normal saline (128 mM NaCl, 2 mM KCl, 1.8 mM CaCl 2 , 4 mM MgCl 2 , 5 mM HEPES and ~35.5 mM sucrose, the pH was adjusted to 7.24 with 1 N NaOH, and osmolality was 300 mOsm kg −1 ) the thoracical cuticle was removed with fine iris scissors to expose the VNC. The specimen was then rinsed thoroughly with saline. After mounting the preparation onto the stage of an upright Zeiss Axio Examiner epifluorescence microscope, the VNC was viewed with a ×40 water-immersion lens (Zeiss W Plan Apochromat ×40/1.0 NA, DIC VIS-R), the UAS-6xmCherry-expressing MN1-4 (genotype: w;23H06-ADZ UAS-6xmCherry/+;30A07-DBD/+) were viewed through a TRITC filter set. The ganglionic sheath and debris hampering access to MN1-4 somata were removed by repeated application of 1% protease type XIV (from Streptomyces griseus, Sigma-Aldrich, P5147) through a patch pipette with a manually broken tip that was then also used to remove loosened debris 52 . Before recording, the specimen was washed with saline for 5 min through a gravitation perfusion system at around 2 ml min −1 , the bath volume was around 300 μl. Recordings were performed using patch pipettes (borosilicate glass capillaries; outer diameter, 1.5 mm; inner diameter, 1 mm, without filament, World Precision Instruments, PG52151-4) pulled with a PC-10 vertical puller (Narishige) filled with internal patch solution (140 mM K-gluconate, 2 mM Mg-ATP, 2 mM MgCl 2 , 11 mM EGTA, 10 mM HEPES. pH was adjusted to 7.24 with 1 N KOH, osmolality was adjusted to 300 mOsm kg −1 with glucose if necessary) that approached the two MNs from opposite sides. Patch pipette tip resistance was between 5 and 6 MΩ with these solutions. Double patch recordings were performed by connecting each patch pipette to a separate Axopatch 200B amplifier (Molecular Devices). Data were filtered at 5 kHz through a lowpass Bessel filter, digitized through an analogue/digital converter (Digidata, 1440), and signals were recorded using pClamp10.7 software (both Molecular Devices). Output gain was 10×. For double recordings, one MN was approached with a patch pipette. After giga seal formation, pipette capacitance artifacts were cancelled manually and whole-cell configuration was established at a holding potential of −70 mV. After setting whole-cell capacitance compensation, correction and prediction values as well as series resistance compensation using the respective dials of the amplifier (only used to judge on recording quality as all of these compensations are disabled in current clamp mode), the recording was disconnected from external influence by using the I=0 setting of the amplifier, making it possible to revert to bath mode without disturbing the already established recording of the first MN. Then, recording of the second MN was established identically. To monitor both recordings, separate channels of the digitizer and in the pClamp10.7 software were used. Quality parameters were as follows: giga seal, >5 GΩ; membrane potential of −70 mV was held with a holding current smaller than ±100 pA; series resistances of >15 MΩ were not accepted to ensure good control when applying current injections; only if these criteria were met for both MNs, the recordings were switched to current clamp mode. Resting membrane potential was around −60 mV without current injection. To determine coupling strength and rectification parameters of gap junctions between MNs, depolarizing and hyperpolarizing current was injected into one MN while the other was monitored simultaneously, and vice versa. Slow tonic firing was induced in one or both MNs at rates of between 3 and 8 Hz, as is observed during flight behaviour, by small somatic current injections while monitoring the respective other MN. For input-output relationships, firing was induced by 1,000 ms square pulse current injections up to 1 nA in 0.1 nA increments.

Intracellular dye filling
Intracellular dye fills with neurobiotin were used to detect dye coupling between electrically coupled MNs. For intracellular dye fills, FMRP RNAi was targeted to MN1-5 (w;UAS-FMRP RNAi /23H06-ADZ UAS-CD4-td-GFP;3 0A07-DBD/+), because reductions in FMRP have been shown to increase both neurobiotin uptake into Drosophila neurons during iontophoresis and, consequently, also increased dye uptake into neurons that are electrically coupled through shakB-encoded gap junctions 28 . After removing the legs and wings, 2-3-day-old adult male Drosophila were pinned in a Sylgard-coated lid of a 35 mm Petri dish and fixed dorsal side up with two minute pins, one through the head and one through the abdomen. After submerging the specimen in normal saline (128 mM NaCl, 2 mM KCl, 1.8 mM CaCl 2 , 4 mM MgCl 2 , 5 mM HEPES and ~35.5 mM sucrose, pH was adjusted to 7.24 with 1 N NaOH, and osmolality was 300 mOsm kg −1 ), it was opened along the dorsal midline up to the neck connectives with iris scissors. The cut DLM was then pinned to the sides with one minute pin each to expose the gut, inner organs and VNC. The gut, salivary glands as well as other inner organs were removed to fully expose the VNC. The specimen was then rinsed thoroughly with saline to remove excess debris. After mounting the preparation onto the stage of an upright Zeiss Axio Examiner epifluorescence microscope, the VNC was viewed using a ×40 water-immersion lens (Zeiss W Plan Apochromat 40x NA 1.0, DIC VIS-R), the UAS-CD4-td-GFP-expressing MN1-5 were viewed through a FITC filter set. The ganglionic sheath was removed focally using a broken patch pipette filled with 1% protease type XIV (from Streptomyces griseus, Sigma-Aldrich, P5147) to allow access to MN somata 52 . Dye fills were performed using sharp glass microelectrodes pulled from filamented borosilicate glass capillaries (Sutter BF100-50-10; tip resistance, ~40 MΩ) with a Sutter P-97 Flaming Brown microelectrode puller. The tip was filled with a 50/50 mixture of TRITC-Dextran 3000 lysin fixable (Invitrogen, 3308) and neurobiotin (Vector Labs, SP-1120) dissolved in 2 M potassium acetate, and the shaft was then filled with 2 M potassium acetate leaving an air bubble between the dye-filled tip and the shaft to avoid dilution of the dye. MNs were impaled by a short buzz (which makes the electrode tip vibrate for a set amount of time, here around 40 ms) with a remote buzz connected to an Axoclamp 2B intracellular amplifier (Molecular Device) that was also used for dye filling. MNs were filled iontophoretically by positive current injection (between 0.5 and 1 nA) in bridge mode until the MN was judged to be filled, after around 10 min. The microelectrode was then removed, and the specimen was fixed with 4% paraformaldehyde at room temperature for 50 min, without shaking. This was followed by three rinses with PBS and three washes for 20 min PBS and six washes for 20 min with 0.5% PBS-Triton X-100 (PBT, Triton X-100, Sigma-Aldrich, T8787), with shaking. The preparation was then incubated with Streptavidin coupled to Alexa 647 (Thermo Fisher Scientific, S-21374) in 0.3% PBT at room temperature in the dark for 2 h with shaking. Streptavidin Alexa 647 was then removed, the preparation was rinsed several times with PBS and then washed three times for 20 min with PBS at room temperature in the dark with shaking. This was followed by an ascending ethanol series: 50%, 70%, 90%, 100% for 10 min each. The preparation was then mounted in methylsalicylate, topped with a high-precision (170 ± 5 μm) cover slip and sealed with clear nail polish. The dye fill was then visualized using the Leica TCS SP8 confocal laser microscope with a helium-neon laser. Alexa 647 was excited at 633 nm and emission was detected between 650 and 680 nm with a photomultiplier tube. Images were taken with a ×40 oil objective (NA 1.3) with a 1.75 digital zoom at a resolution of 1,024 × 1,024 pixels and a z-step size of 1 μm.

Optogenetics
For optogenetic stimulation, UAS-channelrhodopsin was either expressed under the control of a DLM-MN-specific Split-GAL4 line (GMR23H06-ADZ attP49; GMR30A07-DBD attP2; see Extended Data Fig. 5a,b for expression patterns specifically in MN1-5), or under the control of ChaT-GAL4 in cholinergic neurons. We tested multiple different optogenetic transgenes. UAS-CsChrimson for red light activation caused few MN spikes followed by silence, probably caused by a depolarization block, and was therefore not suitable for optogenetic stimulation during ongoing flight. By contrast, blue light stimulation of channelrhodopsin XXL reliably increased MN firing frequencies during flight, and therefore was the most effective transgene for optogenetic activation experiments during flight. As a control, the blue light stimulation was applied in an identical manner to flies without expression of UAS-channelrhodopsin XXL (DLM-MN split-GAL4 and ChaT-GAL4 crossed to w1118 flies, respectively). Blue light stimulation alone had no effect on MN firing frequencies during flight. The light stimulus was applied though fibre optics directed to the ventral thorax from underneath to avoid stimulation of cholinergic neurons in the brain. The light stimulus was presented constantly so that optogenetic activation of either cholinergic pre-motor neurons (Extended Data Fig. 3) or MN1-5 (Extended Data Fig. 5c,d) was not patterned. Moreover, channelrhodopsin XXL activates within milliseconds after blue light stimulation but has very slow inactivation kinetics with an off time constant of 76 ± 12 s (ref. 48). Thus, increased firing rates during flight outlast the blue light stimulus (Extended Data Figs. 3 and 5c).

In vivo calcium imaging from DLM fibres during flight
For in vivo calcium imaging in DLM fibres, the genetically encoded calcium indicator UAS-GCaMP8f 49 was expressed in muscle under the control of Act88-GAL4 50 . All earlier GCaMP versions exhibit signal decay kinetics that are too slow to measure the time course of decay of myoplasmic calcium signals as induced by a MN action potential. Calcium imaging was conducted using a Hamamatsu CMOS (C11440-42U) camera mounted to a Zeiss Axioscope 2FS fluorescence microscope equipped with a ×10 lens that allowed enough working distance for wingbeats under the microscope. Images were acquired at a frame rate of 98 Hz with 16 bit image depth at a 512 × 512 pixel spatial resolution. Myoplasmic calcium responses in DLM fibres were imaged in vivo through the dorsal cuticle either during tethered flight, or after electrical stimulation of selected MNs (Supplementary Video 4 shows myoplasmic calcium responses of DLM fibre 6 in response to electrical stimulation of the MN5 that innervates DLM fibres 5 and 6; Fig. 1a). Myoplasmic calcium signals were acquired using Hamamatsu photonics software (TOKUPIC v.1.0) and presented as percentage change (∆F/F) as previously described 21 .

Computational modelling and equations
Conductance-based MN model and gap junctional circuit model. Individual MNs are described by a single-compartment model with spike generating Na + and K + currents, I Na and I K , based on ref. 39 and further refined on the basis of in situ current clamp and voltage clamp recordings from our laboratory (see below). For I Na, the halfactivation voltage of the activation gate is −33 mV, the half-activation voltage of the inactivation gate is −39.14 mV and the reversal potential E Na is 55 mV. For I K (referred to as I shab in the equations below) the half-activation voltage is −42.14 mV and the reversal potential E K is −72 mV. The maximal sodium conductance g Na is 0.4312 μS and the maximal potassium conductance g shab is varied between 0.13768216 μS and 0.34496 μS throughout the paper (see the 'Shab-induced bifurcations' section; a list of single-neuron-model parameters is provided in Supplementary Table 2). For the gap junctional circuit model, identical single-neuron models were coupled by linear non-rectifying gap junction currents, I gap, with properties based on dual patch current clamp recordings from MN pairs. The current-balance equation of the i th neuron reads where v is the membrane voltage and C m is the membrane capacitance. The leak current I L with the leak conductance g L , the potassium current I shab , and the sodium current I Na are defined as Gates of the voltage-dependent ion channels follow first-order kinetics The bidirectional, non-rectifying gap junctions with linear charge transfer (Fig. 2e,f) are described by the following current Parameters common to all simulations are provided in Supplementary  Table 2. For simulations with homogeneous coupling (Fig. 3b-d,j (left)) a coupling strength of g = 43.5 pS ij gap for all i ≠ j was used. To keep the cases as comparable as possible, the mean coupling strength in the simulations with heterogeneous coupling (Fig. 3h,j) was chosen to be the same as in the homogeneous case. The ratios of coupling coefficients for MN1-4 were chosen to be the same as measured in vivo (Fig. 2f).  53 . For stochastic simulations, the Heun method was used with a time step of 3 μs. For noisy simulations, a zero-mean white-noise current was added to the current balance equation (using the brian2 xi variable) with a noise strength of σ = 0.949 pA ms (Fig. 3d,h,j). For the deterministic simulations in Fig. 3b,c, the classical Runge-Kutta method (RK4) was used with a time step of 100 μs.
In addition, Q represents e k T B , which is the elementary charge e divided by the Boltzmann constant k B and the temperature T in Kelvin. Throughout the simulations, Q was approximated as 39.2/V.
The model parameters common to all models are summarized in Supplementary Table 2; the activation curves of the gates and the kinetic time scales are provided in Supplementary Table 3. Some parameters and equations were reported inconsistently in the original publication. Because the authors 39 provided their code, we were able to extract the parameters used for their simulations, which are reported here. These updated parameters were validated using additional current-clamp and voltage-clamp recordings from our laboratory. Furthermore, to produce different onset bifurcations, the Shab conductance was changed (see below) and the input current was adjusted to move the system close to onset. The values differing between models are summarized in Supplementary Table 4.
For almost all simulations, neurons were initiated in random phases, by drawing timepoints from within one ISI and initializing the simulation with the corresponding state vector of a periodically spiking uncoupled neuron. For the simulations in Fig. 3c, the same initial conditions were used for each coupling strength. For the example in Fig. 3b, the initial conditions were not random, but hand-picked to yield a sequence that matches one of the preferred sequences in Fig. 1f.
Gap junctions that adhere to equation (1) affect the intrinsic properties of neurons because they act as effective leak and capacitance changes. For example, gap-junction currents can influence firing frequencies depending on the synchronization state of the network. In Fig. 3c, the coupling is homogeneous g g = i j gap , gap for all i ≠ j between pairs but the coupling strength g gap is varied. The effect of coupling on firing frequencies is strongest in the intermediate coupling range. In the weak coupling range, the effect on frequency is small because the current flow through gap junctions is small. In the strong coupling range, the neurons synchronize and, consequently, almost no gapjunctional current flow occurs. As frequency homeostasis is not part of this model, we accept biologically implausible firing frequencies in the intermediate range. In biological neurons, a wide range of mechanisms-for example, based on additional ion channels, energy consumption or ionic concentrations-could stabilize the frequencies 55 .

Shab-induced bifurcations.
A local bifurcation analysis of the fixpoints in the model was shown previously 39 . The important quantity for predicting network states of the CPG is the PRC (Fig. 3e), which is closely related to the spike-onset bifurcation of the neuron and, therefore, the neuron's excitability class. The bifurcation diagram in Extended Data Fig. 9a adds non-local and codim-2 bifurcations that organize transitions between excitability classes and their PRCs such as big SNL and small SNL 37,56 . Depending on the Shab channel density, g shab on the ordinate, different onset bifurcations that change the resting state into the spiking mode occur. In the middle range of g shab levels, spiking commences through a well-known SNIC bifurcation. The SNIC region is encapsulated by two SNL points. At the lower end, for g shab smaller than the small SNL (sSNL) point, spiking commences through a small homoclinic loop (HOM, green line in Extended Data Fig. 9a) at lower input currents, I in , than the saddle-node bifurcation (blue line in Extended Data Fig. 9a). At the upper end of the SNIC range, a series of bifurcations is traversed that eventually leads to spike onset through a fold of limit cycles (FLC, bifurcation line not shown) and a subsequent subcritical Hopf bifurcation, here termed the Hopf regime for simplicity. Bifurcations leading from SNIC to the Hopf excitability class are the big SNL, the neutral saddle loop (not shown) and also a cusp 56 . For simplicity, in the main text, when using the term SNL, we refer to sSNL. The Bogdanov-Takens point in between creates an Adronov-Hopf line (red line in Extended Data Fig. 9a) that eventually turns back to create the excitation block at higher input levels.
The consequences for the phase relationships of coupled neurons can be understood from the neuron's PRC at different bifurcations (Extended Data Fig. 9a). A neuron's ability to produce arbitrarily low firing rates depends on its excitability class (determined by the spike onset bifurcation). So does the PRC shape [57][58][59] . The insets in Extended Data Fig. 9a show prototypical PRCs for the three onset bifurcations, which can be achieved at different levels of Shab-channel density. The exact parameters are summarized in Supplementary Table 4. A critical transition with repercussions for the CPG network states is the symmetry breaking in the PRC at the sSNL bifurcation as described in the next section 59 .
An FLC/Hopf onset bifurcation as spike onset can be ruled out as a model for the biological MN1-5 network based on the following grounds: the region in which tonic spiking commences through the FLC/Hopf bifurcation has a finite, non-zero frequency at the rheobase. For the present model (Extended Data Fig. 9a (top arrow)), this is at around 47 Hz and therefore above the dynamic range of the in vivo measures f-I curves from the MNs, which operates between 0-40 Hz (3-15 Hz during normal flight; Fig. 1e).
The bifurcation analysis was performed using AUTO-07P 54 .

Mechanism of splay state generation.
Our theoretical analysis predicts that dynamics close to the SNL point favour splay states 59 . This prediction is based on coupling functions, which are obtained from the biophysical conductance-based models using a phase-reduction (see the next section). The coupling function G(ψ) combines a neuron's PRC and a synaptic transfer term, which, for the system at hand, captures the voltage perturbations caused by gap-junctional coupling (equation (4)). The coupling function describes how one neuron's phase is shifted based on the phase difference to another neuron. Of particular interest is the odd part of the coupling function, G odd (ψ) = G(ψ) − G(−ψ). It reflects the reciprocity of coupling between two neurons and extracts how coupling affects the phase difference in the cells' firing. Phase relationships of tonically firing neurons can be read from G odd (Fig. 3e (bottom)).
Crucially, stable phase relationships depend on a stable fix point in this odd part of the coupling function, that is, a zero crossing with a negative slope. The shape of the coupling function in turn depends on the bifurcation type of the neuron (Extended Data Fig. 9a). For example, for a PRC of a neuron with SNIC dynamics, a small phase distance, Δφ, (neuron 1 fires shortly before neuron 2) is further decreased, whereas a large Δφ (neuron 2 fires shortly before neuron 1) is further increased. This results in a stable fixpoint at phase 0, so that the neurons show synchronized firing (Fig. 3e (bottom)). By contrast, with a PRC at the SNL point, a small Δφ is increased, whereas a large Δφ is decreased. This results in a stable fixpoint at Δφ = 0.5 for each coupled neuron pair. For a network of >2 electrically coupled neurons, this causes a frustrated state because antiphase locking at Δφ = 0.5 cannot be achieved for all units at the same time. In a small network, such as the MN1-5 network with its five coupled neurons, the splay state represents a low-frustration solution that determines the most likely network state 60 .
These analytical predictions based on coupling functions were corroborated by simulations of the full biophysical model. Near the SNL point, networks of five gap-junction-coupled neurons indeed showed a splay state (Fig. 3b), whereas neurons with a SNIC onset synchronized instead (Fig. 3h). These simulations show that, for our parameters, the phase reduction is valid and that coupling functions based on the phase reduction accurately capture the network's synchronization behaviour. On the basis of both coupling functions and simulations, we predict that the combination of coupling and intrinsic excitability determines the splay state in the biological MN1-5 network.
Phase-reduced circuit model with repulsive coupling. To generate coupling functions, we reduce the biophysical, conductance-based neurons to approximate phase oscillators 61 . This is warranted because the cholinergic inputs drive MN1-5 into a tonically spiking regime. The i th phase oscillator is defined by its intrinsic, mean firing rate, f i , and its phase response curve, Z(φ i ), obtained from the biophysical model (Extended Data Fig. 9a). The PRCs are calculated based on a direct perturbation approach 62,63 . The network equation reads: where g(φ i ,φ j ) is the phase-dependent perturbation received through the gap junctions. It depends on the phases of the presynaptic and postsynaptic neurons. If the spike waveforms are uniform, the interaction term, based on equation (1) reads To understand the influence of the PRC on phase relationships in small networks, first consider the phase difference between two coupled neurons i and j, isolated from the network. Define their phase difference ψ = φ i − φ j . Then, using approximate averaging theory, the slow evolution of the phase difference simplifies 33 to odḋ Here, the small frequency detuning is ν = f i − f j and the averaged coupling function is defined as Stable fixpoints of equation (3) determine constant phase relations for neurons with similar intrinsic frequencies as illustrated in Extended Data Fig. 9b,c. In Fig. 3e the odd part, G odd (ψ), and the fixpoints are shown for different choices of Z(φ). The stable in-phase fixpoint directly translates to a synchronous network state also for network sizes of N > 2. However, this is not observed in the CPG recordings (Fig. 1a). For the network to show frustration, pairs of neurons must be phase-repellent.
The top spike trains in Extended Data Fig. 9b illustrate an example that generates stable antiphase fixpoints. To explain stable phase relationships graphically, the model was simplified, such that the shape of the coupling is a scaled version of the PRC (corresponding to gap junction coupling with delta spikes). Moreover, the mean PRC was absorbed into the average firing rate. The stability can be analysed using the iterated map of the phase differences. The return map reads ψ n+1 = ψ n + Z odd (ψ n ). Extended Data Fig. 9c uses iterated maps to show how a PRC with a negative slope at φ = 0.5 leads to a stable fixpoint in antiphase (left). An example with a positive slope at φ = 0.5 is shown in Extended Data Fig. 9b (bottom) and Extended Data Fig. 9c (right) and shows in-phase synchronization.
In conclusion, reciprocally, yet weakly coupled neurons, such as those connected through non-rectifying gap junctions, require an asymmetric coupling function (that is, with a substantial odd part) to show stable phase locking. For antiphase synchronization to be stable and the network in equation (2) to show frustration, the odd part of the coupling function must have a stable fixpoint at phase 0.5. In the low-frequency limit, this can be achieved with a PRC that has a negative slope at φ = 0.5, such as the PRC at the SNL point.
Calculating model coupling coefficients. Coupling coefficients for the model were calculated analytically. For this, we looked at pairs of gap-junction-coupled neurons (not the whole network of five neurons) and it was assumed that the voltage is sufficiently negative for all voltage-gated ion channels to be closed. In this case, the relationship between coupling coefficient (CC) and gap junction strength (g gap ) reads This relationship was verified in simulations that mimicked the experimental approach to determine coupling coefficients 29 . Specifically, two neurons coupled through gap junctions received steady current inputs below the firing threshold. The first neuron, n 1 , received a baseline current I base = −150 pA, the other neuron, n 2 , received an additional perturbation current I base + ΔI (ΔΙ from −5 pA to 5 pA). A third, uncoupled neuron n 3 was simulated with I base input. Simulations were run until all neurons reached steady state, resulting in steady state voltages v 1 , v 2 and v 3 . Then, the coupling coefficients were calculated as CC = . Here, a low I base ensures that the resulting CC is steady over the range of ΔΙ currents and matches the analytical results.
The splayness index. When spike times of all N neurons in a network are available, splayness can be quantified by comparing the neuron's phase differences to those that would arise in a perfect splay state.
To this end, we first discretized the time from the first spike of the last neuron that started spiking until the last spike of the first neuron that stopped spiking with a time step of 1 ms. If there were breaks in the recording (where the fly stopped and started flying again), those times were omitted. Again, the breaks were defined from the last spike of the first neuron that stopped spiking until the first spike of the last neuron that started spiking. This way, all phases of all N neurons were sampled with a 1 ms time step at K time points τ k = k ms for k = 0, 1, ..., K − 1.
For each neuron i, its phases were interpolated between its spikes t n,i as: To calculate the splayness, these phase differences were compared to the phase differences of the most splayed state (splay state) and to those of the least splayed state (sync state): k k With this definition, a perfectly splayed network yields s = 1 and a perfectly in-phase synchronized network yields s = 0.
Analysing pairs of neurons (synchronization index and phase histogram dips). Many experimental recordings include only two of the five MNs (Fig. 3d,h). In these cases, the splayness index cannot be used. Here we relied on a pairwise synchronization measure instead. Specifically, we interpolated the phases of the two neurons in question according to equation (5). Then, a classical Kuramoto order parameter was evaluated: A low synchronization index does not necessarily imply a splay state, as there are other desynchronized network states with equally low synchronization values. We therefore investigated alternative ways to gain information about splayness in the full MN1-5 network when in vivo recordings only comprised two neurons. The shape of pairwise phase histograms provides such information. Specifically, we observed that networks with high splayness often display a dip in pairwise histograms around zero. To test whether this dip is a good indicator for splayness, additional simulations were performed. With increasing noise strength, the dip around zero decreases, but it is more persistent than other histogram characteristics of splayed-out networks (Extended Data Fig. 6a). Furthermore, there was a strong negative correlation between the fraction of the MN4-MN5 pairwise histogram within (−0.1, 0.1) and the full-network splayness in experimental recordings for which all five neurons were available (Extended Data Fig. 6d). Together, this led us to conclude that the dip around zero is an appropriate indicator for splayness.

Comparing synchronization indices.
To compare synchronization indices (Fig. 3d,h), Mann-Whitney U-tests were performed using scipy. stats of SciPy v.1.7.1. Tests were uncorrected because only a few comparisons based on our hypotheses were performed.

Reporting summary
Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.

Data availability
The electrophysiology and imaging data have been deposited at the Zenodo open data repository (https://doi.org/10.5281/ zenodo.7737730). Computational modelling data and additional