The HRC is the dominant model of motion computation in flies and other insects. In this paper we describe several generalizations of the HRC, but it is helpful to first review this canonical model. The HRC comprises three stages of processing. First, two different temporal filters (here, a low-pass filter $\left(f\left(t\right)\right)$ and a high-pass filter $\left(g\left(t\right)\right)$) are applied to each of two spatially filtered visual input streams (Figure 1A, ‘Materials and methods’). These four filtered signals are then paired and multiplied (Figure 1A). Finally, the HRC takes the difference between the two multiplied signals to obtain a mirror anti-symmetric motion estimator (Figure 1A). Because the HRC combines its two input channels via a multiplication operation, the average output of the HRC depends only on 2-point correlations in the visual stimulus. We thus refer to the HRC as a 2-point correlator, and we will return to this mathematical characterization of the HRC repeatedly throughout this work.

Figure 1

Download asset Open asset

No motion estimator is perfect for every stimulus, and this paper explores the hypothesis that evolution has tuned Drosophila's motion estimator for visual experiences that are likely to result from ordinary behavior in natural environments (Appendix 1). We assessed the accuracy of the HRC and other motion estimators by approximating naturalistic motion as the rigid translation of natural images (Clark et al., 2014), with a velocity distribution that mimicked Drosophila's natural behavior (Figure 1B, ‘Materials and methods’) (Katsov and Clandinin, 2008). We spatiotemporally filtered the input signals to simulate the responses of two neighboring photoreceptors (‘Materials and methods’). We quantified the performance of each model as the mean squared error between the input velocity and model output. However, we report each model's accuracy as the correlation coefficient between its output and the true velocity (Figure 1C), an intuitive metric that is equivalent to the mean squared error for correctly scaled model outputs (‘Materials and methods’). In isolation, the local HRC was weakly correlated with the velocity of motion (Figure 1D). Although the HRC's performance can be improved by averaging over space and time (Dror et al., 2001; Clark et al., 2014), this study explores how alternate nonlinear processing can improve motion estimation accuracy without sacrificing spatial or temporal resolution (Clark et al., 2014).

Researchers can probe a fly's motion estimate by measuring its behavioral optomotor turning response (Hassenstein and Reichardt, 1956; Götz and Wenking, 1973; Buchner, 1976; Reichardt and Poggio, 1976). We previously measured optomotor responses from flies walking on a spherical treadmill by recording their turning responses to various visual stimuli (Figure 1E) (Clark et al., 2014). We emphasized binary stimuli called gliders (Hu and Victor, 2010) (Figure 1F), which enforce spatiotemporal correlations to interrogate the fly's motion estimation algorithm. For example, 2-point gliders contain only 2-point correlations (first two stimuli, Figure 1F). Drosophila turned in response to these stimuli (Clark et al., 2014) (black bars, left, Figure 1G), and the HRC correctly predicted that flies would respond to both positive and negative 2-point correlations (gray bars, left, Figure 1G). On the other hand, 3-point gliders contain 3-point correlations without 2-point correlations (last four stimuli, Figure 1F). These stimuli generated motion responses in flies (black bars, right, Figure 1G) that the HRC could not explain (gray bars, right, Figure 1G). Thus, behavioral responses to glider stimuli show that the HRC is an incomplete description of fly motion estimation and provide a useful benchmark for evaluating alternate models.

In this study, we tune our models to optimize motion estimation accuracy, rather than to fit the behavioral data, for two main reasons. First, we want to explore the hypothesis that Drosophila's glider responses follow from performance optimization within biologically plausible circuit architectures. Second, we seek models that will generalize well across visual stimuli, and the measured glider responses under-constrain possible motion estimation models. It's useful to illustrate our procedure with a simple example. The HRC does not account for 3-point glider responses because it is insensitive to 3-point correlations. Nevertheless, 3-point correlations are present in natural stimuli (Clark et al., 2014; Nitzany and Victor, 2014), and their use might facilitate accurate motion estimation. We can explore this hypothesis by summing the HRC with a motion estimator designed to respond specifically to 3-point correlations. For instance, the mirror anti-symmetric ‘converging’ 3-point correlator multiplies one high-pass filtered signal with two low-pass filtered signals (Figure 1H) and mimics the converging structure present in certain glider stimuli (last two stimuli, Figure 1F). We tune the model for motion estimation accuracy by choosing the weights of the HRC and the converging 3-point correlator to minimize the mean squared error (‘Materials and methods’). The resulting model is more accurate than the HRC (Figure 1I) and it predicts that flies should respond to glider stimuli in the observed directions (Figure 1J, ‘Materials and methods’). Nevertheless, this simple model underestimates 3-point turning magnitudes (Figure 1J), indicating a discrepancy between the fly's motion estimator and this performance-optimized model.

In this study, we apply this same basic model building procedure to a series of increasingly general model architectures. There are four benefits to this approach. First, each model incorporates a type of computation that was neglected by earlier models. Thus, we can compare model accuracies to quantify how important various computations are for naturalistic motion estimation. Second, each model has a distinct biological interpretation in terms of Drosophila's motion estimation circuit (Table 1). This allows us to enumerate many directions for future experimental and computational research. Third, this set of models reveals several distinct principles of accurate naturalistic motion estimators, yet no single model illustrates every principle. Finally, by comparing the glider predictions of each model to behavioral data, we can gain insight into which principles underlie Drosophila's known glider responses.

The HRC correlates pairs of photoreceptor signals (Figure 1A). We previously assumed that each photoreceptor's response was generated from incoming contrast signals through linear spatiotemporal filtering. However, real photoreceptors are linear only over a limited range of inputs (Laughlin, 1981; Juusola and Hardie, 2001) (Table 1). Our first model thus modifies the HRC by allowing the photoreceptor responses to become nonlinear (Figure 2A). More specifically, we consider models in which a static nonlinearity transforms the filtered contrast signals before a standard HRC is applied to the two input streams (Figure 2A, ‘Materials and methods’). Since the nonlinearity occurs before the HRC, we refer to this model as the front-end nonlinearity model. By nonlinearly transforming the contrast signals, the front-end nonlinearity model is able to reshape natural sensory statistics. In particular, linear photoreceptor signals inherit complex non-Gaussian statistics from their natural inputs (Figure 2B), but front-end nonlinearities (Figure 2C) can produce transformed signals with alternate statistics (Figure 2D, ‘Materials and methods’). Thus, optimal front-end nonlinearity models should reshape natural statistics into those that best suit the HRC. Previous studies have already demonstrated example front-end nonlinearity models that improve naturalistic motion processing by the HRC (Dror et al., 2001; Brinkworth and O'Carroll, 2009). Here we provide new theoretical insight into these improvements and their consequences for glider responses.

Figure 2 with 1 supplement see all

Download asset Open asset

Although the statistics of natural images are complicated, the mean squared error between the HRC's output and the velocity of motion depends only on a few statistical quantities. Since the HRC is a 2-point correlator, the mean velocity signal decoded by an HRC is determined by the second-order statistics of the image ensemble (Dror et al., 2001). The variance of the motion signal comes from the square of a quadratic signal, and thus the noise statistics of the HRC depend on the fourth-order statistics of the image ensemble (Appendix 2). If the image ensemble is spatially uncorrelated, the situation simplifies further and the correlation between the estimated and true image velocity is determined entirely by the standardized fourth central moment of the input streams, a quantity known as kurtosis (Appendix 3). A larger kurtosis results in a larger error in the motion estimate. Note that some authors use ‘kurtosis’ to refer to the ‘excess kurtosis’, which shifts kurtosis values such that the Gaussian distribution has zero excess kurtosis. This shift is not relevant for our purposes. Because large positive contrasts are relatively probable, naturalistic inputs are highly kurtotic (kurtosis = 9.6 for the spatiotemporal filtering in our simulations) and are thus expected to hinder HRC performance (Figure 2B).

The Gaussian, uniform, and symmetric Bernoulli distributions have much lower kurtosis values (kurtosis = 3.0, 1.8, 1.0, respectively, Figure 2D). In fact, the symmetric Bernoulli distribution has the lowest kurtosis of any probability distribution (DeCarlo, 1997). When we transformed the HRC's inputs to have these statistics (‘Materials and methods’), we found that each nonlinearity substantially improved the accuracy of the HRC (Figure 2E). The contrast equalizing nonlinearity, which produces uniform outputs, performed best and also plays a prominent role in efficient coding theory (Laughlin, 1981). It is interesting that contrast equalization improved the accuracy of the HRC more than binarization (Figure 2E), even though it produced outputs with greater kurtosis. The reason for this is that natural images are spatially correlated, and the accuracy of the HRC over a general image ensemble depends on the ensemble's spatial correlation structure (Appendix 2). Binarization attenuated spatial correlations more strongly than contrast equalization over the natural image ensemble (Figure 2—figure supplement 1), and spatial correlations can enhance the performance of the HRC (Appendices 4, 5). Designing a nonlinearity that optimally sculpts the correlation structure of natural images is not simple and goes beyond the scope of this study.

Each front-end nonlinearity model is sensitive to a variety of higher-order correlations (Appendix 6). We thus tested whether accurate front-end nonlinearity models would predict Drosophila's glider response pattern. However, each front-end nonlinearity model performed poorly at this task (Figure 2F,G). None of the three models predicted that Drosophila would invert its response to positive and negative 3-point gliders (Figure 2G), even though they predicted that the 3-point glider responses would be nonzero. The simplest explanation for this observation is that the front-end nonlinearity models responded to fourth-order correlations that are common to the stimuli, rather than the third-order correlations that defined the glider stimuli and primarily drove the experimental response (Clark et al., 2014). Mechanistically, this result follows from the fact that the nonlinearities that reduced kurtosis (Figure 2C) were not strongly asymmetric around zero contrast (Appendix 6). The binarizing front-end nonlinearity model also failed to predict that Drosophila would respond less to negative 2-point glider stimuli than positive 2-point glider stimuli (Figure 2F). Since this effect was correctly predicted by the standard HRC (Figure 1G), this observation shows that accurate front-end nonlinearity models can distort the processing of 2-point correlations. Although the front-end nonlinearity model did not explain the phenomenon of fly glider perception, future work should investigate whether its merits make it functionally relevant for motion processing in other contexts or species.

Instead of a front-end nonlinearity, Drosophila could use an alternative non-Reichardtian motion estimation strategy that reflects natural sensory statistics, without necessarily requiring nonlinear preprocessing. Previous computational analyses show that motion estimation strategies that distinguish light and dark information can enhance motion processing with natural inputs (Fitzgerald et al., 2011; Clark et al., 2014; Nitzany and Victor, 2014), and recent experiments indicate that flies use separate channels to process the motion of light and dark edges (Joesch et al., 2010; Clark et al., 2011; Behnia et al., 2014; Clark et al., 2014; Meier et al., 2014; Strother et al., 2014) (Table 1). Our next model explores the hypothesis that Drosophila segregates ON and OFF signals in order to facilitate naturalistic motion estimation (Clark et al., 2014) (Figure 3A, ‘Materials and methods’). There are four ways to pair the ON and OFF components of the two filtered signals that enter the HRC's multiplier. For example, one possibility is to pair the ON component of the low-pass filtered signals with the OFF component of the high-pass filtered signal. Since each pairing restricts the HRC's multiplier to a single quadrant of the Cartesian plane, we refer to these four signals as HRC-quadrants. If the quadrants are summed with equal weights, then this model is mathematically identical to the HRC (Hassenstein and Reichardt, 1956; Clark et al., 2011). Unequal weighting coefficients enable the motion estimator to prioritize some quadrants over others, and here we select quadrant weightings that minimize the mean squared error between the model output and velocity (Figure 3B, ‘Materials and methods’). More generally, we refer to any model that linearly combines the four HRC-quadrants as a weighted 4-quadrant model. The precise manner in which the four HRC-quadrants might map onto circuitry remains unclear; we do not suggest there exists separate circuitry for each quadrant. For instance, studies have identified only two motion-processing channels in the Drosophila brain, which might suggest that the fly only uses a subset of the quadrants (Eichner et al., 2011; Joesch et al., 2013; Maisak et al., 2013). On the other hand, each channel appears imperfectly selective for light vs dark signals (Behnia et al., 2014), which in principle enables these two channels to access all four quadrants (Table 1).

Figure 3 with 2 supplements see all

Download asset Open asset

We began by examining how well individual quadrants predicted the velocity of motion. The four quadrants provided motion signals of strikingly different quality (first four red bars, Figure 3C). The most accurate quadrant correlated negative low-pass filtered signals with negative high-pass filtered signals ((− −) bar, Figure 3C). This isolated quadrant already outperformed the full HRC. The quadrant that correlated negative low-pass filtered signals with positive high-pass filtered signals also performed relatively well ((− +) bar, Figure 3C), whereas the quadrants that involved positive low-pass filtered signals performed poorly ((+ +) and (+ −) bars, Figure 3C). This shows that negative signals emanating from the low-pass filter better facilitate motion estimation, and the HRC's uniform weighting of all four quadrants is computationally detrimental.

We next considered all subsets of two, three, or four quadrants. The best subsets for each number of predictors were nested, and the quadrants were incorporated in the order (i) (− −); (ii) (− +); (iii) (+ +); (iv) (+ −). Although all four quadrants enhanced the accuracy of the weighted 4-quadrant model, the benefit of each added quadrant decreased with the number of quadrants (Figure 3C). It is possible to reparameterize the weighted 4-quadrant model in a form that isolates the contributions of various higher-order correlations to the model's accuracy (Appendix 7). Interestingly, this parameterization showed that nearly all the accuracy of the weighted 4-quadrant model can be obtained by supplementing the HRC with a set of odd-ordered correlations that account for the asymmetry between positive and negative low-pass filtered signals (Figure 3—figure supplement 1, Appendix 8). Principal component analysis (PCA) did not reveal this simple interpretation of the model's computation (Appendix 9).

The performance-optimized weighted 4-quadrant model also offered an interesting interpretation of Drosophila's glider response pattern. First note that the model preserved the HRC's response pattern to 2-point glider stimuli (compare left subpanels of Figure 3D and Figure 1G). More interestingly, the model predicted behavioral responses to 3-point glider stimuli that matched the experimentally observed turning directions, and even the response magnitudes were similar between the model and the data (right, Figure 3D). Nevertheless, the model's predictions were imperfect. The primary qualitative discrepancy was that the model failed to predict that positive 3-point glider stimuli would generate smaller turning responses than negative 3-point glider stimuli. The simplest interpretation for this experimental result is that flies might incorporate both 3-point correlations and 4-point correlations into their motion estimation strategy. In particular, since the positive and negative 3-point glider stimuli have inverted 3-point correlations and matched 4-point correlations, third-order and fourth-order correlations would have the same sign for one parity and opposite signs for the other parity. This observation makes it easier to understand the glider predictions of the weighted 4-quadrant model. The optimized model does a good job accounting for the direction and approximate magnitude of the glider responses because it draws heavily on second-order and odd-order correlations, but it fails to predict the 3-point glider magnitude asymmetry because it finds little added utility in higher-order even correlations (Figure 3—figure supplement 1, Appendix 8). This failure stems from architectural limitations in the weighted 4-quadrant model (Figure 3—figure supplement 2), so it is important to consider alternate model classes.

The previous section suggested that the segregation of light and dark signals by Drosophila's motion estimation circuitry might enhance naturalistic motion estimation in a manner that also generates the observed glider responses. In this section, we introduce three hierarchical models to investigate other features of Drosophila's circuit that might have functional consequences for the processing of natural stimuli and gliders (Table 1). We refrain from modifying the temporal filtering of the motion estimator, and we focus on its nonlinear architecture.

The first of these models recasts the HRC and the weighted 4-quadrant model in a more general architecture. This model is the class of mirror anti-symmetric models that apply a 2-dimensional nonlinearity to the low-pass filtered signal from one point in space and the high-pass filtered signal from a neighboring point in space (Figure 4A). Since the observed glider responses indicate that flies use higher-order correlations of both even and odd order, we model this 2-dimensional nonlinearity as a fourth-order polynomial (‘Materials and methods’). The HRC corresponds to the special case of this nonlinearity that multiplies the two inputs (left, Figure 4B). To emphasize how the model class in Figure 4A generalizes the HRC, we refer to it as the non-multiplicative nonlinearity model. In comparison, the weighted 4-quadrant model corresponds to a different nonlinearity that separately scales a pure multiplication in each quadrant of the Cartesian plane. Compared to the HRC, the optimized forms of both the weighted 4-quadrant model and the non-multiplicative nonlinearity model substantially attenuated positive low-pass filtered signals (middle and right, Figure 4B), though the non-multiplicative nonlinearity shows less attenuation. This model architecture provides enough flexibility to generate the glider response pattern (Figure 4—figure supplement 1).

Figure 4 with 2 supplements see all

Download asset Open asset

The non-multiplicative nonlinearity model relaxes some restrictions of the 4-quadrant model. This is prudent because the exact nonlinear transformations implemented by neural circuits in the Drosophila brain remain poorly understood. For example, T4 and T5 are the first direction-selective neurons in the fly brain (Maisak et al., 2013), but the mechanism by which they become direction-selective is not yet known. Furthermore, neurons upstream of T4 and T5 imperfectly segregate light and dark information (Behnia et al., 2014) and show overlap between the two motion pathways (Silies et al., 2013), suggesting that ON/OFF segregation may not be crisply realized. We will discuss this model's estimation accuracy and glider performance in the next section.

Drosophila's motion processing circuitry suggests two more generalizations of the non-multiplicative nonlinearity model. First, note that the non-multiplicative nonlinearity model inherits the HRC's assumption that each nonlinear unit only acts upon the low-pass filtered signal from one point in space and the high-pass filtered signal from the neighboring point (Figure 4A). In contrast, the converging 3-point correlator (Figure 1H) shows that the accuracy of motion estimation can sometimes be enhanced by nonlinearly combining both low-pass filtered signals (Figure 1I). Moreover, connectomic evidence conflicts with the non-multiplicative nonlinearity model's constraints, because each T4 cell receives synaptic connections from both the Mi1 cell and the Tm3 cells (T4's two major input channels) at overlapping points in space (Takemura et al., 2013). The unrestricted nonlinearity model removes this restriction of the non-multiplicative nonlinearity model by allowing a 4-dimensional nonlinearity to act on all four filtered signals (Figure 4C). Here, we again model this nonlinearity as a fourth-order polynomial (‘Materials and methods’). The unrestricted nonlinearity allows the motion estimator to nonlinearly combine multiple temporal channels from the same point in space. Recent experiments indicate that the Mi1 and Tm3 cells alone are insufficient to account for the motion processing of the T4 channel (Ammer et al., 2015). Future work might generalize the unrestricted nonlinearity model to include three or more temporal channels at each point in space.

The models presented so far operate only on a pair of neighboring photoreceptors, and the final generalization incorporates a third point in space. Averaging EMDs over space improves the accuracy of whole-field motion estimation (Dror et al., 2001), but Drosophila's neural circuitry suggests that it might adopt a more sophisticated strategy to combine signals across space. In particular, single T4 cells receive synaptic inputs from Mi1 cells and Tm3 cells from more than two retinotopic columns (Takemura et al., 2013). This arrangement could allow the circuit to incorporate higher-order correlations that are distributed across three or more spatial input channels. To explore whether this possibility has computational significance, we generalized the unrestricted nonlinearity model to provide unrestricted access to six temporal channels distributed across three points in space (Figure 4D). We refer to this model as the extra input nonlinearity model. We approximate its 6-dimensional nonlinearity as a fourth-order polynomial (‘Materials and methods’).

Having introduced the rationale behind the non-multiplicative, unrestricted, and extra input nonlinearity models, it is straightforward to examine their performance as motion estimators. First note that the polynomial non-multiplicative nonlinearity model was a better motion estimator (Figure 4E) than the weighted 4-quadrant model (Figure 3C). This implies that some useful signatures of naturalistic motion are not made accessible by simply segregating ON and OFF motion signals. Interestingly, this performance improvement is largely due to 3-point correlations, and models that exclude fourth-order polynomial terms still outperform the weighted 4-quadrant model (Figure 4—figure supplement 2). Third-order correlations are only useful for motion estimation because of light–dark asymmetries in natural stimulus statistics (Fitzgerald et al., 2011; Clark et al., 2014), so this result implies that ON/OFF segregation provides an imperfect way to account for the complexity of light–dark asymmetries found in the natural world. The non-multiplicative nonlinearity model also made novel use of low-order correlations to improve its motion estimate (Appendix 10).

The three models are hierarchical because the non-multiplicative nonlinearity model is a special case of the unrestricted nonlinearity model, which is itself a special case of the extra input nonlinearity model. Thus, we expect each model to perform at least as well as its predecessor, but it is possible that some circuit elaborations will not introduce useful computational cues. Nevertheless, we found that that the unrestricted nonlinearity model performed better than the non-multiplicative nonlinearity model, and the extra input model performed better than the average of two neighboring unrestricted nonlinearity models (Figure 4E). Therefore, both models incorporated novel computational signatures with relevance for visual motion estimation. Although the relative improvements were fairly small, it's worth noting that the improvement from spatial averaging is also small, and it is possible that the fly brain builds an accurate motion estimator by combining a large number of weak predictors of motion.

Each of these three generalized models predicted 2-point glider responses (Figure 4F) that closely resembled the standard HRC (left, Figure 1G). Each model also correctly predicted the experimental turning directions to each of the 3-point glider stimuli (Figure 4G). The magnitudes of the 3-point glider turning responses did not unambiguously favor any of the three hierarchical models (Figure 4G) or the weighted 4-quadrant model (right, Figure 3D). Each model did better on some stimuli and worse on others. Nevertheless, the predicted glider responses did make several interesting points. First, the extra input nonlinearity model predicted a clear asymmetry between positive and negative 3-point gliders (Figure 4G). This shows that some of the even-ordered correlations found in 3-point glider stimuli have relevance for naturalistic motion estimation. Second, the observation that each model provides qualitatively similar glider response patterns illustrates that animals could use multiple nonlinear mechanisms to access ethologically relevant higher-order correlations. Future experiments should directly assess the functional relevance of the different models in the hierarchy. Finally, the qualitative agreement between all of these predictions and the experimental data supports the general hypothesis that glider responses could reflect underlying nonlinear mechanisms that facilitate motion estimation in natural environments.

In this paper, we sequentially introduced several models in order to isolate specific ideas about the relationships between Drosophila's behavior, its motion estimation circuit, and the statistical demands of accurate motion estimation in natural environments. The front-end nonlinearity model explored an interesting candidate principle for visual motion estimation, but it conflicted sharply with fly behavior (Figure 2G) and excluded the conceptual insights offered by other models. For example, the front-end nonlinearities we considered eliminated the asymmetry between light and dark contrasts (Figure 2D), removing the need for separate ON and OFF processing. However, the remaining models embodied ideas that are complementary rather than exclusive, and these models should not be thought of as competitors. Instead, we will show here that the final, most general model incorporates the variety of conceptual points that were initially illustrated by specific models.

The structure of the non-multiplicative nonlinearity models can be directly plotted (Figure 4B), but is not easy to visualize the 6-dimensional nonlinearity that defines the extra input nonlinearity model. We therefore need an alternate technique to illustrate its computations. We proceed by leveraging three ideas. First, a wide variety of visual motion estimators can be expanded as an abstract series of multipoint correlators (e.g., see Poggio and Reichardt, 1980, Fitzgerald et al., 2011, and Appendices 6, 7, 11), and it is straightforward to pictorially represent a multipoint correlator (e.g., see Figure 1A,H, and more to come). In the extra input nonlinearity model, this expansion is immediate because we have already parameterized its 6-dimensional nonlinearity as a polynomial. Importantly, this expansion should be considered at the algorithmic level (Marr and Poggio, 1976), and we do not suggest that the wiring of brain circuits will reflect a large number of higher-order correlators. To the contrary, a large number of higher-order multipoint correlators may be implemented implicitly by high-dimensional nonlinearities suggested by Drosophila's visual circuitry. Second, we note that certain multipoint correlators can be recombined into a 2-dimensional non-multiplicative nonlinearity that facilitates easy comparisons with the HRC and weighted 4-quadrant models (e.g., see Figure 4B). Taken together, these two points mean that we can represent the performance-optimized extra input nonlinearity model in terms of non-multiplicative nonlinearity models and multipoint correlators, each of which are easy to represent graphically.

This graphical representation could be unwieldy because of the shear number of higher-order correlators in the model. Thus the third and final point is that we need a way to identify a relatively small number of terms that substantially improve the accuracy of motion estimation and illustrate the conceptual content of the model. To achieve this, we used lasso regression (Tibshirani, 1996) to identify models with fewer multipoint correlators that still enabled accurate motion estimation (‘Materials and methods’). This analysis revealed that fewer than half of all multipoint correlators were needed to account for the full accuracy of the extra input nonlinearity model (rightmost bars, Figure 5A). In fact, the accuracy of naturalistic motion estimation increased rapidly as the few correlators were sequentially added (left bars, Figure 5A), and a model that used 16 out of the 209 possible predictors was already able to produce 74% of the gain offered by the full extra input nonlinearity model (red bar, Figure 5A).

Figure 5 with 1 supplement see all

Download asset Open asset

The leading 16 predictors compactly illustrated how the extra input nonlinearity model recapitulates the conceptual advances offered by the other models (Figure 5B). Four of the predictors combine to implement a mirror-symmetric non-multiplicative nonlinearity model that acts on the first and second points in space (first term, Figure 5B). The dominant contribution to the nonlinearity is the HRC's multiplier, but an additional third-order term breaks the symmetry between positive and negative low-pass filtered signals (Figure 5—figure supplement 1). Thus, the extra input nonlinearity model approximately correlates neighboring points in space, as the HRC would suggest, but it differentially weights positive and negative low-pass filtered signals, like the weighted 4-quadrant model. It also replicates the main insight from the non-multiplicative nonlinearity model: the best treatment of asymmetric light and dark information need not be as simple as pure ON/OFF segregation. The model used another eight predictors to construct two more non-multiplicative nonlinearity models, one that surveyed the second and third points in space and another that surveyed the first and third points (first and second terms, Figure 5B, Figure 5—figure supplement 1). These components make the previously highlighted conceptual points and add the observation that spatial averaging improves estimates.

The final four predictors implemented two mirror anti-symmetric multipoint correlators (third and fourth terms, Figure 5B). In particular, two predictors went towards implementing a converging 3-point correlator that spanned the first and third spatial points (third term, Figure 5B). This estimator made the model's asymmetric treatment of light and dark signals more nuanced than permitted by the non-multiplicative nonlinearity model, and it also incorporated motion signals that combine multiple temporal signals from the same point in space. This latter point was the main conceptual motivation for the unrestricted nonlinearity model. Finally, the last two predictors implemented a 4-point correlator that combined temporal signals from three distinct points in space (fourth term, Figure 5B). This component reinforces the conceptual motivation for the extra input nonlinearity model and gives a concrete example of a computationally relevant higher-order correlator that is distributed across three points in space. It's interesting that the leading fourth-order correlator spanned three spatial points, because the extra input nonlinearity model was the first performance-optimized model that generated a substantially asymmetric response to positive and negative 3-point gliders (Figure 4G).

This paper set out with the goal of exploring whether the statistical demands of naturalistic motion estimation could provide a useful lens for interpreting features of Drosophila's behavior and neural circuitry that push beyond the canonical HRC. Although we have considered several interesting classes of visual motion estimators, the space of possible motion estimators is much larger (Figure 5C). For instance, these models have not explored the impact of temporal filter choice on naturalistic motion estimation. Nor have they assessed the possibility of more than two temporal filters, which is be suggested by anatomical (Takemura et al., 2013) and physiological (Ammer et al., 2015) experiments. More generally, the neural circuits contributing to Drosophila's motion estimator are still incompletely known, and the extent to which the fly brain's biological complexity reflects computational sophistication remains an open question. Theoretical considerations will be critical for resolving that question and pinpointing the most relevant principles underlying visual motion estimation.

We simulated the linear responses of neighboring photoreceptors to naturalistic motion using methods similar to previous work (Clark et al., 2014). We began with a database of natural images (van Hateren and van der Schaaf, 1998). We converted each natural image to a contrast scale, $C\left(\overrightarrow{x}\right)=\left(I\left(\overrightarrow{x}\right)-I_0\right)/I_0$, where $C\left(\overrightarrow{x}\right)$ is the contrast at the spatial point $\overrightarrow{x}$, $I\left(\overrightarrow{x}\right)$ is its intensity, and $I_0$ is the average intensity across the image. Since we only consider horizontal motion, we emulated the spatial blurring of Drosophila's photoreceptors in the vertical dimension by filtering across rows with a Gaussian kernel (FWHM = 5.7°). We then took the central row of each filtered image to represent a one-dimension variant of the natural image, denoted $c\left(x\right)$. We applied reflective boundary conditions to generate images that covered 360° and down-sampled each resulting image to 1° pixels by averaging. Photoreceptor blurring from signals in the horizontal dimension depends on the velocity of motion. In particular, we model the response of the ith photoreceptor as

where $T$ is a causal exponential kernel (timescale = 10 ms), $M$ is a Gaussian kernel (FWHM = 5.7°), $x_i$ is the location of the ith photoreceptor, and $\nu$ is the velocity of motion.

Each naturalistic motion comprised a randomly selected one-dimensional natural image, an offset to set the initial location of the photoreceptors, and a velocity drawn from a zero-mean normal distribution with a standard deviation of 90°/s. In this manner, we simulated the responses of three horizontally adjacent photoreceptors (spaced by 5.1°) to 5 × 10⁵ naturalistic motions (each with duration = 800 ms, time step = 5 ms). We then explicitly enforced left-right symmetry in the naturalistic ensemble by pairing each naturalistic motion with a new simulated motion, in which the natural image is reflected, the velocity is inverted, and the offset is chosen such that $\left\{V_1\left(t\right),V_2\left(t\right),V_3\left(t\right)\right\}$ in the new naturalistic motion is exactly $\left\{V_3\left(t\right),V_2\left(t\right),V_1\left(t\right)\right\}$ from its partner. The final symmetric ensemble thus consists of 10⁶ naturalistic motions.

The HRC applies two temporal filters to its photoreceptor inputs. We denote the kernels of the low-pass and high-pass filters as f and g, respectively, such that the output of a local HRC is

where * denotes convolution (Figure 1A). We consider the HRC's velocity estimate for a given naturalistic motion as its value at the final time point of the simulation. We model the filter kernels as

and

where $\tau$ = 20 ms and $g\left(t\right)$ is comparable to lamina monopolar cell responses (Clark et al., 2011; Behnia et al., 2014). We built the alternate motion estimators considered in this work from the same four filtered signals, $\left\{\left(f\ast V_1\right),\left(g\ast V_1\right),\left(f\ast V_2\right),\left(g\ast V_2\right)\right\}$, always considering the estimator's output at the final time point as its velocity estimate. Thus, none of our models modified the spatial or temporal processing of the HRC, reflecting our emphasis on how nonlinear processing might be tuned for naturalistic motion estimation. The global output of an array of HRCs would be obtained by pooling signals across space. Here we focus on spatiotemporally local strategies for motion estimation and at most pool motion signals across two neighboring motion detectors.

We evaluate motion estimators by the mean squared error between their output and the true velocity. To minimize the mean squared error of the HRC, we scale its output by $r^{\left(R\right)}{\sigma }_{\nu }/{\sigma }_R$, where $r^{\left(R\right)}$ is the correlation coefficient between the HRC's output and the velocity of motion, ${\sigma }_{\nu }$ is the standard deviation of the velocity distribution, and ${\sigma }_R$ is the standard deviation of the HRC's output. Once the HRC is scaled in this manner, its mean squared error is

More generally, this equation rewrites the mean squared error of any optimally scaled motion estimator in terms of its correlation coefficient with the velocity. All motion estimators considered in this paper are optimally scaled, and we find the correlation coefficient to be more intuitive than the mean squared error. We thus always report the performance of each motion estimator in terms of the correlation coefficient between the true and estimated velocity.

We fit the linear weighting parameters in the models of Figures 1I, 3A, 4A,C,D to maximize the estimation accuracy over a simulated ensemble of naturalistic motions. The formulas provided in subsequent sections of the ‘Materials and methods’ will cast each motion estimation scheme as a linear combination of a variety of motion predictors,

where the $x_i$ are nonlinear combinations of $\left\{\left(f\ast V_1\right),\left(g\ast V_1\right),\left(f\ast V_2\right),\left(g\ast V_2\right),\left(f\ast V_3\right),\left(g\ast V_3\right)\right\}$ that depend on the model architecture, and the w_i are associated weighting coefficients. We chose the weights to minimize the mean squared error between the true and predicted velocity, which is the standard scenario considered by ordinary least-squares regression. The same weights maximize the correlation coefficient between the true and predicted velocity, and we typically present model accuracies as correlation coefficients.

We used twofold cross-validation to protect against over-fitting. In particular, we randomly divided the ensemble of naturalistic motions into a training set of 500,000 symmetrically paired examples and a testing set of the remaining 500,000 examples. We determined the weighting coefficients by minimizing the empirical error over the training set, and we reported accuracies over the test set. To estimate error bars for each model's accuracy, we computed twenty random divisions of the naturalistic motion ensemble and calculated the standard deviation of the estimation accuracy.

We generated 25 random instantiations of each glider stimulus considered by our previous experimental work (Figure 1F, duration = 3 s, update rate = 40 Hz, pixel size = 5°) (Hu and Victor, 2010; Clark et al., 2014). We evaluated the response of each model to these stimuli by averaging the outputs of 60 identical local motion estimators (each separated by 5.1°) over the last two seconds of visual stimulation. Glider predictions were equal and opposite for the left and right variants of the stimuli, so we pooled leftward and rightward stimuli in all figures (Figure 1F shows the rightward variants). We scaled each model's output such that the average response to the positive 2-point glider was 1. All figures associated with glider responses show the mean and standard error of each model's response across the 25 glider instantiations.

The model in Figure 2A replaces the linear photoreceptor signals, $V_1$ and $V_2$, with nonlinear photoreceptor signals

where $h$ is some nonlinear function. Thus, the motion estimate from the front-end nonlinearity model is

To implement the contrast equalizing nonlinearity, we replaced values of $V_i\left(t\right)$ by their rank-order (scaled and shifted to range between −1 and +1). Note that all $V_i\left(t\right)$ were sorted together (i.e., including all spatial points, temporal points, and simulated naturalistic motions). When multiple $V_i\left(t\right)$ had the same value, they were given the same rank. To implement binarizing nonlinearities, we again sorted the $V_i\left(t\right)$ and found the values corresponding to the threshold locations. For example, to calculate the binarizing nonlinearity with two steps (Appendix 4): (i) we found the $V_i\left(t\right)$ values corresponding to the 25th and 75th percentiles; (ii) signals below the 25th percentile or above the 75th percentiles were assigned the value of −1; and (iii) signals between 25th and 75th percentiles were assigned the value of +1. To implement the Gaussianizing nonlinearity, we again rank-ordered the $V_i\left(t\right)$ (scaled to range between 0 and 1) and applied the inverse Gaussian cumulative distribution function to these ranks. The HRC is the special case of this model where the front-end nonlinearity is linear.

The weighted 4-quadrant model in Figure 3A separately correlates bright and dark signals. Mathematically, it is

where $w_{ab}^{\left(Q\right)}$ are adjustable weights that parameterize the model,

${\left[x\right]}_{+}$ is equal to $x$ when $x$ is positive and zero otherwise, and ${\left[x\right]}_{-}$ is equal to $x$ when $x$ is negative and zero otherwise. The HRC is the special case of this model where $w_{++}^{\left(Q\right)}=w_{+-}^{\left(Q\right)}=w_{-+}^{\left(Q\right)}=w_{--}^{\left(Q\right)}$.

The non-multiplicative nonlinearity model in Figure 4A replaces the HRC's multiplication step with a more flexible two-dimension nonlinearity. In particular, it is

where we approximate the nonlinearity, $\eta$, as a fourth-order polynomial

and $w_{ij}^{\left(N\right)}$ are adjustable weights that parameterize the model. We include terms up to fourth order in this model to ensure that it is flexible enough to describe the published glider response data. In particular: (i) the second-order terms accommodate responses to 2-point glider stimuli; (ii) the third-order terms accommodate parity-inverting responses to 3-point glider stimuli; and (iii) the fourth-order terms enable the model to respond with unequal magnitude to positive and negative parity 3-point glider stimuli (Figure 4—figure supplement 1). Thus, this model has 14 parameters. The HRC is the special case of this model where only $w_{11}^{\left(N\right)}$ is nonzero.

Here we model the 4-dimensional nonlinearity in Figure 4C as a fourth-order polynomial of the four filtered signals in the HRC. In general, this motion estimator is

where $w_{ijkl}^{\left(S\right)}$ are adjustable weights that parameterize the model, and we set $w_{0000}^{\left(S\right)}=0$ because this term has no utility for naturalistic motion estimation. Thus, this model has 69 parameters. The HRC is the special case of this model where $w_{1001}^{\left(S\right)}= -w_{0110}^{\left(S\right)}\ne 0$, and all other parameters are zero.

Here we model the 6-dimensional nonlinearity in Figure 4D as a fourth-order polynomial of the six filtered signals in two neighboring HRCs. In general, this motion estimator is

where $w_{ijklmn}^{\left(E\right)}$ are adjustable weights that parameterize the model, and we set $w_{000000}^{\left(E\right)}=0$ because this term has no utility for naturalistic motion estimation. Thus, this model has 209 parameters. The average of two neighboring HRCs is the special case of this model where $w_{100100}^{\left(E\right)}= -w_{011000}^{\left(E\right)}=w_{001001}^{\left(E\right)}=-w_{000110}^{\left(E\right)}\ne 0$, and all other parameters are zero.

Lasso regression augments the squared error with an L₁ penalty on nonzero weighting coefficients that favors sparse solutions (Tibshirani, 1996). We used lasso regression to identify subsets of predictors that might enable accurate motion estimation (Figure 5A). Once we identified a predictor subset using lasso regression, we refit the nonzero model weights using ordinary least squares regression (i.e., without the weight penalty).

The pattern of light that stimulates the retina encodes information about the relative motion between the retina and its visual environment. The manner in which this information is encoded depends on the geometry of the photoreceptor array, the statistics of self-motion, and the statistics of the visual environment. The principal goal of this paper is to illustrate several ways that the brain's nonlinear processing of visual motion signals might be tuned to reflect specific features of the natural visual environment. We thus begin by enumerating some computational signatures of visual motion in natural environments, thereby exposing a diversity of stimulus features that visual system nonlinearities might aim to extract.

In the real world, animals encounter visual environments that are intricately structured and far from random (Appendix figure 1A) (Ruderman and Bialek, 1994; van Hateren and van der Schaaf, 1998; Geisler, 2008). When an animal rotates with constant angular velocity through the environment, the spatiotemporal response profile of the photoreceptor array encodes the velocity of self-motion through the slope of oriented streaks in space-time (front face, Appendix figure 1B) (Adelson and Bergen, 1985). Thus, a visual system with a dense array of noiseless photoreceptors could extract the angular velocity of an arbitrary image by computing the ratio of temporal and spatial derivatives (Potters and Bialek, 1994). The statistics of the image ensemble become relevant once multiple interpretations of the sensory world become logically consistent with the photoreceptor data. In particular, the optimal motion estimator depends on the statistics of the image ensemble when photoreceptors have noise (Potters and Bialek, 1994; Fitzgerald et al., 2011), and a nonzero spacing between photoreceptors introduces ambiguity via aliasing (Potters and Bialek, 1994). In these cases, the animal can use prior information regarding the sensory environment and its motion in order to weigh the plausibility of each sensory interpretation.

Appendix figure 1

Download asset Open asset

Full field motion transforms spatial features (top face, Appendix figure 1B) into temporal features (side face, Appendix figure 1B) in a manner that depends upon the velocity of motion. Consequently, one can also think about the visual signatures of motion in terms of spatiotemporal correlations between photoreceptors. The luminance contrast encoded by the $i$th photoreceptor is $C_i\left(t\right)=\left(I_i\left(t\right)-I_0\right)/I_0$, where $I_i\left(t\right)$ is the luminance intensity seen by the $i$th photoreceptor at time $t$ and $I_0$ is the average luminance intensity over the visual field. Thus, the average contrast is zero, and the simplest correlation function corresponds to the product of two spatially separated contrast signals. Measured over an ensemble of natural images, this 2-point correlation function had a global maximum at zero spatial offset (Appendix figure 1C). Consequently, the velocity of motion is encoded by the peak of the temporal cross-correlation function between two neighboring photoreceptors, which occurs at the temporal offset that equals the photoreceptor spacing (5.1° for Drosophila) divided by the velocity of motion (Appendix figure 1D). Natural images also contain many higher-order correlations (Ruderman and Bialek, 1994; Geisler, 2008). For instance, the nonzero skewness of natural images implies that the third-order correlation that multiplies the contrast at one point with the squared contrast at a neighboring point also has a peak at zero spatial offset (Appendix figure 1E). Correspondingly, the peak of the temporal 3-point correlation function between neighboring photoreceptors encodes the velocity of motion (Appendix figure 1F). This argument generalizes to nth-order correlation functions when the ensemble of natural images has a nonzero nth moment. Note that this argument does not necessarily imply that a motion estimator would benefit from the incorporation of all nonzero correlation functions, because the velocity signals provided by one correlation function could be redundant with those provided by others.

Importantly, photoreceptor correlation functions also encode velocity information away from their peaks. For example, the velocity of motion influences the widths of the temporal cross-correlation functions between pairs of photoreceptors (Appendix figure 1D,F). To see this, note that the values of the temporal correlation functions at zero temporal offsets are velocity independent, whereas the peak locations are closer to zero for larger speeds (Appendix figure 1D,F). This implies a more rapid falloff for higher speeds. This fundamental effect occurs because nearby points are more correlated in natural environments and photoreceptors rapidly survey distant points when the speed of motion is high.

The description above illustrates how visual motion becomes encoded in photoreceptor correlations. A central goal of research in visual motion estimation is to understand how neural circuits invert (or decode) that encoding of velocity. Just as a broad class of functions can be represented as a power series, a broad class of motion estimators can be represented as a Volterra series (Poggio and Reichardt, 1980; Fitzgerald et al., 2011). Each term in the Volterra series can be interpreted as a multipoint correlator that decodes velocity information from a specific correlation function (Fitzgerald et al., 2011). For example, the HRC and the motion energy model are 2-point correlators that decode velocity from 2-point correlations, whereas the Bayes optimal motion estimator capitalizes on a wider variety of correlation functions (Potters and Bialek, 1994; Fitzgerald et al., 2011). Because multipoint correlators relate intuitively to measurable properties of the image ensemble, we will find that decomposing a motion estimator in terms of multipoint correlators is often illuminating. Moreover, we will use multipoint correlators as a common basis to compare the computations performed by mechanistically distinct models.

In this section we derive an expression for the accuracy of a general 2-point correlator in terms of the statistics of naturalistic motion.

We consider a general 2-point correlator that temporally correlates visual signals from the spatial points i and j. Mathematically, this estimator has the form

where the 2-point kernel, $k_{i,j}^{\left(2\right)}\left(t_1,t_2\right)$, defines the correlator by specifying how each 2-point correlation contributes to the motion estimate. We model the response of the ith photoreceptor as

where $T$ is a temporal integration kernel, $M$ is the photoreceptor's spatial acceptance profile, ${\theta }_i$ is the location of the ith photoreceptor, $c\left(\theta \right)$ is the spatial contrast pattern of the visual world, and $\nu \left(t\right)$ is the time-dependent velocity. This formula simplifies to the formula in the ‘Materials and methods’ when $\nu \left(t\right)$ is time-independent. If $T$ is an invertible linear filter, then a more convenient representation of the photoreceptor signals is

where $U_i=T^{-1}\ast V_i$, $\mathcal{C}=M*c$, and $*$ is the convolution operator (Potters and Bialek, 1994). We can rewrite the 2-point correlator in this representation as

where

is the 2-point kernel that converts correlations in the $U$ variables to a velocity estimate.

Recall that we quantify the performance of visual motion estimators based on the mean squared error between the true and estimated velocities

where ${\sigma }_{\nu }$ = 90°/s is the standard deviation of the velocity distribution. For estimators that are scaled to minimize their mean squared error (‘Materials and methods’), this formula can be rewritten as

where

is the correlation coefficient between the estimated and true velocities. Thus, minimizing the mean squared error is mathematically equivalent to maximizing the correlation coefficient if all motion estimators are correctly scaled. We find the correlation coefficient to be a more intuitive error metric than the mean squared error, so many of our results will be presented in terms of correlation coefficients.

The numerator of the correlation coefficient is determined by the second-order statistics of the image ensemble,

where Δ_ij is the angular separation between the ith and jth photoreceptors, and

is the 2-point correlation function over the ensemble of spatially filtered natural scenes. Note that ${\mathcal{C}}^{\left(2\right)}$ is independent of $x$ because reasonable image ensembles are translationally invariant. Also note that the 2-point correlation function of filtered natural images is related to the correlation function of unfiltered images by

where $C^{\left(2\right)}\left(\mathrm{\Delta }\right)$ is the correlation function of unfiltered images, and we've assumed that $M$ is a symmetric function. We model $M$ as Gaussian with FWHM of 5.7°, so $M*M$ is also Gaussian with FWHM of $\sqrt{2}\times {5.7}^{\circ }={8.1}^{\circ }$.

On the other hand, the denominator of the correlation coefficient is determined by fourth-order statistics of the image ensemble,

where

is the 4-point correlation function of the ensemble of filtered natural images. Notice that the second argument of ${\mathcal{C}}^{\left(4\right)}$ in Equation 12 lacks the additive factor of ${\mathrm{\Delta }}_{\mathit{ij}}$ because $U_i\left(t-t_1\right)$ and $U_i\left(t-t_3\right)$ correspond to the same point in space. As above, ${\mathcal{C}}^{\left(4\right)}$ is related to the unfiltered 4-point correlation function through a fourfold application of the photoreceptor spatial acceptance filter.

The preceding analysis shows that only the second-order and fourth-order statistics of the natural image ensemble contribute to the correlation coefficient between an arbitrary 2-point correlator and the true velocity. The same quantities also determine the mean squared error. Thus, the second-order and fourth-order statistics of the image ensembles are the critical determinants of a 2-point correlator's motion estimation accuracy. Note that both the HRC and the motion energy model fall into this important class of visual motion estimators, so our analysis is also important for understanding visual motion estimation by vertebrates.

In this section, we apply the results of Appendix 2 to the special case of normally distributed velocities and spatially uncorrelated image ensembles. This calculation reveals an important role for kurtosis in motion estimation, and we discuss how nonlinearities in the early visual system could cope with highly kurtotic naturalistic inputs.

In this section, we assume that the velocity is time-independent (i.e., $\nu \left(t\right)=\nu )$ and normally distributed

We also assume that the image ensemble is spatially uncorrelated. By this, we mean that the luminance contrast at each point in space is statistically independent of the luminance contrast at all other points in space. Thus, the second-order correlation function is

where ${\sigma }_C$ is the standard deviation of the luminance contrast, and $\delta \left(\mathrm{\Delta }\right)$ is the Dirac delta-function. The fourth-order correlation function is

where $k_4$ is the excess kurtosis of the contrast distribution. The excess kurtosis is zero for normally distributed contrasts. It can either be positive or negative for other contrast distributions. Note that we define the kurtosis of a probability distribution to be its fourth central moment normalized by the square of its second central moment. Thus, the kurtosis of a normal distribution is 3. We caution readers that some other sources use ‘kurtosis’ to refer to the excess kurtosis.

With these assumptions, the signal term represented by Equation 9 is

and the noise term represented by Equation 12 is

where we've assumed that the 2-point correlator is mirror anti-symmetric,

in order to ignore contributions from static signals. This mirror-symmetry assumption holds for the HRC and the motion energy model. Since the denominator of the correlation coefficient is set by $\sqrt{\langle {\left(v_e^{\left(2\right)}\left(t\right)\right)}^2\rangle }$, both the signal and the noise are proportional to ${\sigma }_C^2$. Thus, the only remaining dependence on the image ensemble is through the excess kurtosis. Note that

Thus, the correlation coefficient is maximized by making $k_4$ as small as possible.

In conclusion, if the image ensemble is spatially uncorrelated (at second and fourth-order), then the image ensemble only affects the correlation coefficient between the velocity and a 2-point correlator through its kurtosis. The best accuracy is achieved when the kurtosis is minimized. In reality, the assumption that the image ensemble is spatially uncorrelated is clearly wrong. Natural images are strongly correlated, and even if they weren't, they'd become correlated once they are filtered by the photoreceptors' spatial acceptance filter. Nevertheless, Figure 2E empirically shows that introducing several front-end nonlinearities that decrease the kurtosis also improve the accuracy of naturalistic motion estimation. Thus, kurtosis provides a useful guide for the design of neuronal nonlinearities. On the other hand, Figure 2D,E demonstrate that it's too simplistic to assume that the kurtosis is the only relevant factor for the accuracy of a 2-point correlator. As we'll discuss in the next section, spatial correlations in the image ensemble also affect the accuracy of 2-point correlators.

When we applied a contrast-equalizing or binarizing nonlinearity to naturalistic inputs before evaluating the HRC, we found that both nonlinearities substantially improved the accuracy of the HRC (Figure 2E). Interestingly, contrast equalization improved the accuracy of the HRC more than binarization (Figure 2E), even though it produced outputs with greater kurtosis. The reason for this is that natural images are correlated (Appendix figure 1), and the accuracy of the HRC over a general image ensemble depends on the ensemble's spatial correlation structure (Appendix 2). Binarization attenuated spatial correlations more strongly than contrast equalization over the natural image ensemble (Figure 2—figure supplement 1), which leads us to hypothesize that correlations present in the natural image ensemble might benefit the HRC's performance. In Appendix 5 we will provide theoretical support for this idea. Here we begin with a less mathematical argument that also supports our hypothesis.

A comparison between the estimation performance of binarizing and equalizing front-end nonlinearities was complicated by the fact that the models produced outputs that differed in both their point statistics and their correlation structures. To gain more direct insight into how spatial correlations affect motion estimation performance, it would be helpful to compare front-end nonlinearity models that differ only through their output correlation structures. We implemented this comparison using a family of binarizing front-end nonlinearities that undergo multiple steps between +1 and −1 (Appendix figure 2A). Although these nonlinearities are not physiologically realistic, they are conceptually useful because they each produced a stimulus ensemble that minimized the kurtosis yet achieved distinct correlation structures (Appendix figure 2B). These nonlinearities thus allow us to assess directly whether spatial decorrelation of inputs degrades the motion estimation performance of the HRC. We found that each binarizing front-end nonlinearity model outperformed the original HRC (Appendix figure 2C). However, we found that the magnitude of the improvement decreased with the number of steps (Appendix figure 2C). Since spatial cross-correlations also decreased as a function of the number of steps (Appendix figure 2B), these results support our hypothesis that the correlations present in natural visual inputs aid the functionality of the standard HRC.

Appendix figure 2

Download asset Open asset

The HRC correlates two signals that are offset in space and differentially delayed in time. One intuition that researchers often apply to this computation is that the correlation operation effectively detects times when two signals that are offset in space and time are equal. However, a motion estimator that strictly obeyed this intuition would be agnostic to the spatial correlation structures present in the input signals, and our results show that the HRC is not (see also Appendix 5). Instead, the HRC also generates motion signals when its two input channels are imperfectly aligned, and these signals depend strongly on the correlation structure of the inputs (Appendix figure 1D). Our results thus show that the HRC's ability to detect imperfect coincidences contributes significantly to its performance as a motion estimator, as was suggested intuitively in Appendix 1.

In this section, we apply the results of Appendix 2 to the special case of normally distributed velocities and normally distributed image ensembles. This model formalizes how spatial correlations in the natural world affect the accuracy of motion estimation by 2-point correlators and shows how spatial decorrelation can adversely affect estimation accuracy. For example, we'll show that the simplest HRC is unable to extract motion signals from high frequency components of the image ensemble, yet those components still lead to variability in the motion estimator. Thus, this HRC works best when the image ensemble is correlated in a manner that avoids high-frequency components in the signal, and spatial low-pass filtering at the photoreceptor level can help to eliminate the high-frequency image components that hurt the HRC's accuracy.

Here we use the same velocity distribution that we used in Appendix 3 (i.e., Equation 14). However, we now allow the two point correlation function to have arbitrary structure

where $S_k$ are the Fourier coefficients for ${\mathcal{C}}^{\left(2\right)}\left(\mathrm{\Delta }\right)$, and we have noted that the image ensemble is $2\pi$-periodic. Note that $S_k$ is called the power spectrum of the image ensemble, and uncorrelated ensembles correspond to the special case where $S_k$ = constant. With these assumptions

By evaluating the integral, we find that this velocity expectation is

Thus, if we define

then

Each frequency component of the image ensemble linearly contributes to the correlation between the 2-point correlator's response and the velocity. The weight of each frequency component is determined by the structure of the 2-point correlator and the width of the velocity distribution.

We compute the fourth-order moment of the image ensemble using Wick's theorem for Gaussian moments, which says

This immediately implies that

Once again, it's convenient to rewrite this expression in the Fourier domain

With these assumptions

We evaluate the expectations over velocity by noting that each has the form

for some spatial offset $\mathrm{\Delta }$ and temporal offsets $\left\{{\delta }_1,{\delta }_2\right\}$. Thus, if we define

then

Power spectrum components contribute to the 2-point correlator's variance quadratically.

Putting these pieces together, the expected squared error achieved by a 2-point correlator is a quadratic function of the power spectrum

We're interested to know whether spatial correlations can enhance the accuracy of the 2-point correlator. This will be the case unless a uniform power spectrum minimizes $ϵ$. Note that every physically meaningful power spectrum is non-negative

Thus, the minimum of $ϵ$ either occurs at an extremum point or on the boundary of admissible solutions. If the minimum occurs on the boundary, then a subset of the $S_k$ are exactly equal to zero. In particular, the power spectrum would not be constant, which implies that the image ensemble would be spatially correlated. At an extremum point, we must find

for every $k$. A uniform power spectrum can only satisfy this condition if

where $\beta >0$ is the (constant) value of each power spectrum component. This is generally not the case, so correlations exist that would help typical 2-point correlators.

For example, the simplest HRC, which replaces the low-pass and high-pass filters with pure time delays, is

where $A$ is a constant with units of °/s. For this model,

Substituting this expression into the above formulas, we find

and

where ${\mathrm{\Delta }}_0$ is the spacing between adjacent photoreceptors. Note that

On the other hand,

This does not approach zero, even for large values of $k_1$. Therefore, ${{{\displaystyle \sum }}^{\text{}}}_{k\prime =0}^{\infty }{\mathrm{\Gamma }}_{kk\prime }$ diverges and ${\gamma }_k\ne \beta {{{\displaystyle \sum }}^{\text{}}}_{k\prime =0}^{\infty }{\mathrm{\Gamma }}_{kk\prime }$. In this model, high frequency components lack signal but contribute noise. It's helpful if these frequency components are absent from the image ensemble. Future work should more fully investigate the role of spatial correlations in naturalistic motion estimation.

The response of the front-end nonlinearity model to a 3-point glider stimulus is determined by the higher-order correlations that it detects in the stimulus. Furthermore, we argued in Appendix 1 and Figure 1I that higher-order correlations can contribute to the accuracy of visual motion estimators. We now describe how front-end nonlinearities provide pair-correlation mechanisms with access to certain types of higher-order correlations.

We suppose that the front-end nonlinearity, denoted $h$, has a power series expansion:

Then the cross-correlation function between two non-linearly transformed input streams, denoted $y_1$ and $y_2$, is

where $V_1$ and $V_2$ are linear photoreceptor signals. This substitution explicitly demonstrates that the front-end nonlinear transformation enables pair correlation mechanisms to incorporate higher-order correlations of the form $\langle V_1^n\left(t\right)V_2^m\left(t+\tau \right)\rangle$. The choice of nonlinearity specifies the expansion coefficients, $h_n$, which in turn determines the pattern of higher-order correlations that the pair correlator incorporates into its velocity estimate. For example, sensitivity to odd-ordered correlations demands that $h_n$ be large for some even values of $n$. These expansion coefficients would manifest themselves in the structure of the front-end nonlinearity as asymmetries between positive and negative contrasts, but strong asymmetries were not needed to eliminate kurtosis in natural image ensembles (Figure 2C). Inversely, one could use this equation to determine whether a set of expansion coefficients exist that would implement a desired series of multipoint correlators. The preceding argument implies that strongly asymmetric front-end nonlinearities would be needed to account for the 3-point glider responses.

In this Appendix, we rewrite the weighted 4-quadrant model in a basis that isolates its dependence on 2-point correlations, on higher-even-ordered correlations, and on two types of odd-ordered correlations. In Appendix 8, we'll discuss the motion estimation performance of the weighted 4-quadrant model in this basis in order to gain insight into why performance-optimized weighted 4-quadrant models also predict 3-point glider responses that resemble Drosophila behavior.

The weighted 4-quadrant model supposes that the input signals are segregated into four separate streams: