Introduction

It has long been known that the spinal cord is capable of generating locomotor movements even in the absence of detailed instructive input from the brain [1]. The core dynamical features of spinal locomotor circuitry are left-right alternation and rostro-caudal propagation of cyclical activity over a range of frequencies. These core features are most obvious during swimming in aquatic vertebrates [2], but they are also apparent during limbed locomotion [3, 4, 5, 6, 7, 8].

Despite decades of work, the mechanisms by which rhythm is generated for coordinated locomotion are not yet fully understood. One line of research has focused on the intrinsic properties of individual excitatory spinal interneurons, whose mutual interactions serve to orchestrate rhythmic oscillations among downstream interneurons and motor neurons [9, 10, 11]. This perspective is supported by the observation of rhythmic bursting of excitatory interneurons in response to tonic input [12, 13], but it is potentially challenged by the observation that targeted disruptions of sources of phasic inhibition interfere with rhythmogenesis [14, 15, 16, 13]. Another line of work has suggested that rhythmic oscillations may be generated as a network-level phenomenon via random recurrent connectivity, without requiring specialized single-cell properties [17]. While this model is supported by the observation of “rotational” population dynamics rather than purely biphasic dynamics [17], the possible roles of particular cell types and their characteristic connectivity patterns are not addressed. In addition, the observation of recruitment of speed-selective interneuron populations [18, 11] and the functional roles of cell-type-specific intersegmental connectivity within spinal cord [19] have not yet been fully accounted for in either the cellular- or the network-level perspective.

Here, we construct a hierarchy of increasingly detailed models in which spinal central pattern generation occurs in a distributed manner, and the essential ingredients for producing coordinated locomotion are cell-type- and speed-specific connectivity motifs rather than specialized cellular properties or random connectivity. In all of these models, we represent individual neurons in a highly simplified way, so the dynamical properties of the network are entirely driven by emergent population dynamics due to connectivity patterns among units. First, we show that letting connectivity be determined by intersegmental phase relationships in a single-population, inhibition-dominated network is sufficient to account for left-right alternation and rostro-caudal propagation. Next, we show that dividing our units into fast- and slow-preferring populations is sufficient to account for variable-frequency control with constant phase lag via speed-dependent recruitment. Finally, in our most detailed model we show that the strength and modularity of recurrent excitation facilitates faster locomotion, but that there is a trade-off between maximum speed and controllability at intermediate speeds.

Together, these results point to an updated model of the spinal locomotor network as a distributed pattern generator, in which rhythm generation and other locomotion features emerge from collective population dynamics via patterned connectivity and speed-dependent recruitment.

Results

To investigate the degree to which population-specific connectivity patterns might account for the phenomenology of the spinal network described above, we developed a family of models at different levels of detail. In order to focus our investigation on connectivity patterns and the emergent population dynamics that result from them, we modeled the dynamics of individual neurons in a highly simplified way:

where r_i(t) is interpreted as the firing rate of unit i at time t, τ_i is the membrane time constant, D_i is the tonic drive, W_ij are recurrent weights from other units in the circuit, and [.]₊ denotes rectification. The axonal time delays Δ_ij are proportional to the number of segments between units i and j.

The three subsections that follow describe a hierarchy of increasingly detailed models, all of which are described by this basic dynamical equation. This approach enables us to focus mainly on basic computational properties within the relatively abstract, high-level model, whereas the lower-level models focus on implementation of distinct cell-type populations for obtaining additional features.

Phase relationships determine connectivity in a distributed pattern generator

Recent work mapping cell-type-specific connectivity patterns in zebrafish [20, 21, 19], mouse [22, 23], and other organisms has revealed a substantial presence of both excitatory and inhibitory connectivity traversing multiple segments. This is somewhat difficult to interpret in light of classical models of the spinal locomotor circuit, in which dynamical single-cell properties such as bursting or synaptic fatigue generate oscillations within each segment, and these oscillators are then chained together with descending excitation. We hypothesized that spatially patterned connectivity alone may be sufficient to drive coordinated locomotion even in the absence of dynamical single-cell properties. To test this, we constructed a high-level model in which dynamical single-cell properties are absent, and the basic features of left-right alternation and segment-to-segment propagation are driven by connectivity patterns. We took inspiration from “moving bump” models of brain circuits, for example in models of head-direction circuits, which have shown that inhibition-dominated, asymmetric connectivity can give rise to activation sequences [24, 25, 26]. Such models are also motivated by recent work in zebrafish demonstrating the importance of strong phasic inhibition for generating rhythmic locomotor activity [13]. This high-level model featured a single, homogeneous population of units described by Equation (1), where all units have identical membrane time constant τ and common excitatory input drive D.

The connections W_ij between units in this model were all set to be negative or zero, so that units in the model receive excitation only from the tonic drive, not from each other. Previous related work on threshold-linear recurrent networks with binary inhibitory weights and tonic input has shown that, despite their relative simplicity, such networks are capable of supporting rich dynamics [27]. Due to the (piecewise) linearity of the activation function, the magnitude of the nonzero weights has no effect on the dynamics, so we set W_ij = −1.

Given these simplifying choices, the nontrivial question in designing this model is which pairs of units should inhibit one another. The key principle that this model led us to is that the connectivity between pairs of units should encode the desired phase relationships between those units. A version of this idea is already present in classical models, where left and right populations of neurons within each segment laterally inhibit one another, ensuring that they fire exactly out of phase. Here, we apply this same principle to model the longitudinal connectivity of the spinal locomotor network. According to this principle, a mid-body unit should inhibit contralateral units in the same segment and in nearby segments, as well as inhibiting ipsilateral units more distally (Figure 1A). To make this precise, we can summarize the desired phase relationships for a pair of units i and j as follows:

Figure 1.

Where N is the number of segments. In this equation, s_i denotes the segment number of unit i, and a_i = 0 if unit i is in the left hemisegment or 1 if it is in the right hemisegment.

In order to realize these desired phase relationships, we set the connection strength from neuron j to neuron i as W_ij = f (ϕ_ij), where ϕ_ij ∈ [0, 1) is calculated from Equation (2), and

This equation describes local contralateral inhibition together with more-distal ipsilateral inhibition. Because ϕ_ij = 0.5 for two units that are out of phase, it is necessary to choose ϕ_l < 0.5 and ϕ_u > 0.5. Further, in order for propagation to proceed head-to-tail rather than tail-to-head, it is necessary for the connectivity to be asymmetric, with a relatively larger window of disinhibition downstream and a relatively smaller window of disinhibition upstream. Hence, we chose ϕ_l = 0.3 and ϕ_u = 0.8 to define the window of inhibition. Finally, we can note that the most-distal ipsilateral projections described by Equation (3) (of length ≳ N/2) are not actually necessary, since the units targeted by such inhibition already receive contralateral inhibition from nearby active units on the opposite side. Indeed, such long-range projections are not typically found in the zebrafish [28, 29, 30]. Hence, in this model and in the others presented below, all connections between segments more than 13 segments apart were set to 0 (Figure 1B). This truncation has essentially no effect on the dynamics, since the units that it would target are already inhibited by contralateral inhibition originating from nearby segments.

Driving this network with tonic input, and setting the axonal delay Δ_ij = 0 for simplicity, we observed the key features of coordinated locomotion: phasic bumps of activity in each segment, with smooth propagation of the bump from head to tail and strict alternation between left and right within each segment (Figure 1C). While this high-level model is a highly simplified and abstract representation of the spinal locomotor network, it illustrates several key ideas that may be relevant for characterizing locomotor circuitry in biological organisms and will continue to be present in the more elaborate models developed in the following subsections. First, rhythm can be generated by the network as a whole, rather than by dynamical single-neuron properties or single-segment oscillations. Second, the network’s recurrent connectivity is dominated by inhibition, illustrating that inhibition is capable of and perhaps necessary for sculpting the dynamics. Finally, the pattern of inhibitory connectivity is determined by the desired phase relationships between pairs of units. In the following sections, we will build upon this model by introducing multiple interneuron populations and allowing for more-heterogeneous connectivity between units.

Fast and slow speed modules implement frequency control via recruitment

While the single-population model described above was able to implement coordinated locomotion, it did so at a single, fixed frequency. This frequency was set by the only timescale in the model: the membrane time constant. Inspired by observations of speed-module structure in spinal circuitry in zebrafish [31] and mouse [32], whereby interneurons are recruited at either fast or slow locomotion speeds but not both, we made a minimal change to our initial model that introduces an additional timescale, replacing each unit by two units: a “fast” unit with a 1-ms membrane time constant, and a “slow” unit with a 10-ms time constant, as observed experimentally [33]. Trivially, in the absence of connectivity between the fast and slow modules, this model amounts to two copies of our earlier model and is capable of operating at two speeds: a fast speed if tonic input is provided only to the fast-preferring units, and a slow speed if tonic input is provided only to the slow-preferring units [34]. In each case, the period of left-right alternation and head-to-tail propagation would be roughly proportional to the corresponding membrane time constant (though perhaps not perfectly proportional if axonal delays are included).

In order for the two parts of the network to oscillate at a single global frequency and to interpolate continuously between slow and fast speeds, connectivity between the two modules is required. To begin, we adopted the simplest choice and made the connectivity weights the same for between-module projections and within-module projections, with the same spatial connectivity patterns for between- and within-module projections (Figure 2A-B). We found that the circuit exhibited coordinated locomotion characterized by head-to-tail propagation and left-right alternation at a single global frequency that was intermediate between the characteristic frequencies of the fast and slow modules. Further, the locomotion frequency could be controlled in a graded manner by driving the two modules with different inputs. Locomotion was fast when the fast units received most of the tonic drive, and it was slow when the slow units received most of the tonic drive, with smooth interpolation of intermediate frequencies as the ratio of tonic drives was varied (Figure 2C).

Figure 2.

Analyzing the model in greater detail revealed that the relative levels of tonic drive to the two populations strongly modulated the frequency (Figure 3B, solid line; Figure 3C), but that there was essentially no effect on frequency if the two drives were co-modulated in the same way (Figure 3B, dash-dot line; Figure 3D). We next computed the mean peak amplitudes of the fast and slow populations individually, as well as the amplitude averaged over the two populations (Figure 3A), which, as will be made clear in the following subsection when excitatory units are introduced, we take as a rough proxy for local excitatory activity and hence the amplitude of physical movements. In this case, we found that the average amplitude exhibited very little variation if the two drives were varied in opposite ways (Figure 3A, solid line; Figure 3C), but that the amplitude varied strongly if the two drives were co-varied in the same way (Figure 3A, dash-dot line; Figure 3D).

Figure 3.

Together, these results show that, by varying the drives to fast and slow units either together or differentially, a model with fast and slow populations can exhibit independent control of locomotion frequency and amplitude, as observed in zebrafish [35]. This is made possible by selective recruitment of fast and slow speed modules, such that the fast population is active at faster locomotion frequencies, while the slow population is active at slower locomotion frequencies, with a smooth crossover between these two regimes (Figure 3E). Such recruitment has been observed experimentally [31] but has not been accounted for by previous models. In Figure S1, we further show that the disengagement of the slow module as frequency increases can occur due to inhibition from the increasingly active fast module, without requiring changes in the tonic drive to the slow population, which is also observed experimentally [31].

To establish that the model exhibits coordinated locomotion across all frequencies, we additionally computed the left-right phase difference between the pairs of units on either side of each segment, finding that the phase differed by half of a period across all locomotion frequencies (Figure 3F). We then computed the phase difference between pairs of units in adjacent segments, finding that this phase was approximately equal to 1/N across all locomotion frequencies, where N is the number of segments (Figure 3G), so that the length of the spinal cord exhibits approximately one period of oscillation. This constant phase lag relationship is a core feature of locomotion that has been observed during swimming in a variety of animals [36]. While it does not appear by default in models that chain individually oscillating segments together via longitudinal excitation, constant phase lag emerges in our model as a natural consequence of the fact that the desired phase relationships between segments are built into the circuitry via long-range inhibitory projections. Because the phase relationships are determined by the connectivity, which is the same at all speeds, the phase relationships do not depend on locomotion frequency.

These results show that, by including fast and slow subpopulations and coupling these sub-networks together, the model succeeds in producing coordinated locomotion across a broad range of frequencies, with independent control of frequency and amplitude, via selective recruitment of the two subpopulations. While this model succeeds in producing much of the phenomenology of the spinal locomotor network, it does not yet fully address the diversity of excitatory and inhibitory cell types that are known to play a role in the biological circuitry [19]. We turn our attention to this question in the following section.

Excitatory and inhibitory cell types

We incorporated excitatory and distinct inhibitory populations by replacing each unit in the above model with four units, leading to an eight-population model (since each of the four types has fast and slow subtypes) (Figure 4A). One of these populations (corresponding to V2a interneurons in the zebrafish) consisted of excitatory units with descending, ipsilateral projections. The other three populations consisted of inhibitory units, essentially breaking the inhibitory population from the earlier model into three populations that have identical activity (since they all receive the same inputs) but differ in their projection targets. One population of inhibitory units (corresponding to V1 interneurons in the zebrafish) had ascending, ipsilateral projections; another (corresponding to V2b interneurons in the zebrafish) had descending, ipsilateral projections; and a third (corresponding to dI6 and V0d interneurons in the zebrafish) had contralateral projections. The tonic drive was provided equally to all units.

Figure 4.

The outgoing spatial connectivity of each of these cell types is illustrated in Figure 4B. Each cell type projects equally to all of the units within each segment that it targets, so that all units within each segment receive the same inputs. As in the two-population model above, we set the spatial connectivity patterns for the inhibitory units according to the desired phase relationships between units, with short-range contralateral inhibition and intermediate-range ipsilateral inhibition. For the excitatory units, we assumed that the projections are descending only in order to facilitate head-to-tail propagation. All units within each speed module had the same membrane time constant and axonal conduction velocity.

Before analyzing the full model in detail, we decoupled the two speed modules from one another and began by studying the effects of various single-cell and cell-type-specific connectivity properties on the characteristic oscillation frequency of an individual speed module. Unsurprisingly, the locomotion frequency depended strongly on the membrane time constants, with smaller values of these parameters leading to faster frequencies (Figure 5A-B). Because the characteristic frequencies of the two speed modules set the upper and lower limits of locomotion frequency once the modules are coupled together in the full model, it is likely advantageous for an organism to have values of these parameters that differ strongly in fast- and slow-preferring neurons. This agrees with observations from the zebrafish, where the membrane time constants and axonal delays differ for fast- and slow-preferring excitatory and inhibitory interneurons [37, 33]. For our subsequent simulations, we fixed these parameters for the fast and slow modules at the experimentally determined values indicated in Figure 5A-B [37, 33].

Figure 5.

We next investigated the effect of connectivity properties on locomotion frequency for the decoupled speed module. We found that the frequency was modulated by more than a factor of two as the excitatory projection distances were varied (Figure 5C). This agrees with observations from zebrafish, where the extents of intersegmental projections have been shown to differ for fast- and slow-preferring excitatory interneurons, with fast-preferring V2a interneurons projecting more distally than slow-preferring V2a neurons [20]. For our subsequent simulations, we fixed these parameters for the fast and slow modules at the experimentally determined values illustrated in Figure 4B [20].

Having shown that the spatial extent of excitatory projections has a strong effect on locomotion frequency, we next asked whether varying connectivity properties would also modulate the range of possible frequencies in the full model with two coupled speed modules. Varying the global strength of excitatory projections had a strong effect on the range of possible frequencies, with stronger excitation facilitating faster locomotion (Figure 6A). In particular, whereas the purely inhibitory model with experimentally determined membrane time constants and axonal conduction velocities realizes a maximum frequency much lower than that observed in larval zebrafish (approximately 20 Hz, black line in Figure 6A), the inclusion of excitatory interneurons facilitates maximum frequencies of over 50 Hz, which is in the same range as peak swim speeds in the zebrafish [13]. Thus, while excitatory interneurons are not necessary for producing coordinated locomotion in our model, they do facilitate faster locomotion, suggesting that this may be a fundamental role for feedforward excitation in the spinal network.

Figure 6.

Given that connectivity within and between speed modules has been shown in zebrafish to be modular, with stronger projections within modules than between modules [12], we asked what would be the effect of varying modularity in the model. We defined modularity as the difference between intra-vs. inter-module connection strength divided by the sum of these quantities, such that modularity of 1 corresponds to fully decoupled modules, while modularity of 0 corresponds to identical connection strengths within vs. between modules.

Varying the modularity of all four populations together, we found that there was essentially no effect on the maximum or minimum possible frequencies. Further, the model lost the ability to produce locomotion at intermediate frequencies as modularity was increased (Figure 6B). However, when we varied modularity among only the excitatory or only the inhibitory populations, we observed much more significant changes in the maximum frequency (Figure 6C-D). These changes occurred in opposite directions, suggesting that the lack of an observed change in frequency range when both types of modularity were varied together (Figure 6B) was due to cancellation between these two effects.

Together, these results show that the strength of feedforward excitation and the modularity of excitatory connectivity have a strong effect on the range of possible locomotion frequencies. There is a trade-off, however, in that the model loses the ability to smoothly interpolate between fast and slow frequencies in cases where the excitatory connectivity becomes too strong or too modular (Figure 6A,C). This requirement that excitation not be too strong is in accord with experimental observations from zebrafish [13], which have shown that peak excitatory post-synaptic currents are much weaker than peak inhibitory post-synaptic currents in V2a interneurons, consistent with the possibility that excitation may be globally weaker than inhibition in the spinal circuitry. Further, the fact that the model exhibits a frequency range similar to that of the zebrafish for parameters that are close to the critical values where smooth frequency control becomes impossible suggests that the spinal locomotor circuit faces a trade-off between speed and controllability, and that its excitatory connectivity may be configured in a way that optimizes this trade-off.

Having established the roles played by single-cell and connectivity properties of different cell-types in the eight-population model, we fixed these parameters and analyzed the behavior of the model over the range of possible tonic drives to fast and slow populations (Supplemental Figure S2). Similar to the two-population model, the eight-population with coupled fast and slow speed modules model exhibited head-to-tail propagation with constant phase lag (Supplemental Figure S2G), left-right alternation (Supplemental Figure S2F), independent control of speed and amplitude (Supplemental Figure S2C-D), and frequency-dependent recruitment of fast and slow populations (Supplemental Figure S2E).

We next investigated the effects of perturbing the model by partially ablating (i.e. attenuating the outgoing activity of) each interneuron population (Figure 7). At all locomotion speeds, we found that ablating excitatory units decreased locomotion frequency. This is in agreement with experiments in zebrafish, where ablation of excitatory V2a interneurons had the same effect [38]. Further, we found that ablating inhibitory units with ascending ipsilateral projections decreased locomotor frequency, while ablating inhibitory units with descending projections increased locomo-tor frequency across all locomotion speeds. This is also in agreement with experiments in zebrafish, where ablation of inhibitory V1 interneurons slowed swimming [39], while ablation of V2b interneurons led to faster swimming [29, 21]. Finally, we found that ablating the contralaterally projecting inhibitory units led to a modest increase in frequency, but that coordinated locomotion was lost when the degree of ablation became too great. The impact on frequency was most obvious at fast speeds, with a more modest impact at slow speeds. This is consistent with recent experiments in zebrafish, which found the impact of attenuating contralateral inhibitory projections from dI6 neurons on coordination was most obvious at fast speeds [13]. Similar results were found in Xenopus, where silencing contralaterally projecting inhibitory interneurons can eliminate rhythm generation [15] or leads to an increase in swim frequency [40].

Figure 7.

In addition to testing the effect of ablations on locomotion frequency, we also observed their effects on locomotion amplitude, which we define here as the average amplitude across all units in the network (Figure 7B). In contrast to the effects on frequency, here we observed that the effects were more selective at particular locomotor frequencies. Specifically, the amplitude increased strongly upon ablating inhibitory units with ascending ipsilateral projections during low-frequency locomotion, while the amplitude decreased strongly upon ablating inhibitory units with descending ipsilateral projections during low-frequency locomotion. Ablating the excitatory or contralateral inhibitory units, in contrast, did not have a strong effect on amplitude.

Together, these results show that, where comparisons with experimental data are possible, perturbations to our model lead to effects on locomotion frequency that generally agree with experimental observations. This agreement provides support for the possibility that the basic mechanisms underlying variable-frequency locomotion in our model—namely cell-type-specific connectivity patterns and speed-module recruitment—may also be at play in the spinal locomotor network. To our knowledge, the effects of ablations of interneuron populations on locomotion amplitude, such as those shown in Figure 7B, have not been systematically tested. Testing these predictions with future experiments would provide a further test of our model and the mechanisms that underlie it.

Discussion

In this study, we began by postulating that cell-type-specific connectivity alone could be sufficient for producing the main phenomenological features of the spinal locomotor circuit, without requiring dynamical single-cell properties. We found that coordinated locomotion could be achieved in an inhibition-dominated network in which connectivity is determined by desired phase relationships and variable-speed control is implemented by recruitment of frequency-selective populations. Further, while structured excitatory connections were not necessary for producing coordinated locomotion or frequency control, they were useful for increasing peak locomotor frequency, albeit at the cost of losing some control at intermediate frequencies. Together, this family of models shows that network-level interactions are sufficient to generate coordinated, variable-speed locomotion. It further provides new interpretations of intersegmental excitatory and inhibitory connectivity, as well as the basic, recruitment-based mechanism of speed control.

A main conclusion of our models is the importance of intersegmental, inhibition-dominated connections for achieving coordinated locomotion, where patterns of ipsilateral and contralateral inhibition are established by desired phase relationships. Similarly, a very recent network-level model based on the mouse locomotor circuit has proposed that these are key features for obtaining co-ordinated locomotion in that context as well [41]. The congruence of these results suggests that long-range, cell-type-specific connectivity patterns may be primary drivers of locomotor dynamics across species and fit cross-species observations of spinal interneuron diversity based on conserved molecular markers. While these results suggest that intrinsic cellular mechanisms are not necessary to generate coordinated locomotor rhythms, they are certainly sufficient in certain circumstances, as demonstrated by previous lesion and pharmacological studies [9].

Another conclusion of our models is that both projection strength and modularity among fast and slow excitatory units increase the maximum possible frequency, but at the cost of losing some control at intermediate frequencies. Recent studies in zebrafish suggest the possiblity that the spinal circuit may overcome this limitation by having multiple subtypes of excitatory interneurons with different degrees of modularity. Specifically, V2a neurons with descending-only axons exhibit a greater degree of speed-dependent recruitment (suggesting a higher level of modularity), while those with bifurcating axons exhibit less [13]. Moreover, V2a neurons with bifurcating axons fire more reliably, compared to more sparsely firing descending-only V2a neurons, and they form stronger connections to motor neurons [37, 13], consistent with a hierarchical organization distinguishing interneuron rhythmogenesis from motor neuron recruitment. Including motor neurons and these distinct subtypes into a model and testing their effects will be an interesting direction for future work.

Given that the models that we have presented favor simplicity over realism, attempting to capture as much phenomenology as possible with a minimal number of tunable parameters, they are undoubtedly missing features that may be important for describing more-nuanced functional aspects of locomotion in aquatic vertebrates. For instance, we did not include in our model variations in intersegmental projection distances for dI6 and V0d neurons related to speed [42], nor did we include excitatory commissural interneuron classes, including V0v [43] and V3 interneurons [44]. Moreover, other recent models of the locomotor circuit in zebrafish have highlighted the importance of electrical synapses for rapidly initiating swim bouts and characterizing early stages of development [45, 46, 47]. Including these components would be a worthwhile extension of the models presented here.

A somewhat unique aspect of our approach has been to develop a hierarchy of models to describe the same neural circuit at varying levels of detail. (See also the related approach in Ref. [46], which developed a series of models describing the same circuit at different stages in development.) Having such a family of models that fit within the same modeling framework enabled us to (i) use the higher-level models to better motivate the choices made in our lower-level models and (ii) draw connections between models at different levels to gain additional insight into the functional roles of particular populations. Following this approach, we thus obtained greater insight into the neural mechanisms underlying behavior than would be possible from any one model individually. We expect that this general approach could be useful more broadly for characterizing neural circuits and their relation to behavior.

Methods

Simulation details

The multi-segment model has 30 segments and 2 sides (left and right), for a total of 60 hemisegments. Each hemi-segment contains one unit corresponding to each neuron type. In the 1-population model, there is only one inhibitory neuron type; in the 2-population model, there are both fast and slow types of inhibitory neurons; in the 8-population model there is a fast and slow type of each of excitatory, ascending ipsilateral inhibitory, descending ipsilateral inhibitory, and contralateral inhibitory neuron type.

The state of each unit is described by its firing rate as a function of time. The time course is calculated by Equation 1, where h_i is the activity of unit i, and [·]₊ = ReLU() is the rectified-linear function. The membrane time constants are set to 1ms (10ms) for the fast (slow) module, or 1ms in models with no speed modules.

The numerical simulation was performed using Euler’s method with timestep equal to 0.1ms. For all simulations considered in this work, the simulation was run for 600 ms (or 6000 timesteps). The firing rate of each unit is set to a small random rate sampled uniformly from [0, 0.01] at t = 0 in order to break symmetry.

For the 1-population model there were no synaptic time delays included. In the 2- and 8-population models, synaptic time delays are equal to:

Where τ_delay is the base delay amount, and Δ_S is the distance (measured in number of segments) between the pair of units which the connection is between (Δ_S=0 when the units are in the same segment). For the 2- and 8-population models, τ_delay is set to 0.2 ms for the fast population and 0.5 ms for the slow population, matching roughly the experimentally measured axonal conduction velocities detailed in Ref. [37]. Note that the synaptic time delays are determined by the identity of the source unit, not the target unit.

Units in the model are driven by a constant tonic input. This drive is varied separately for the fast and slow populations. In addition, each unit receives recurrent input from the other units according to the interneuron connectivity W_ij, which represents the strength of the connection from unit j to i. The connectivity matrices, W_ij, are shown in Figures 1B, 2B, and 4B.

The base value of inhibitory projections is set to −0.5. The base strength for the excitatory connections is set to 0.5f_E, where the multiplicative factor f_E takes a value between 0 and 1 and determines the strength of excitation relative to inhibition. Except when this parameter is explicitly varied (see Figure 6), we set f_E = 0.4.

In the 2-population and 8-population models, we introduce the speed mixing factor, f_sm. The connection strengths for all fast-to-fast and slow-to-slow connections were multiplied by (1 − f_sm), whereas all fast-to-slow and slow-to-fast connection strengths were multiplied by f_sm. This generates a modularity of m = 1 − 2f_sm. In the case of the 8-population model, we can apply a global speed mixing factor f_sm or a separate speed mixing factor for the excitatory population f_sm,E and for the inhibitory population f_sm,I. The results of varying these parameters are shown in Figure 6.

Ablation was introduced as an overall factor, f_A, that multiplied the connection strengths of all connections sourced from the ablated population. This then implies that

Analysis methods

The output of the simulation is a collection of time series giving the firing rate of each unit as a function of time. By inspection, it was clear that the time series settled into a sensible oscillation after an initialization period. To ensure the observations corresponded to the stable oscillatory mode, all analysis was performed on the times series from t=100ms to t=600ms (i.e. we cut out the first 100ms).

For each unit’s firing rate time series, r(t), the amplitude is defined as the difference between the maximum and minimum values of the time series. That is

The amplitude of the simulation is then defined as the mean amplitude across all units (or, if amplitude is reported for a particular population, across all units within that population). Errors are given by the standard deviation in amplitude across all units.

To extract the frequency of the time series r(t), we calculated the period of each time series. To find the period, we considered the autocorrelation spectrum defined by:

Here r(t) = 0 whenever t is outside the domain t ∈ [100ms, 600ms]. In all cases the maximum autocorrelation occurs at k = 0. In cases with a single dominant frequency as k increases away from zero, the autocorrelation drops to global minimum (which we define as k = k_min) then rises to a local maximum at a value of k equal to the period (see Figure S3). This is then followed by a series of local minima and maxima of lesser size. Crucially, this local maximum can be isolated as the global maximum if we only consider k > k_min. Hence, we define the period of r(t) by

and we define the frequency of r(t) by f = 1/T. From there frequencies of the different unit types were averaged within each hemi-segment weighted by the amplitude. This amplitude-averaged frequency was then used to compute an unweighted mean frequency and standard deviation across all hemi-segments.

We calculate the phase ϕ of r(t) as

Where arg(z) gives the argument (or phase) of the complex number z and is the Fourier coefficient of r at the measured frequency, ν. Notice that we have normalized our phases to be in the range (−1/2, 1/2] instead of (−π, π]. The phase difference Δϕ_ij between phases ϕ_i and ϕ_j is calculated with

to avoid erroneously large phase differences for phases near and . Each phase difference at each location is averaged via a standard amplitude-weighted mean. A circular mean and circular standard deviation are used to then average the amplitude averaged phase differences across all hemi-segments.

Failures of coherent oscillation

Two issues arise that can prevent the extraction of a single well-defined frequency from a completed simulation: (i) There are multiple dominant frequencies in the time series for some units, or (ii) the frequencies do not agree between the fast and slow populations.

In the case of having multiple dominant frequencies, the autocorrelation spectrum no longer follows the easily interpretable shape described above. Rather, the time series and autocorrelation appear like the example in Figure S3E and S3F. A key observation is that the global minimum coincides with the first local minimum only in cases with a single dominant frequency. Therefore, to systematically find those time series with multiple dominant frequencies, we compare the frequency found using the global minimum as k_min with the frequency found using the first local minimum as k_min.

In some cases (especially at high modularity), a well-defined frequency occurs within the fast population that differs significantly from a well-defined frequency within the slow population. This also demonstrates a failure for the system to oscillate coherently.

To systematically determine whether a simulation resulted in a well-defined frequency, we compare four frequency values:

1. The mean over the fast population of frequencies found using the global minimum

2. The mean over the slow population of frequencies found using the global minimum

3. The mean over the fast population of frequencies found using the first local minimum

4. The mean over the slow population of frequencies found using the first local minimum

Only if all values agree (to within a small tolerance) do we consider the frequency well-defined for the whole population.

Supplemental Figures

Supplemental Figure 1:

Supplemental Figure 2:

Supplemental Figure 3:

Acknowledgements

We are grateful for discussions with Martha Bagnall. This work was supported by the National Institutes of Health BRAIN Initiative (U01-NS136458).