Model network constructionModels of neuron and synapse

In the proposed neuron network models, simulation for each neuron is based on the leaky integrate-and-fire (LIF) model with conductance-based synapses described in the previous study (Su et al. 2017; Huang et al. 2019). We chose a spiking neuron model with biophysically realistic synapse models over the simpler firing-rate-based models because a spike-based model reveals several challenges posed by real nervous systems and can be used to explore more detailed neural mechanisms than the rate-based models do. For example, dynamical instability caused by sparse spike inputs cannot be investigated by rate-based models because their synaptic inputs are always continuous. Hence the spiking neuron model is a reasonable choice.

In the present models, the LIF model is given by

$$_\frac_}=-_(_-_)-__$$

(1)

where \(_\) is the membrane capacitance (= 0.1nF for all neurons), \(_=_/_\) is the membrane conductance and is set to make the membrane time constant \(_\) equal to 15 ms for each neuron, \(_\) (= − 70 mV) is the resting potential, \(_\) is the synaptic current elicited by spike input from neuron j. The synapses are conductance-based, and three types are modeled in the present study: glutamatergic (via NMDA receptor), cholinergic (Ach), and GABAergic (GABAA). The synaptic currents induced by cholinergic, AMPA, and GABAergic receptors are given by

and as for NMDA receptors, the synaptic current is

$$_=\frac__(_-_)}^]}^_}}$$

(3)

where \(_\) is the synaptic conductance (equivalent to synaptic weight), \(_\) is the gating variable, \(_\) is the reversal potential, which is set to 0 mV for all excitatory synapses (including NMDA, and Ach) and is set to -70 mV for the inhibitory GABAA synapses. [Mg2+] is the extracellular magnesium concentration. The gating variable \(_\) describes the activation of channels due to spike input and is given by

$$\frac_}=-\frac_}+_\delta (t-_^)$$

(4)

for GABAA and Ach receptors, and by

$$\frac_}=-\frac_}_}+_\delta (t-_^) , \frac_}=\alpha _(1-_-\frac_}_})$$

(5)

for NMDA receptors, where \(\tau\) is the time constant (5 ms, 20 ms, and 100 ms for GABAA, Ach, and NMDA receptors), \(\delta\) is the delta function, \(_^\) is the time of the k-th presynaptic spike from neuron j, \(\alpha\)(= 0.6332) is a factor for adjusting the increment of the NMDA gating variable, or receptor activation rate.

All the neuron and synaptic receptor-related parameters described in this section are fixed in the present models. The value of the membrane time constant (\(_\)=15 ms) is determined based on the electrophysiological study of the DM1 projection neurons of the fruit fly antennal lobe (Gouwens and Wilson 2009). Other parameters are adopted from the typical values used by other studies or by us in the neural network models of fruit flies or other species (Wang 2002; Lo and Wang 2006; Huang et al. 2019).

The network models

In the present study, we constructed two classes of the HD circuit models with one (R class) based on Su et al. 2017 and Han et al. 2021 and another (Delta class) based on Kakaria and de Bivort 2017 and Pisokas et al. 2020. Each class of models included three types of neurons, EPG, PEN, and R (for the R class only) or Δ7 (for the Delta class only). EPG neurons, or the compass neurons, support attractor dynamics and form localized activity, or an activity bump, that encodes the head direction (Clandinin and Giocomo 2015; Seelig and Jayaraman 2015; Turner-Evans and Jayaraman 2016; Kim et al. 2017). PEN neurons form the shifter circuits, and they can move the activity bump clockwise or counterclockwise when the body rotates in the absence of visual input, i.e., in darkness (Fig. 1a, right). The global inhibition, which is essential for the formation of activity bump, is provided by the GABAergic ring (R) neurons (in the R class models) or by the Δ7 neurons (in the Delta class models). Each neuron class (EPG, PEN, R, and Δ7) exhibits a specific innervation pattern at the neuropil level and can be further divided into several types (EPG-L1, EPG-L2, for example), which innervate different subregions of the neuropils. Based on the connectomic data, 1–5 neurons were identified for each type (Scheffer et al. 2020). Therefore, we included three identical neurons for each type in the present models.

Instead of using the exact circuit structures proposed in the original studies, we revised the circuits of each model based on the recently released connectome of the central complex (Scheffer et al. 2020; Hulse et al. 2021). Several changes were made: (1) In the original R class and Delta class models, the local recurrent excitation was hypothesized to be carried out by the EPG↔PEG loops. PEG is a set of excitatory neurons projecting from PB to EB and the Gall. However, the connectomic data indicate very weak PEG → EPG connections. Instead, there are strong self-recurrent connections within the EPG neuronal group. Therefore, we removed EPG↔PEG loops and replaced them with EPG↔EPG loops in both classes of models (Fig. 1b and c). (2) In the original Delta class models, each Δ7 neuron innervates all 18 glomeruli in the PB with two axonal and 16 dendritic connections. However, the connectomic data (Scheffer et al. 2020) do not support a broad dendritic innervation in PB. With a careful analysis (Supplemental Fig. 1), we determined that dendrites of each Δ7 neuron only innervate 6–8 PB glomeruli (Fig. 1d). (3) A recent study (Hulse et al. 2021) revealed a total of 18 EPG neuron types, EPG-L1 to L9 and EPG-R1 to R9. Except for EPG-L1 and R2, each of the remaining 16 EPG neurons provides inputs to one of the 16 PEN neurons and receives PEN inputs from neighboring PEN neurons, forming perfect clockwise or counterclockwise shifter circuits (Fig. 1e, left). In contrast, the EPG-L1 and R1 neurons break this pattern by not projecting to any PEN neuron (Fig. 1e, right). Because it is not clear whether these two “atypical” neurons participate in the head-direction function, we created two variants of the models, one with EPG-L1 and R1 neurons (named E18) and one without (E16) for each of the R and Delta model classes. Therefore, we investigated four basic models in the present study: R-E18, R-E16, Delta-E18, and Delta-E16 (Fig. 1f). The connection tables and innervation tables of the four models are shown in supplemental Figs. 2 and 3, respectively. We also investigated the Hybrid model, which combined R-E16 and Delta-E16.

Previous studies reported that the PEN class neuron consists of two subtypes: PEN_a and PEN_b (also named PEN1 and PEN2, respectively) (Green et al. 2017; Hulse et al. 2020). The functional differences between both subtypes are unclear. In the present study, we regarded the PEN neurons as the PEN_a subtype because PEN_a dendrites innervate EPG neurons directly in PB, while PEN_b neurons only contact EPG neurons indirectly via PEG neurons (Turner-Evans et al. 2020), which the present models do not need.

As mentioned earlier, each Δ7 neuron only innervates a subset of the PB glomeruli. However, the whole population of Δ7 neurons tiles the PB and forms a complete and homogenous coverage, except for three atypical Δ7 neurons, L4R6_R, L6R4_L, and L7R3_L. These three neurons do not follow the order of the innervation pattern exhibited by other Δ7 neurons (Supplemental Fig. 4). We have tested the Delta models by including the three atypical Δ7 neurons, but the models failed to work completely. Therefore, in the present study, we excluded L4R6_R, L6R4_L, and L7R3_L from the two Delta models (Delta-E18 and Delta-E16).

Note that although the network structure of each model was determined by the connectomic data, we kept the synaptic weight (wij in Eqs. 2 & 3) of each connection as a tuning parameter. The weight of each synapse in the models is a product of two variables: factor and base. The base value (K) defines the basic weight of synapses between each pair of neuron classes (PEN to EPG, for example) and is a tunable parameter in the present study (see Fig. 2c and supplemental Fig. 5 for tuning ranges). The factor is used to scale the weight values within each pair of neuron classes and is not tunable (Supplemental Fig. 2). For example, the factor is 2 for PEN-L3→EPG-R8 but is 1 for PEN-L3→PEG-R7, and all PEN to EPG connections share the same base value, which ranged from 2 to 25 when we swept through the parameter space and looked for usable parameters. The connectomic data provide information about the number of synapses formed between each pair of neurons (Scheffer et al. 2020; Hulse et al. 2020). Although it is tempting to use such information to indicate the synaptic weights, the synaptic numbers revealed in the connectomic data are highly variable, even within the same types of connections. For example, the mean synaptic numbers between EPG and PEN_a neurons vary from 1 (averaged over three EPG_L7 neurons → one PEN_a_L2 neuron) to 45 (averaged over two EPG_L1 neurons → one PEN_a_L2 neuron). All tested models failed if we simply set the synaptic weights proportional to these numbers.

Fig. 2

Robustness of the models. a The protocol for the robustness test (also shown in Fig. 1g). b Examples of one successful trial (left, R-E16) and three failed trials (right three, front left to right: R-E16, R-E16, Delta-E18). The successful trial produced a stable activity bump that tracked visual cue movement and body rotation, while the failed trials did not. c We swept through four-dimensional parameter space and marked the parameter sets that led to successful trials (orange and purple dots). The panel shows, for each model, a three-dimensional slice in the four-dimensional sweeping: (top left) R-E16, KR→EPG = 14, (top right) R-E18, KR→EPG = 14, (bottom left) Delta-E16, KDetla→EPG = 14, and (bottom right) Delta-E18, KDetla→EPG = 14. K is the weight base, which needs to be multiplied by the factor shown in supplemental Fig. 2 to become the synaptic weight. d Robustness (defined by the number of usable parameter sets) as a function of the Δ7/R → EPG synaptic strength for all models. The R models were more robust than the Delta models e The widths (FWHM) of the bumps in all four models. The red dashed line shows the experimental bump width data from (Kim et al. 2017). Due to the innervation patterns of the Δ7 neurons (Fig. 1d), the mean bump widths of the Delta models were significantly larger than the R models (t-test, *** = p < 0.001). An EB wedge is 22.5° (or 0.125π) wide. f Example trials of the R-E16 and Delta-E16 models with heterogeneity in the synaptic weights of the EPG → R/Delta and Delta/R → EPG connections

Simulation protocolsEstimation of calcium activity and the activity bump

In the present study, the outputs of the models are spike trains, which are different from the calcium activities observed in most experimental studies. To compare the modeled and observed EPG activity in each EB wedge, we followed the method used by Su et al. 2017. Specifically, for each EB wedge w, we calculated the population mean activity rw(t) of the EPG neurons by

$$_(t)=\frac_^_^\delta (t-\tau -^)^d\tau$$

(6)

where N is the number of EPG neurons projected to the wedge w, n is the number of total spikes generated by these N neurons, \(^\) is the time of the i-th spike, \(\delta\) is the delta function, and \(d\) (= 721.5 ms) is the exponential kernel indicating a 500 ms half-life, mimicking the hundred seconds of the half-life of the calcium indicator GCaMP6f (Chen et al. 2013). We further estimated the activity bump by performing Gaussian fitting to rw (t) over w for each timestep t. The resulting Gaussian function was treated as the representation of the activity bump. The peak position, height, and width (FWHM, Full width at half maximum) of the Gaussian function were used in the subsequent analyses.

Visual cue and body rotation

As demonstrated in previous studies (Seelig and Jayaraman 2015; Green et al. 2017), the activities of EPG neurons in EB and PEN neurons in PB encode the azimuthal position of visual stimuli. For simplicity, we mapped the 16 EB wedges to the 360° horizontal visual space as described in the previous study (Su et al. 2017). We activated the corresponding EB wedge to simulate the presentation of the visual cue. The present study used two types of stimulations, visual cue and body rotation, depending on the purpose of the tests. For the visual cue stimulation, we activated the corresponding EB wedges by stimulating the PEN neurons, which innervate the EB wedges. The activation is done with a spike rate of 50 Hz and maximum conductance of 2.1 ns on the cholinergic receptors. Note that the same effect can be achieved by stimulating the EPG neurons. However, both methods are equivalent from the modeling perspective due to the direct connections from PEN to EPG neurons. The second type of stimulation represents body rotation and is used to move the activity bump in the darkness without any visual input. This is done by activating unilateral PB via spike inputs (2,210 Hz, with a conductance of 0.3 ns on the NMDA receptors) to PEN L2-L9 or R2-R9 neurons, which constitute the shifter circuits that move the activity bump clockwise or counterclockwise, respectively.

Robustness test

In the robustness test, we swept through a wide range of synaptic weight bases (K) to determine the working parameter space of each model. The test included four types of synapses, KEPG→PEN, KPEN→EPG, KDelta→EPG (or KR→EPG), KEPG→Delta (or KEPG→R), and the bases ranges from 5 to 25 ns with a step of 1 for the first two types of synapses and 1–20 ns with a step of 1 for the last two. Therefore, a total of 176,400 sets of weights were tested for each model. We ran one trial for each set of weights, and each trial consisted of a 10 s period of visual cue stimulation followed by a 10 s period of body rotation stimulation. The cue moved 45 degrees per second in the visual cue period, while no visual cue was presented in the body rotation period. The 10 s body rotation period was further separated into a 5 s counterclockwise rotation (stimulating PEN R2-R9) followed by a 5 s clockwise rotation (stimulating PEN L2-L9) (Fig. 1gi). For a given parameter set, a trial was considered as failed if it met any of the following four conditions. First, the amplitude of the activity bump was smaller than 1 spike/s for more than 10 ms (diminished bump condition). Second, the width of the activity bump was larger than 360° for over 10 ms (spread bump condition). Third, the bump failed to move in the body rotation period (11-20 s) (immovable bump condition). Fourth, the Gaussian fitting failed to return any result for more than 5 ms (no bump condition).

Persistency test

The persistency test is separated into two sub-tests: the static and motion persistency tests.

The static persistency test examined how well the activity bump maintains a static location after the visual cue stimulus is turned off. The test started with a 1 s static visual cue stimulation followed by a 9 s total darkness (visual cue off) without any body rotation stimulation (Fig. 1gii). A trial was considered as successful if no “diminished bump” and “spread bump” conditions occurred (see Robustness Test above for the definition of the two conditions). We performed 1000 trials for each model with randomly selected usable parameter sets for each trial. For successful trials, we further analyzed how steady the bump was in the 9 s dark period by computing the standard deviation of the bump peak from the original visual cue location.

The motion persistency test used the same test protocol as the robustness test (Fig. 1giii). In the body rotation period, we expect that an HD model’s bump first shifts linearly in the counterclockwise direction for 5 s and then reverses the direction for another 5 s. To evaluate how persistent the bump exhibits this linear movement, we performed linear regression for the trajectory of the bump peak separately for the first 5 s and the last 5 s of the body rotation period. We used the mean \(^\) values of the two regressions as the measure of persistence. The mean \(^\) is between 0 and 1, with 0 representing no detectable bump movement while 1 representing perfect linear movement.

Speed test

The speed test consisted of only a visual cue period. In each trial, the visual cue moved for one round (360°) in the counterclockwise direction, followed by the clockwise direction for another run. (Fig. 1giv). To evaluate how fast the activity bump can follow the moving cue, we tested the following cue speeds: 0.25, 0.28, 0.312, 0.35, 0.42, 0.5, 0.625, 0.83, 1.25, and 2.5 \(\pi rad/s\). For the last two speeds, we allowed additional 3 and 7 rounds in each direction, respectively, so that the trials lasted long enough for the analysis. We randomly selected seven usable parameter sets for each model and simulated 1000 trials for each parameter set in each speed condition. We followed the criteria stated in Robustness Test except that the “diminished bump condition” is changed from 10 to 5 ms. We classified each trial as successful or failed and then calculated the success rate by dividing the number of successful trials by the total number of trials. Finally, we computed the mean success rate over the seven parameter sets for each speed condition and each model. The seven parameters are manually picked to be uniformly distributed in the parameter space to ensure that the results reflect the overall performance from an unbiased sampling of the useable parameters.

Dynamical characteristics test

In the dynamical characteristics test, we tested whether the activity bump “jumps” to the new location accordingly when the visual cue suddenly shifts from the original location to the new location. The bump jumping usually occurs when the two cue locations are far from each other (Kim et al. 2017). To this end, we selected two EB regions that were 180° apart, and the visual cue switched between the two regions twice in each trial (Fig. 1gv). To identify the trials that exhibited successful bump jumping, we checked whether a discontinuity of the active bump occurred immediately following each cue switch event (5 s ~ 5.5 s and 8 s ~ 8.5 s) by checking if “no bump condition” was detected during these two switching periods. The success rate is calculated by dividing the number of successful trials by the total number of trials.