A theoretical circuit for heading integration

We consider a population of ${\displaystyle N}$ ‘compass neurons’ with an activity that encodes the direction of the insect as an angular variable θ. We represent the activity of this population by a vector where each element corresponds to the activity of one neuron, ${\displaystyle \mathbf{a}(\theta )=[{\mathbf{a}}_1(\theta ),{\mathbf{a}}_2(\theta ),...,{\mathbf{a}}_N(\theta )]}$. We take ${\displaystyle N=8}$ neurons, consistent with data from many insect species which possess an eight-column organisation, with each column encoding a different direction (Honkanen et al., 2019; Stone et al., 2017; Pisokas et al., 2020).

Each neuron’s activity is updated depending on its current firing rate and the inputs it receives, both from other neurons in the circuit and externally. We formulate this update rule as follows:

where ${\displaystyle \mathbf{W}}$ is the circuit’s weight matrix, representing the connections between neurons; ${\displaystyle \varphi }$ is the neural activation function, that converts the total neural input into an output firing rate; and ${\displaystyle \mathbf{u}(t)}$ is the external input that encodes the insect’s angular velocity, via the PEN population of neurons (Green et al., 2017; Turner-Evans et al., 2017; Sayre et al., 2021; Pisokas et al., 2020; Turner-Evans et al., 2020; Hulse et al., 2021).

To simplify our derivations we allow the neural activity values to be both positive and negative, interpreting these values as being relative to a baseline neural firing rate. Similarly, we allow the weights to be both positive and negative, a common simplification in computational models (Li et al., 2023; Cornford et al., 2020; Kriegeskorte and Golan, 2019). This simplification will be addressed in section ‘Comparing the predicted circuit with biological data’ where we compare our model with experimental data.

Mathematical principles for neural heading integration

A circuit capable of performing heading integration must fulfil the requirements outlined in Table 1. The first two requirements allow us to establish the family of possible path integration circuits. We then use the principle of noise minimisation to determine which circuits perform best.

Table 1. Requirements for a heading integration circuit.

Circular topology	Rotational symmetry	Noise minimisation
The activity should have the same topology as the variable it is encoding to prevent discontinuities. To encode heading, the activity should have the topology of a 1D circle.	The heading integration circuit should work similarly, irrespective of the direction in which the insect travels. There should not be a bias for any direction.	The circuit should minimise the noise of the neural representation so the insect can navigate as precisely as possible.

Open in a new tab

Determining the optimal circuit

As we consider networks with ${\displaystyle N=8}$ neurons, consistent with multiple insect species (Honkanen et al., 2019; Stone et al., 2017; Pisokas et al., 2020), the possible single harmonic circuits are ${\displaystyle f=\{1,2,3,4,5,6,7\}}$, where the zeroth harmonic is discarded because it only represents the baseline neural activity. This gives the following activities and weights, derived from Equation 5 and Equation 14

The circuits for the first four harmonics are plotted in Figure 3, but not all of these circuits are valid. In particular, circuits with even spatial frequencies have multiple neurons with identical activity values because they share a common divisor with ${\displaystyle N=8}$, as explained in detail in section ‘Degenerate circuits’.

For ${\displaystyle f=4}$ there are only two unique activity values, ${\displaystyle {\mathbf{a}}_n=\pm \cos (\theta )\mathrm{\forall }n}$, so this circuit can only encode one dimension, not a circular topology (see Figure 3D). For ${\displaystyle f=2}$ and ${\displaystyle f=6}$ (not plotted), there are four unique activity values, allowing the angle to be properly encoded. However, these circuits are degenerate because they are composed of two independent subcircuits, each encoding one direction (see Figure 3B). Because the weights connecting the two subcircuits are all zero, the activity in each subcircuit is independent of the activity in the other, and so the overall activity cannot be constrained to have the required circular topology.

All circuits with odd harmonic frequencies are equivalent. For example, as shown in Figure 3, the connections in the circuit for ${\displaystyle f=3}$ are the same as those for the ${\displaystyle f=1}$ circuit after the neuron identities are permuted. As detailed in section ‘Equivalence under permutation’, the frequencies ${\displaystyle f=\{1,3,5,7\}}$ always give the same circuit because the odd frequencies are coprime with the number of neurons ${\displaystyle N=8}$.

Since the activities and weights of all the non-degenerate circuits ${\displaystyle f=\{1,3,5,7\}}$ are the same as the base harmonic ${\displaystyle f=1}$ up to a permutation, we choose the lowest harmonic ${\displaystyle f=1}$, which gives us the following activity and weights:

Comparing the predicted circuit with biological data

Our theory proposes an optimally noise-resistant circuit for heading integration, and its corresponding activity. The prediction that heading should be encoded as a sinusoidal activity bump is consistent with previous theoretical models (Touretzky et al., 1993; Wittmann and Schwegler, 1995; Hartmann and Wehner, 1995; Zhang, 1996; Vickerstaff and Di Paolo, 2005; Haferlach et al., 2007; Stone et al., 2017), as well as experimental evidence in both the locust and fruit fly (Heinze et al., 2009; Turner-Evans et al., 2017). We note, however, that data from the fruit fly shows a more concentrated activity bump than what would be expected from a perfect sinusoidal profile (Seelig and Jayaraman, 2015; Turner-Evans et al., 2017), and that calcium imaging (which was used to measure the activity) can introduce biases in the activity measurements (Siegle et al., 2021; Huang et al., 2021). Thus the sinusoidal activity we model is an approximation of the true biological process rather than a perfect description.

Importantly, our theory proposes that the optimally noise-resilient heading integration circuit should have synaptic weights that follow a sinusoidal pattern, even though such weights are not necessary for producing sinusoidal activity as discussed in section ‘Minimising noise propagation in the circuit’. This is the main prediction of our theory and we sought to validate it using connectivity data from the locust and fruit fly.

However, before we can directly compare our model with biological circuits, we must address a number of modelling simplifications:

• In the model we considered a population of 8 neurons, but in insects there are eight neural columns, each with several neuron types (EPG, PEG, PEN, Delta7). Of these we model only the EPG (known as the compass neurons) which encode the integrated head direction.

• The synaptic connections between neurons in our model can be both positive and negative, while biological neurons follow Dale’s law, meaning that a single neuron can either have only positive or only negative synapses.

• The neural activity in the model was centred around zero and could have both positive and negative activity (firing rates), while real neurons only have positive firing rates.

Therefore, we simplified the biological connectivity to produce an equivalent circuit that could be directly compared with our model prediction. The neural population in our theoretical model corresponds to the biological EPG neurons, as these encode the integrated heading. We considered the other 3 neuron types that are part of the compass circuit (PEG, PEN, and Delta7 – see Kakaria and de Bivort, 2017; Pisokas et al., 2020) as just implementing connections between EPG neurons in accordance with biological constraints. We counted the number of different paths between EPG neurons, accounting for the sign of the connections (whether the path passed through an inhibitory Delta7 neuron), and used the net path count as a proxy for connectivity strength. This process is explained in detail in section ‘Path counting’.

We then computed the average connectivity profile, i.e., how each neuron connected to its neighbours around the ring, and compared this profile to the closest fitting sinusoid. Because the neuron gains, absolute synaptic strength and biophysical properties of the neurons are unknown, the units of the net path count are not necessarily equivalent to our abstract connection strength. We therefore fit an arbitrarily scaled and shifted sinusoid to the connection counts:

where ${\displaystyle m-n}$ is the circular distance between 2 neurons while β and γ are constants that are fit to minimise the precision-weighted mean squared error compared with the experimental connectivity profile (see section ‘Fitting weights’).

Analysing the data from Pisokas et al., 2020, we consider the shortest excitatory and inhibitory paths between EPG (also known as CL1a) neurons in the locust, which have lengths of 2 and 3, respectively. There are no direct connections of length 1, and all the paths of length ≥4 must pass through the same neuron-type multiple times. This path counting analysis for the locust is shown in Figure 4A, and the procedure is detailed in section ‘Path counting’. We find that the connectivity profile between neurons for the locust is very close to sinusoidal in shape, supporting our theoretical prediction.

For the fruit fly, we used data from Scheffer et al., 2020, which provides synapse counts between pairs of neurons. These synapse counts were used as a proxy for connectivity strength, a view that has been validated in previous experiments (Liu et al., 2022; Barnes et al., 2022). After identifying the neurons and connections of interest, we grouped the neurons in eight columns following the logic presented in Pisokas et al., 2020, with detailed methodology explained in section ‘Data preprocessing’. We repeated the path counting analysis for the fruit fly with synapse count data (Figure 4B), and found that while the data are noisy, the connectivity profile fits a single sinusoid pattern reasonably well. However, the high variability in the synapse counts makes our hypothesis difficult to differentiate from alternative shapes (see section ‘Fitting weights’).

Taken together, this analysis shows that our theory is consistent with experimental data – using binary connectivity data for the desert locust and synapse count data for the fruit fly.

Learning rules and development

Having validated the connectivity of our theoretical circuit by comparing it to experimental data, we ask whether our circuit lends itself to biological development. Specifically, we show that even though our circuit requires precise connection strengths, this connectivity can be developed naturally by a Hebbian learning rule. Our approach follows from previous research which has shown that simple Hebbian learning rules can lead to the emergence of circular line attractors in large neural populations (Stringer et al., 2002), and that a head direction circuit can emerge from a predictive rule (Vafidis et al., 2022). In contrast to this work, we focus only on the self-sustaining nature of the heading integration circuit in insects and show that our proposed sinusoidal connectivity profile can emerge naturally.

Because the weight matrix is circulant, its eigenvalues are equal to the Fourier spectrum of its first row, which we derived to only have one nonzero value, for ${\displaystyle f=1}$. The network therefore only has a single eigenvalue and so projects activity into a single dimension. This operation is similar to classical dimensionality reduction methods such as principal component analysis, which can be implemented by Hebbian-like learning rules (Dayan and Abbott, 2001). We thus analyse the effects of incorporating a modified Oja’s rule into our model, a classical variant of Hebbian learning where the synaptic strength between 2 neurons grows when both neurons are active simultaneously, and the total synaptic strength is regularised to prevent exploding weight growth,

where ${\displaystyle {\mathbf{W}}_{n,m}}$ is the synaptic connection strength from neuron ${\displaystyle m}$ to neuron ${\displaystyle n}$, ${\displaystyle a_m}$ and ${\displaystyle a_n}$ are the pre- and post-synaptic activities, respectively, and ${\displaystyle \eta {\displaystyle \frac{d\theta }{dt}}}$ is an adaptive learning rate where the plasticity is proportional to the rotational speed, whereas in the classical Oja’s rule it would be constant.

For our analysis we assume that the insect faces all possible directions, and therefore that the neural activity goes around the full circle. We then integrate the weight updates over some long period of time,

where θ is the integration space. Applying the activity from Equation 17, we can find the fixed point of this update rule, when ${\displaystyle \mathrm{\Delta }{\mathbf{W}}_{n,m}=0}$:

Combined with Equation 7 and Equation 17, this result means that if there is sinusoidal activity in the network, the weights will naturally converge to the optimal sinusoidal values by way of our variant of Oja’s rule. This leads to the emergence of sinusoidal activity and weights from self-consistency: initial noisy sinusoidal activity will enforce weights that are close to a sinusoid, and those weights will filter out noise to make the activity even closer to a sinusoid, which will in turn make the weights more sinusoidal. This process will iteratively make the activity and the weights converge to the solution from Equation 17.

An important point in the use of Oja’s rule is that it would tend to concentrate the activity in a single harmonic. In a linear network, the harmonics would compete during learning, leading to one single harmonic emerging and all others being suppressed, as shown in section ‘Linear neurons’. For neurons with a nonlinear activation function, secondary harmonics would emerge, but would remain small under mild assumptions, as shown in section ‘Neurons with nonlinear activations’. Oja’s rule will still cause the weights to converge to approximately sinusoidal connectivity.

The finding that the weights will converge to a sinusoidal connectivity with learning has two interesting consequences from a biological standpoint: First, our circuit can emerge even when starting with only very coarse initial weights, without the need for high precision initial connectivity. Second, this simple plasticity rule allows the system to repair or recover from perturbations in its synapses as shown by simulations in Figure 5. However, this learning rule only applies to the weights that ensure stable, self-sustaining activity in the network. The network connectivity responsible for correctly integrating angular velocity inputs (given by the PEN to EPG connections in the fly) might require more elements than a purely Hebbian rule (Stringer et al., 2002), such as the addition of a predictive component (Vafidis et al., 2022).

Evolution of the eight-column circuit

Having derived and experimentally validated the theoretical circuit, we now address another question: whether there might be a reason that insect head direction circuits have an eight-column architecture (Honkanen et al., 2019; Stone et al., 2017; Pisokas et al., 2020). The derivations leading to Equation 17 are valid for other values of ${\displaystyle N}$, so there is no a priori reason to expect ${\displaystyle N=8}$.

Our reasoning follows recent studies in genetics (Johnston et al., 2022; Dingle et al., 2018), which argue that an observed organism is more likely to have resulted from a simple than a more complex genome: evolution favours simplicity. We note that powers of two are easier to generate with replication dynamics than other numbers, because they just require each cell to divide a set number of times. Other numbers require that, at some point, two cells resulting from a division must behave differently, necessitating more complex signalling mechanisms and rendering this possibility less likely to have been developed by evolution without some other driving factor. We therefore expect ${\displaystyle N}$ to be a power of two unless required to be otherwise.

As we show in section ‘Circuits with different neuron counts’, not all numbers of neurons enable a working circuit. The circuits for ${\displaystyle N=2}$ and ${\displaystyle N=4}$ are degenerate – either producing a single dimensional encoding, or two disconnected circuits that do not enforce the required circular topology. ${\displaystyle N=8}$ is the smallest power of two that could result in a non-degenerate circuit. This hints at the possibility that the eight-column architecture is not a chance evolutionary artefact, but rather that it is the genetically simplest circuit capable of performing heading integration.

Path integration dynamics for heading and position

Here, we make more precise statements about the dynamics of the circuit. The topology of the activity is a circular line attractor, so that any perturbation falls back into a circle and the position of the activity around the circle represents an angle. In our circuit with the dynamics from Equation 1, this circular line attractor is achieved by setting

where ${\displaystyle r}$ is the radius of the attractor. Given that the dynamics are effectively linear we can apply a Fourier transform and write the original dynamics around the line attractor directly in the spatial Fourier basis,

which we can expand to obtain

which only applies to the circular activity where the dynamics are linear. To incorporate the dynamics that force the activity to return to the cycle, we add a new term into the network dynamics which fulfils Equation 22 – increasing ${\displaystyle \Vert {\mathcal{F}}_f[\mathbf{a}]\Vert }$ if ${\displaystyle \Vert \mathbf{a}\Vert <r}$ and decreasing it if ${\displaystyle \Vert \mathbf{a}\Vert >r}$. This increase or decrease can be incorporated as a simple scaling factor on the activity decay term

where ${\displaystyle \alpha (r-\Vert \mathbf{a}\Vert )}$ is a nonlinear smooth function with a single minimum and ${\displaystyle \alpha (0)=1}$, which forces the activity to have the appropriate magnitude, ${\displaystyle \Vert \mathbf{a}\Vert =r}$.

Note that this doesn’t require the neurons to have access to a global activity magnitude signal because each neuron receives a sufficient number of inputs to locally compute the total activity in the network. See the ${\displaystyle f=1}$ circuit in Figure 3: each neuron receives inputs from all others except those tuned to an orthogonal direction, but the activity of these orthogonal neurons can be computed from the correlated neurons next to them. For example, neuron 1 receives input from all neurons except 3 and 7, but the activities of these neurons can be computed from the activities of neurons 2 and 4 or 6 and 8, respectively. Note also that this is not the case for degenerate circuits such as ${\displaystyle f=2}$ in Figure 3, which would require an activation function with access to global information to constrain the total activity.

As the activity always falls back to the circular line attractor, the heading integration is linear around the circle. This implies that any small movement of the animal or the perception of a sensory cue is first projected onto the circle, then linearly integrated. We therefore only consider how network inputs cause the activity to rotate around the network,

where ${\displaystyle {\mathcal{F}}_f[\mathbf{u}(t){]}_{\perp \mathbf{a}}}$ is the projection of the input signal perpendicular to the current activity pattern.

In more intuitive terms, the neurons have a saturating nonlinear activation function where they modulate their gain based on the total activity in the network. If the activity in the network is above the desired level, ${\displaystyle r}$, the gain is reduced and the activity decreases, and when the activity of the network is less than desired level, both the gain and the activity increase. Note that in this scenario transient deviations from the line attractor, which would induce nonlinear behaviour in the circuit dynamics, are tolerable. External inputs, ${\displaystyle \mathbf{u}(t)}$, could transiently modify the shape of the activity, producing activity shapes deviating from what the linear model can accommodate. For example, the shape of the bump attractor could be modified through nonlinearities while the insect attains high angular velocity (Turner-Evans et al., 2017).

Such nonlinear dynamics do not conflict with the theory developed here, which only requires linearity when the activity is projected onto the circular line attractor. In our framework, the linearity of integration at the circular line attractor is not a computational assumption, but rather it emerges from the principle of symmetry.

Ambiguities in multiple harmonics decoding with drift

We consider the case where multiple harmonics are used, and their phases have drifted from each other. We only focus on angular drift, rather than noise in the full activity, because as noted in section ‘Path integration dynamics for heading and position’, any deviations in the overall activity level in the network will dissipate.

For example, we consider the harmonics ${\displaystyle f_1}$ and ${\displaystyle f_2}$. The activity is given by Equation 5

If the phase of the second harmonic drifts by ${\displaystyle \delta \theta }$,

We can calculate the alignment between the activity with drift and the activity without drift as a dot product of the activity vectors,

However, there are other angles where the alignment is better. For example, we can consider an estimate angle, ${\displaystyle \hat{\theta }}$, with its corresponding activity, ${\displaystyle \mathbf{a}(\hat{\theta })}$, which gives

For small ${\displaystyle \delta \theta }$, this is maximised when ${\displaystyle \hat{\theta }=\theta +\delta \theta /2}$. However, since θ is circular, if the drift is ${\displaystyle \pi }$ then there are two possible positions given by ${\displaystyle \hat{\theta }=\theta \pm \pi /2}$.

This implies that an encoding using multiple harmonics does not necessarily offer a unique decoding in the presence of noise.

Equivalent circuits and degeneracies

Path counting

We counted the number of different paths between EPG neurons, accounting for the sign of the connections (whether the path passed through an inhibitory 7 neuron), and used the net path count as a proxy for connectivity strength. The results of this analysis are shown in section ‘Comparing the predicted circuit with biological data’, but we will detail the process here using the locust network from Pisokas et al., 2020, as an example, shown in Figure 6.

We consider the shortest excitatory and inhibitory pathways between EPG neurons, which have lengths of 2 and 3 respectively. There are no direct paths between ${\displaystyle \mathrm{EPG}}$ neurons in the locust circuit. In this case there are two paths of length 2 that implement self-excitatory connections for ${\displaystyle {\mathrm{EPG}}_1}$:

• ${\displaystyle {\mathrm{EPG}}_1\to {\mathrm{PEG}}_1\to {\mathrm{EPG}}_1}$

• ${\displaystyle {\mathrm{EPG}}_1\to {\mathrm{PEN}}_1\to {\mathrm{EPG}}_1}$

And there is one connection of path length 2 that connects ${\displaystyle {\mathrm{EPG}}_1}$ to each of its nearest neighbours:

• ${\displaystyle {\mathrm{EPG}}_1\to {\mathrm{PEN}}_1\to {\mathrm{EPG}}_2}$

• ${\displaystyle {\mathrm{EPG}}_1\to {\mathrm{PEN}}_1\to {\mathrm{EPG}}_8}$

For inhibitory connections there are four paths of length 3 that connect ${\displaystyle {\mathrm{EPG}}_1}$ to the neuron on the opposite side of the ring, ${\displaystyle {\mathrm{EPG}}_5}$. These are:

• ${\displaystyle {\mathrm{EPG}}_1\to {\mathrm{Delta7}}_5\to {\mathrm{PEN}}_5\to {\mathrm{EPG}}_5}$

• ${\displaystyle {\mathrm{EPG}}_1\to {\mathrm{Delta7}}_5\to {\mathrm{PEG}}_5\to {\mathrm{EPG}}_5}$

• ${\displaystyle {\mathrm{EPG}}_1\to {\mathrm{Delta7}}_4\to {\mathrm{PEN}}_4\to {\mathrm{EPG}}_5}$

• ${\displaystyle {\mathrm{EPG}}_1\to {\mathrm{Delta7}}_6\to {\mathrm{PEN}}_6\to {\mathrm{EPG}}_5}$

There are three connections of path length 3 connecting ${\displaystyle {\mathrm{EPG}}_1}$ to ${\displaystyle {\mathrm{EPG}}_4}$:

• ${\displaystyle {\mathrm{EPG}}_1\to {\mathrm{Delta7}}_4\to {\mathrm{PEN}}_4\to {\mathrm{EPG}}_4}$

• ${\displaystyle {\mathrm{EPG}}_1\to {\mathrm{Delta7}}_4\to {\mathrm{PEG}}_4\to {\mathrm{EPG}}_4}$

• ${\displaystyle {\mathrm{EPG}}_1\to {\mathrm{Delta7}}_5\to {\mathrm{PEN}}_5\to {\mathrm{EPG}}_4}$

Finally, there is one path of length 3 connecting ${\displaystyle {\mathrm{EPG}}_1}$ to the neuron perpendicular to it around the ring, ${\displaystyle {\mathrm{EPG}}_3}$:

• ${\displaystyle {\mathrm{EPG}}_1\to {\mathrm{Delta7}}_4\to {\mathrm{PEN}}_4\to {\mathrm{EPG}}_3}$

We then add these connections together to obtain the net path count profile for the network. Because the network is rotationally symmetric, the connections from ${\displaystyle {\mathrm{EPG}}_1}$ are the same as the connections from ${\displaystyle {\mathrm{EPG}}_2}$ but with the neurons indices incremented by 1. The generalised profile is shown in Table 3.

Data preprocessing

The following section details our data preprocessing to produce a central complex circuit for the fruit fly similar to that in Pisokas et al., 2020, using synaptic counts from the fruit fly connectome dataset (Scheffer et al., 2020). Full details are shown in our available code.

We first identified the 6 neuron types which corresponded to the neurons of interest in the central complex (Table 4).

We grouped the two PEN and two EPG populations together for further analysis.

From the connectome data we created a connectivity matrix (Figure 7) containing the number of synapses between all pairs of neurons of the above-mentioned types. This contained 129,473 synapses between 152 neurons in total.

Next, neurons in the same glomerulus were grouped together. The EPG neurons were divided between 18 glomeruli, L1-L9 and R1-R9. The PEN neurons were also divided between 16 glomeruli, L2-L9 and R2-R9. The PEG neurons were divided into 18 glomeruli with only 1 neuron in each, L1-L9 and R1-R9. The Delta7 neurons had 10 unique sub-types (Table 5) which were grouped into eight glomeruli based on their left glomerulus index – i.e., L4R5_R and L4R6_R were grouped in L4.

After this grouping we had a connectivity matrix of 60 neurons in total (Figure 8).

For the EPG, PEN, and PEG neurons we then grouped glomeruli mirrored in the two hemispheres together – i.e., L1 and R1 were grouped into glomerulus 1. This resulted in a connectivity matrix of 34 neurons (Figure 9).

Finally, as in Pisokas et al., 2020, we grouped the PEG9 and PEG1 neurons and EPG9 and EPG1 neurons together. The former because they both output to the same ellipsoid body segment, and the latter because they both receive common input. This resulted in a connectivity matrix with 32 neurons (Figure 10) which we used for the analysis in Figure 4 as this network had the eightfold rotational symmetry compatible with our theoretical model.

Fitting weights

Having calculated the paths or synaptic strengths between ${\displaystyle \mathrm{EPG}}$ neurons in the network, we next computed the average connectivity profile for the network – how each ${\displaystyle \mathrm{EPG}}$ neuron connected to its neighbours a certain distance around the ring. We then compared this connectivity profile to the closest fitting sinusoid to quantify to what degree our prediction of sinusoidal weights was consistent with biological data.

We use the vector

where ${\displaystyle d=((m-n)+4\mathrm{mod}8)-4\in [-\mathrm{4...3}]}$ is the signed circular distance from neuron ${\displaystyle m}$ to neuron ${\displaystyle n}$.

In the networks from Pisokas et al., 2020, the weights are rotationally symmetric, so we minimise the mean squared error between the sinusoidal profile and the biological connectivity profile,

by using least squares.

The network derived from the synaptic count data in Scheffer et al., 2020, is not rotationally symmetric so each value in the average connectivity profile, ${\displaystyle {\omega }_d}$, has a corresponding standard deviation, ${\displaystyle {\sigma }_d}$. In this case we minimise the precision-weighted mean squared error, which emphasises fitting connectivity profile values that are more consistently seen in the network:

The result of fitting sinusoidal profiles to the data is shown in section ‘Comparing the predicted circuit with biological data’. Since the fit for the fruit fly is not as clear as for the locust, we also compared the observed weights in the fruit fly with Gaussian and von Mises curves. Since both the Gaussian and von Mises distributions also have a parameter specifying their width, this requires fitting an extra parameter. For the Gaussian we have

where ${\displaystyle {\beta }_g,{\sigma }_g,{\gamma }_g}$ are parameters to fit. For the von Mises distribution we have

where ${\displaystyle {\beta }_v,{\kappa }_v,{\gamma }_v}$ are parameters to fit. For both the Gaussian and von Mises, we find a good agreement with the weights in Figure 11, and a lower root mean square error (RMSE) than with any combination of harmonics. However, we note that our sinusoidal model has two parameters instead of three, so to make a fair comparison we use the corrected Akaike information criterion (AICc) (Burnham and Anderson, 2004), which is given by

where ${\displaystyle p}$ is the number of distribution parameters, ${\displaystyle N}$ the number of samples, and ${\displaystyle \ln \left(L\right)}$ is the log-likelihood, which in this case is the sum of squared errors, ${\displaystyle \ln \left(L\right)=-N{\text{RMSE}}^2}$. We find that the lowest AICc (corresponding to the best model) corresponds to the harmonic ${\displaystyle f=1}$ (Table 6).

As a complementary approach to evaluate the shape of the distribution, we first fit the Gaussian and von Mises distributions to the best fit ${\displaystyle f=1}$ curve. We then freeze the width parameters of the distributions (${\displaystyle {\sigma }_g}$ for the Gaussian and ${\displaystyle {\kappa }_v}$ for the von Mises) and only optimise the amplitude and vertical offset parameters (β and γ) to fit the data. This approach limits the number of free parameters for the Gaussian and von Mises distributions to two, to match the sinusoid. The results are shown in Figure 11 and Table 6. Both the fixed-width Gaussian and von Mises distributions are a slightly better fit to the data than the sinusoid, but the differences between the three curves are very small.

In simplifying the fruit fly connectome data, we assumed all synapses of different types were of equal weight, as no data to the contrary were available. Different synapse types having different strengths could introduce nonlinear distortions between our net synaptic path count and the true synaptic strength, which could in turn make the data a better or worse fit for a sinusoidal compared to a Gaussian profile. As such, we don’t consider the 2% relative difference in RMSE between the ${\displaystyle f=1}$ sinusoid and fixed-width Gaussian and von Mises distributions to be conclusive.

Overall, we find that the cosine weights that emerge from our derivations are a very close match for the locust, but less precise for the fly, where other functions fit slightly better. Given the limitations in using the currently available data to provide an exact estimate of synaptic strength (for the locust), and due to the high variability of the synaptic count (for the fruit fly), we consider that our theory is compatible with the observed data.

Convergence of Oja’s rule with multiple harmonics