Fly Photoreceptors Encode Phase Congruency

More than five decades ago it was postulated that sensory neurons detect and selectively enhance behaviourally relevant features of natural signals. Although we now know that sensory neurons are tuned to efficiently encode natural stimuli, until now it was not clear what statistical features of the stimuli they encode and how. Here we reverse-engineer the neural code of Drosophila photoreceptors and show for the first time that photoreceptors exploit nonlinear dynamics to selectively enhance and encode phase-related features of temporal stimuli, such as local phase congruency, which are invariant to changes in illumination and contrast. We demonstrate that to mitigate for the inherent sensitivity to noise of the local phase congruency measure, the nonlinear coding mechanisms of the fly photoreceptors are tuned to suppress random phase signals, which explains why photoreceptor responses to naturalistic stimuli are significantly different from their responses to white noise stimuli.


Introduction
In his seminal papers [1][2] Horace Barlow postulated that sensory pathways are tuned to detect, enhance, and efficiently encode the stimuli that are important for survival. However, as Barlow pointed out [1], the signal transformations performed by sensory neurons are difficult to characterize using ordinary physiological investigations. In particular, the responses of sensory neurons contain components from linear and nonlinear transductions that are difficult to separate [3].
Although it has been known for some time that sensory neurons, whether auditory [4][5] [6], olfactory [7] or visual [8][9] [10], respond nonlinearly when driven by stimuli that have 'naturalistic' properties, the nonlinear relationship between the statistical properties of the stimuli and the neuron responses-or in other words, the computations performed by these neurons [11]-have not been fully characterized.
In the visual system, the detection of the boundary or edges of objects is crucial for object segregation, categorization and recognition as well as for motion detection [12]. In this context, the temporal structure of the retinal images is very important [12] since moving spatial edges generated by self-motion or by moving objects produce temporal edges at the photoreceptor level. Therefore, selectively enhancing the salience of these temporal features should facilitate downstream processing of the spatio-temporal visual stimuli for motion detection. As edges correspond to points of local maximum phase alignment of the constituent Fourier components, local phase congruency, which is invariant to illumination and contrast, is an accurate measure of edge saliency that can be encoded by the photoreceptors [13] [14] [15].
There are many clues, which suggest photoreceptors are tuned to distinguish and selectively process the temporal phase correlations present in light stimuli. These phase correlations are biologically relevant [4][16] [17].
It is well known, for example, that photoreceptor responses to naturalistic stimuli are highly nonlinear [8] whereas Gaussian white noise stimuli tend to linearize the response [18].
Naturalistic stimuli exhibit local and global phase correlations caused by the edges, contours and textures present in the natural scene [12,19], which can be described by higher-order statistics [20][21][22][23] and can only be encoded by applying nonlinear transformations to the stimuli [24].
In contrast, Gaussian white noise signals exhibit no phase correlations and are completely characterised by first-and second-order statistics. Because white noise stimuli lack higherorder correlations present in natural scenes, linear encoding ensures that photoreceptors do not generate spurious higher-order correlations that the fly brain would use to distinguish environmental features.
Experimental and modelling studies investigating the role of photoreceptors in the detection of moving point objects [25][26][27] provide further evidence that photoreceptors contribute to the enhancement of image features that are essential for object recognition. The photoreceptors of male houseflies, for example, exhibit surprisingly large responses to moving targets [28] which cannot be explained by linear models derived from photoreceptor responses to white noise stimuli. This highlights the need to use naturalistic-like stimuli with a higher-order statistical structure, in order to characterize the nonlinear encoding mechanisms of photoreceptors [29].
More recent work [16] has shown that temporal processing performed by photoreceptors of male hoverflies (Eristalis tenax) enhances not only the moving target but also relatively static features present in the background, which are important for navigational purposes. Although this work highlights the importance of nonlinear processing in target detection, the relationship between the image features and the nonlinear component of the response is not elucidated.
The current paper provides a quantitative characterization of the relationship between the statistical properties of environmental stimuli and fly photoreceptor responses. The study is based on a nonlinear dynamical model that predicts accurately the responses of individual fly photoreceptors to white noise and naturalistic stimuli, for the entire environmental range of light intensities.
Higher-order frequency response functions analytically derived from the model equations are then used to characterize the nonlinear transformations that enable Drosophila photoreceptors to encode measures of dependence between phase angles of different frequency components of the temporal stimuli, which are invariant to contrast and illumination, such as local phase congruency. We argue that these phase-related measures, which are encoded nonlinearly, may facilitate the identification of behaviourally important features in the natural scenes.
By carrying out a comparative signal-to-noise analysis of the linear and nonlinear components of the response, we show that photoreceptors are tuned to selectively improve the Signalto-Noise Ratio (SNR) of the nonlinear component of the photoreceptor response, which encodes the local phase congruency measure. This explains why the photoreceptor responses to naturalistic stimuli are significantly different from their responses to white noise stimuli.
To validate the results, we carried out electrophysiological experiments using temporal stimuli that allow us to separate the nonlinear component of the photoreceptor responses, which encode local phase congruency, directly from the measured responses. It is shown that our model correctly predicts the measured nonlinear component of the photoreceptor response, elicited by local and non-local phase correlations introduced deliberately in the synthetic stimuli.
A similar analysis carried out using recordings from blind hdc JK910 mutant flies, indicates that the nonlinear transformations underlying the detection of phase correlations in the temporal light stimuli are performed by phototransduction alone, and do not require synaptic interactions between neighbouring neurons.
The results could have important implications beyond fly photoreceptors. A similar nonlinear encoding strategy may well be implemented in the mammalian retina or in other types of sensory neurons.

Electrophysiology
Flies were prepared for intracellular in vivo recordings from blue-green sensitive R1-R6 photoreceptors according to previously described methods [9,30]. To present the temporal stimulus pattern at different light levels, we designed a computer controlled point light source with two converging light paths (S1 Fig). In each path, LED drivers with light feedback (Cairn Research, model OptoLED) assured a linear relation between the light pattern, stored on a computer and the output of high power LED's (Seoul, model Z-Power LED P4, white, 240 lm). Neutral density filters (Kodak Wratten, ND gel filters) were used to generate 5 distinct light intensity levels, from L 0 (bright) to L -4 (very dark). Only one path was active at a time. This allowed modifying the filter setting of the inactive path in real time.
Step changes of light intensity were thus achieved by switching between the two paths implementing different filter settings. Both LED drivers were carefully calibrated to produce exactly the same light output for a given reference signal and filter setting.
The amplified temporal light stimuli (inputs) and voltage responses (outputs), sampled at 2 kHz, were low-pass filtered by analogue low-pass elliptic filters (KEMO Limited, model VBF/ 23) with a 500 Hz cut-off before being used to drive the LEDs or to perform A/D conversion respectively. The measured output of the LEDs was considered to be the input to the photoreceptor. A/D and D/A conversions (12bit resolution) were performed using National Instruments A/D and D/A boards (PCI M-IO 16E4 and PCI 6713). Custom written software was used to interface the NI boards with MATLAB (Mathworks, R7.14). For modelling and analysis purposes the data was down-sampled to 400 Hz, which provided sufficient bandwidth to capture in full the photoreceptor dynamics.

Stimuli
To characterize in full the nonlinear photoreceptor dynamics, we used naturalistic input stimuli that resemble the light fluctuations these cells are subjected to in a fly's natural habitat. Stimulus sequences with an average power spectrum S(f) = 1/f typical for natural images were selected from the van Hateren stimulus collection [29].
To derive the adaptation rules over the full operating range, we concatenated stimulus sequences with different mean intensity levels as shown in Fig 1a, which allowed us to characterize the transient photoreceptor responses to step changes of mean light intensity.
To illustrate further the validity of the model, we generated a synthetic temporal stimulus by scanning line-by-line an image (Fig 1d). The stimulus was used to excite both the final photoreceptor model as well as measure in vivo, fly photoreceptor responses. The model simulation and the electrophysiological recordings were subsequently used to reassemble the corresponding processed images for comparison.
The synthetic stimulus sequence mimics what a photoreceptor would experience as a fly moves through a natural scene, containing shadows and sunlit areas. As light intensities in such a scene can vary up to 10,000-fold, we modified our relative illumination range accordingly in 5 distinct logarithmic levels (L 0 = bright to L -4 = very dark).

Model Development
The photoreceptor model was derived directly from electrophysiological recordings using the nonlinear system identification methodology [31,32] based on the NARMAX (Nonlinear  AutoRegressive Moving Average with eXogenous inputs) model, which is applicable to a wide class of nonlinear dynamical systems. The Wiener (Linear-Nonlinear) and Hammerstein (Nonlinear-Linear) models, that are routinely used in neuroscience, can be viewed as special types of NARMAX models. Notable advantages of NARMAX methodology include the fact that it does not require knowledge of the model structure and that the noise is modelled explicitly to ensure unbiased estimates of model parameters, even if the noise is not additive or white.
The photoreceptor model consists of a polynomial discrete-time NARMAX model with variable input gain complemented by a dynamic gain control model with three adaptation time-scales (S1 File and S5 Fig). The nonlinear model with an appropriately adjusted gain fully characterizes the photoreceptor dynamics for stimuli with constant mean intensity for the entire operating range. The gain control model captures the dynamic relationship between the stimulus and the gain of the photoreceptor while it adapts to changes in mean light intensity. The three control loops reflect the dynamics of different biochemical mechanisms of light adaptation (S1 File) which we characterized in our detailed biophysical model of the photoreceptor [33].

Photoreceptor Response Decomposition using the Generalized Frequency Response Functions
Compared with previous models of fly photoreceptors [8,34], the current model not only predicts the photoreceptor responses over the entire environmental range of stimuli but also models explicitly the relationship between the stimulus intensity and the dynamic gain of the nonlinear NARMAX filter. The model can also be used to characterize the nonlinear transformations performed by photoreceptors by deriving analytically the generalized frequency response functions [35][36][37] (GFRFs) (Fig 2a-2c) of the system. GFRF's, which are extensions to the classical linear frequency response function, characterize the linear and nonlinear relationships between the photoreceptor's input and output spectral components.
The approach is similar to that adopted in their pioneering work by Victor and Shapley [3], investigating the receptive field mechanisms of retinal ganglion cells, and by French et al [34] and Asyali and Juusola [38], in their study of Drosophila photoreceptor. Victor and Shapley [3] applied stimuli composed of sums of sinusoids, to estimate, directly from data, the first and second order kernels of a Wiener series expansion of the response. In the case of the fruit fly, French et al [34] and Asyali and Juusola [38] used steps and white noise sequences, respectively, to elicit photoreceptor responses and fit first and second order Volterra kernels.
In contrast, here we use a dynamical model to derive analytically the GFRFs i.e. the Fourier Transforms of the kernels of the Volterra series associated with the nonlinear system (S1 File).
The main novelty of our analysis is that we use the GFRFs to compute spectral and temporal decompositions of the response, which allow us to provide an analytical interpretation of the role played by the nonlinear transductions at photoreceptor level rather than in the neurons downstream of photoreceptors. This provides a unique insight into the nonlinear encoding algorithms implemented by the fly photoreceptors.
The second-order GFRF, H 2 (jω 1 , jω 2 ) (Fig 2c) can be used to evaluate the magnitude and phase of the second-order output spectrum of the system Y 2 (jω) at a frequency ω, in response to all pairs of input frequencies Uðjo 1 Þ ¼ jUðjo 1 Þje jo 1 t ; Uðjo 2 Þ ¼ jUðjo 2 Þje jo 2 t satisfying ω = ω 1 ± ω 2 . In general, the n th -order frequency response function describes the contributions made by combinations of n input frequencies to the n th -order output spectrum. The total Output frequency response decomposition elucidates differences in photoreceptor coding white noise and naturalistic stimuli. a, Block diagram illustrating the output frequency response decomposition approach. b,c, Plots of the magnitude and phase for the first-and second-order frequency output spectrum of the photoreceptor, subject to an arbitrary stimulus, is given by [39].
where Y n (jω) is obtained by integrating the contributions from all possible combinations of n input frequencies satisfying By applying the inverse Fourier transform, we obtain the equivalent time-domain decomposition of the total system response into first-order (linear), second-and higher-order responses. Eqs (1-3) provide the key to elucidate the role of nonlinear transformations at the photoreceptor level. In our case, given a naturalistic stimulus sequence the total response of the photoreceptor y (t) = y 0 + y 1 (t) + h.o.t. can be approximated just by the first and second-order responses; the relative mean squared error introduced by ignoring the higher order terms is~1.3% for the bright intensity level L 0 and less than 4E-3 for levels L -1 , L -2 , L -3 .

SNR Computations
The photoreceptor response decomposition derived earlier allows us to evaluate separately the improvement in the Signal-to-Noise Ratio (SNR) of the linear and nonlinear components of the response, relative to the SNR of the noisy stimulus incorporating edges. One would expect that the photoreceptor processing selectively enhances the phase-related measure of feature significance, which are encoded nonlinearly.
Given the input signal where u s (t) is the feature-rich pulse sequence and u n (t) is a white noise sequence, the signal-tonoise ratio of u(t) is defined as Slices through the second-order magnitude and phase functions are taken along the integration lines given by f 1 + f 2 = 1Hz and f 1 + f 2 = 31Hz. d, Polar plots and distributions of the Fourier components U(jω) of the white noise (black) and naturalistic phase (red) stimuli before and after weighting by H 2 (jω 1 , jω 2 ) for pairs of input frequencies (ω 1 , ω 2 ) satisfying ω 1 +ω 2 = ω' = 2π rad/s. The total output frequency component Y 2 (jω') for white noise (black) and naturalistic (red) stimuli obtained by integrating H 2 (jω 1 , jω 2 ) U(jω 1 ) U(jω 2 ) along the line ω 1 +ω 2 = ω' = 2π. e, White noise (top, black line) and naturalisticphase (bottom, red line) stimuli. f, Output magnitude spectra |Y 1 (jω)| and |Y 2 (jω)| of the first-and secondorder responses, y 1 (t) and y 2 (t) to white-noise (top, black line) and naturalistic-phase (bottom, red line) stimuli. g, The second-order component y 2 (t) of the model response y(t) for white noise (top, black) and naturalistic-phase (bottom, red line) stimuli. where Pðu s;n Þ ¼ The noise-free photoreceptor output y s (t) is the response to the noise-free pulse sequence u s (t). The output distortion introduced by the noise is defined y n ðtÞ ¼ y sþn ðtÞ À y s ðtÞ where y s+n (t) is the response to u s+n (t).
To characterize the signal-to-noise properties of the first and higher-order responses, we decompose y s (t) and y n (t) into first and second-order responses and compute The photoreceptor's noise reduction performance was characterized in terms of the SNR improvement factor which was computed as an average over 50 repetitions. The 'noise-free' input u s (t) consists of a sequence of positive and negative pulses of amplitude +/-4.75, each lasting 50ms, with a delay of 200ms, superimposed on a constant level of background illumination L 0 . To account for the optics of the photoreceptors lens, the pulses were smoothed by convolving the signal with a Gaussian function (5ms standard deviation), accounting for the sensory neurons receptive field properties.

Photoreceptor Model Predicts Responses to Arbitrary Stimuli over the Environmental Range of Light Intensities
To demonstrate that the estimated model provides an accurate representation of R1-R6 photoreceptors for the entire environmental range of light intensities, the model was validated extensively using intracellular recordings from the photoreceptors of different flies (S1 File). The model predictions match remarkably well the experimental responses to naturalistic (S4 Fig) and white noise stimuli (S6b Fig) as evidenced by the relative mean squared prediction errors summarized in S1 and S2 Tables respectively. To further illustrate visually the prediction accuracy, we generated an input time series by scanning line by line a naturalistic image (Fig 1d). This temporal stimulus was used to stimulate photoreceptors and measure in vivo their responses that were converted back to an image. Fig 1d, shows the original image side-by-side with the images reconstructed from the experimental and model predicted photoreceptor response time series.
The second-order frequency response function explains the difference in encoding naturalistic and white noise stimuli To investigate the link between the phase structure of the stimulus and the response, we derived a synthetic stimulus using the magnitude of a Gaussian white noise stimulus and the phase spectrum of a naturalistic stimulus sequence with 1/f magnitude characteristic. The analytical GFRFs derived for our model were used to compute and compare the linear and second-order responses y 1 (t) and y 2 (t) to the original Gaussian white noise (GWN) stimulus (Fig 2e top) and to the modified GWN stimulus with naturalistic-phase (Fig 2e bottom).
In line with previous experimental studies [30], we found that |Y 1 (jω)|, the magnitude spectrum of the linear component of the photoreceptor response (Fig 2f-inset panels), dominates the magnitude spectrum of the nonlinear component |Y 2 (jω)| (Fig 2f-main panels) for GWN as well as naturalistic stimuli. However, whilst |Y 1 (jω)| essentially remains unchanged, there is a marked increase of |Y 2 (jω)| for naturalistic stimulus compared to the GWN case (Fig 2fmain panels). This is also reflected in the time domain. The amplitude of the second order component y 2 (t) of the response to the modified GWN stimulus (Fig 2g bottom) is significantly larger compared to the second order component corresponding to the original GWN stimulus (Fig 2g top). Specifically, whilst the variance of the y 2 (t) component of the response to the random-phase stimulus represents~2% of the total response, for the naturalistic-phase stimulus with the same magnitude spectrum, the variance of the y 2 (t) component increases to~50% of the variance of the total response. This increase is entirely due to the non-random structure of the phase spectrum. For the white noise stimulus the Fourier components Uðjω k Þ ¼ jUðjω k Þje jθ ω k have phase angles θ ω k that are uniformly distributed in the range [0, 2π) (Fig 2d top). Because the phase of H 2 (jω 1 , jω 2 ) is remarkably flat for frequencies satisfying ω i + ω j = ω (see for example the phase slice along f 1 +f 2 = f = 1 Hz and f 1 +f 2 = f = 31 Hz shown in Fig 2c), the phase-angles of the 'weighted input frequencies' H 2 (jω 1 , jω 2 )U(jω 1 )U(jω 2 ) satisfying ω i + ω j = ω remain uniformly distributed. Consequently, the magnitude of the second-order frequency spectrum of the response Y 2 (jω), defined in Eq (2) for n = 2, will be small as shown in Fig 2f (top). In essence, the uniform distribution of phases of H 2 (jω 1 , jω 2 )U(jω 1 )U(jω 2 ) means that for every complex vector H 2 (jω 1 , jω 2 )U(jω 1 )U(jω 2 ) with phase φ there is a vector H 2 (jω 0 1 , jω 0 2 )U(jω 0 1 )U(jω 0 2 ) with similar amplitude but opposite phase φ' = φ±180; these vector pairs tend to cancel out when the integral Eq (2) is computed. As a result, the corresponding second order temporal response y 2 (t) is small as seen in Fig 2g (top). The second-order response is not zero because the magnitude of H 2 (jω 1 , jω 2 ) is not constant along the constant frequency lines ω i + ω j = ω, as shown in Fig 2c (magnitude slices along f 1 +f 2 = 1Hz, for example).
For the modified GWN stimulus with 'naturalistic' phase spectrum and white noise magnitude spectrum, the phases of the frequency components U(jω i ) and the phases of H 2 (jω i , jω j )U (jω i )U(jω j ) are not distributed uniformly and as a result, the magnitude of Y 2 (jω) does not cancel out (Fig 2d bottom panels).
The results show that the photoreceptor is sensitive to the phase structure of the temporal stimuli, specifically to correlations between the phases of different frequency components of the temporal stimulus. The analysis shows that the photoreceptor responds linearly to white noise and nonlinearly to naturalistic, feature rich stimuli because of the particular shape of the phase function associated with the second-order frequency response function.
Fly photoreceptors are tuned to encode robustly temporal local phase congruency Spatially-localised features, such as edges, are ubiquitous in natural scenes. Edges correlate with contours and textures in a natural scene, which form the basis for higher-level visual processing tasks such as motion detection and object recognition. Not all edge features are characterized by sharp changes in luminance at the boundary of an object. Often edges are blurred or soft-edged like those of shadows [40]. However, all edge-like features exhibit high local phase congruency i.e. the phases of the constituent Fourier components are aligned.
The previous frequency response analysis predicts that points of high local phase congruency, where the difference in phase between different frequency components is zero will elicit very strong second-order responses in fly photoreceptors since in this case the magnitude of Y 2 (jω) is just the integral of the magnitude of H 2 (jω 1 , jω 2 )U(jω 1 )U(jω 2 ) for ω i + ω j = ω i.e. no cancellation occurs.
To test this hypothesis, we designed a synthetic stimulus consisting of a sequence of narrow square pulses (S1 File) superimposed on a Gaussian white noise sequence with mean L 0 (bright stimulus). The variance of the stimulus was designed such that the output of the automatic gain control module is essentially constant i.e. gain adaptation plays no role in these experiments. The photoreceptor responses to this stimulus, predicted by our model, were found to match closely (S3 Table) the in vivo intracellular recordings made from the fly photoreceptors using the same stimulus sequence (S6 Fig). As before, the model responses were decomposed into linear and nonlinear (second-order) components according to Eqs (1-3) using the GFRFs derived for the NARMAX model with the constant gain corresponding to the 'bright' mean intensity level L 0 .
As seen in (Fig 3a and 3b), the pulses are hard to distinguish from the background noise in the synthetic stimulus and the linear component of the model response. In contrast, the nonlinear component of the response encodes their location quite precisely by large negative peaks (Fig 3c). Even a simple threshold decoder can be used to extract the encoded 'message'-position of the pulses-from the nonlinear component of the response. A similar decoder applied to the linear response would generate a significant number of false positives.
Remarkably, the nonlinear response appears very similar to the output (Fig 3d) of an established algorithm for computing local phase congruency [15], which is widely used to detect edges in computer vision. Given that the phase congruency measure is invariant to changes in intensity and contrast, it provides arguably the most sparse and efficient representation for edge-and line-like features [14].
A major practical implementation issue is that, being a normalized quantity, phase congruency is highly sensitive to noise. Although the algorithm used to compute phase congruency implements noise reduction techniques [15], it does not detect all the real pulses (Fig 3a) whilst spurious detections still occur (see missing or extra red 'peaks' in Fig 3d compared with pulse locations indicated in Fig 3a).
From this point of view, nonlinear transductions at photoreceptor level encode robustly local phase congruence because the nonlinearity is tuned to reject white noise signals. To demonstrate this, we compute the SNR improvement factors (S1 File) for the linear y 1 (t) and nonlinear y 2 (t) components of the response to stimuli consisting of a pulse sequence with added white noise having different variance levels. As seen in Fig 3g, for y 2 (t) the SNR improvement factor Q 2 is significantly higher (almost five fold improvement) than Q 1 computed for y 1 (t) (two fold improvement).
The nonlinear response encodes robustly the phase correlations buried in noise because the phase of H 2 (jω 1 , jω 2 ) is almost constant along the integration paths ω 1 + ω 2 = ω, which ensures that the phase shift introduced by the second-order frequency response function is independent on the input frequencies.
Artificially changing the phase of the second-order frequency response function, makes the second-order response noisier, reduces the amplitude of the response around the steps and introduces spurious peaks in places where the local phase congruency is low, as seen in Fig 3e. This provides strong evidence that the nonlinear transductions in fly photoreceptors are optimized to enhance behaviourally important higher-order statistical correlations in the natural scenes whilst being largely insensitive to random-phase stimuli.
The reverse-engineered algorithm implemented by photoreceptors is remarkable for its simplicity and, to the best of our knowledge, provides an entirely new approach for computing local phase congruency. While state-of-the art conventional algorithms based on wavelet filter banks are complex, computationally expensive and sensitive to noise [15], the photoreceptor algorithm implements a single nonlinear filtering operation to encode local phase congruency.
To illustrate the practical applicability of the photoreceptor-inspired edge detection algorithm, we computed edge maps for the image shown in Fig 3h by applying a threshold decoder to the standard local phase congruency map of the image (Fig 3h) and to the nonlinear component of the photoreceptor responses to time-series of pixel intensity values along each line of the image. Visually at least, the edge map generated using photoreceptor algorithm (Fig 3m) is 'cleaner' than the edge map generated using the standard local phase congruency algorithm (Fig 3l).
Fly photoreceptors encode non-local phase correlations between the spectral components of the input Natural images exhibit not only local but also global phase correlations. It has been argued that both local and global higher order statistics of natural images play an important role in texture Fly Photoreceptors Encode Phase Congruency or symmetry discrimination [41,42]. To test the sensitivity of photoreceptors to non-local phase correlations we designed a synthetic stimulus (Fig 4a), which exhibited quadratic phase coupling (QPC) at 10Hz (i.e. phase(f 1 ) + phase(f 2 ) = phase(f 3 ) where f 1 + f 2 = f 3 = 10 Hz).
Specifically, the QPC stimulus was constructed by computing the Fourier spectrum of a Gaussian white noise signal, modifying the phases of the spectral components to satisfy the above conditions whilst keeping the magnitude function unchanged and finally applying the inverse Fourier transform (S1 File). As seen in Fig 4d, the resulting phase-modified signal has the same Fourier magnitude spectrum as the original white noise signal. However, whilst the phases of the QPC input frequencies U(jw) are still uniformly distributed between 0 and 2π, they are clearly correlated as seen in Fig 4b. The responses of the photoreceptor model to the white noise and QPC stimuli were decomposed into linear-and second-order responses (Fig 4c). Subsequently, the Fourier spectrum was computed separately for each component of the photoreceptor response.
Because the linear response is not sensitive to the phase structure of the stimulus, the Fourier magnitude spectra of the linear responses to the two stimuli sequences are identical (see Fig 4e and 4b). In contrast, as expected, the second-order response to the QPC stimulus shows a significant increase in the magnitude of the 10 Hz output frequency compared with the second-order response to the white noise stimulus.

Experimental Validation of the Photoreceptor Model Predictions
The higher order visual processing neurons have to extract the nonlinear codes generated at photoreceptor level from the overall responses. A simple approach to extract the second-order response to a given stimulus is illustrated in S9 Fig. Essentially, the sum of even-order responses to a given stimulus is the average between the photoreceptor response to the stimulus and the response to the out-of-phase (inverted) version of the stimulus. Since in our case the higher-order responses greater than two are negligible, this method generates the second-order component of the photoreceptor response. Using this approach it was possible to demonstrate experimentally that the nonlinear computations performed by the photoreceptor are indeed those predicted by the model-based analysis.
To extract the second-order responses to a given temporal stimulus directly from the experimental recordings we constructed stimulus sequences by alternating the original stimulus sequence with its inverted version (S1 File). The experimental responses to the two versions of the stimuli were averaged and compared with the model predictions. Experiments were carried out using white noise stimuli, stimuli consisting of a pulse sequence superimposed on white noise as well as stimuli exhibiting quadratic phase coupling at 10Hz. Fig 5a shows that the second-order components of the photoreceptor responses (seven repetitions) to the stimulus consisting of the pulse sequence superimposed by Gaussian white noise, extracted directly from the experimental recordings and the model predicted y 2 component.
As predicted, the nonlinear response is significant around points of maximum local phase congruency; the amplitudes of the negative peaks in the nonlinear response at the location of the pulses represent more than 25% of the corresponding peak amplitudes of the total response. The close match between the model predicted and the experimentally derived y 2 component around the negative excursions triggered by the embedded pulses, is further illustrated in S10 Fig. The linear-and the second-order components of the response account for 91% and~9% of the overall variance of the total response, respectively. The prediction error variance corresponding to the y 2 component extracted directly from experimental data (S10 Fig) represents~14% of the total y 2 variance.
On the other hand, the magnitude spectrum of the y 2 component of the photoreceptor response to the QPC stimulus (Fig 5b) shows a~5 fold increase in magnitude at 10Hz compared to the magnitude spectrum of the y 2 component of the response to the original GWN stimulus, as predicted by model.

Fly Photoreceptors Encode Phase Congruency
Given that the photoreceptor model was derived using experimental recordings from a different fly, these results demonstrate further the validity of our model and, more importantly, that the nonlinear encoding of phase correlations is a generic information processing strategy of fly photoreceptors.

Investigating the role of the retinal network
Photoreceptor responses in wild flies are modulated by feedback from two classes of interneurons, i.e. large monopolar cells (LMC) and amacrine cells (AC), and from axonal gap-junctions, which pool the responses from six photoreceptors [43][44][45]. It is therefore natural to ask to what extent the processing capabilities demonstrated earlier are due to temporal processing by the photoreceptor alone. To elucidate this question, we measured photoreceptor responses to the original multi-level naturalistic stimulus in blind hdc JK910 mutants [46,47] that lack histamine in their photoreceptors. Fly photoreceptors use neurotransmitter histamine to communicate visual information to interneurons [48]. In the histamine deficient mutants, the lamina interneurons fail to receive and transmit visual information and their feedback synapses can no longer modulate photoreceptor output [47]. Essentially, these mutant flies are blind. By comparing intracellular recordings from photoreceptors of wild type flies to those of hdc JK910 mutants, one can test how the lamina network affects adaptation and information processing in photoreceptors.
As seen in S11 Fig, mutant photoreceptor responses have dramatically reduced contrast sensitivity for bright (L 0 ) stimuli and their capability to quickly adapt to the mean illumination is significantly impaired. However, the light adapted responses of histamine photoreceptors to naturalistic stimuli having a mean luminance L-1 are very similar to wild-type responses. As we are interested to investigate the nonlinear properties of isolated photoreceptors exhibiting a normal response range, we inferred and validated (S12 Fig) a mutant photoreceptor model using measured responses to naturalistic stimuli sequences with mean luminance L-1 (S1 File). The second-order GFRFs (Fig 6a), derived for this model, are very similar to those of wild-type photoreceptors. In particular, the phase constancy along integration lines is preserved in mutant flies. As a consequence, the mutant photoreceptor responses to the two classes of synthetic stimuli are very similar to the wild type responses (Fig 6b), providing strong evidence that the nonlinear transformations, underlying the detection of high-order phase correlation in the temporal light patterns, are performed by the photoreceptor alone, independently of neighbouring neurons.

Discussion
We have demonstrated that R1-R6 photoreceptors perform nonlinear transformations that encode biologically relevant, higher-order statistical features that are represented in the Fourier phase spectrum of temporal stimuli. In particular, we have shown that photoreceptors encode points of maximum local phase congruency, which occur at the location of an edge or line, as well as long-range phase correlations, which characterize symmetry and texture properties of natural images [41].
An important conclusion of our analysis is that the nonlinear transductions in fly photoreceptors are tuned to maximize the response to combinations of spectral components that are congruent in phase or are phase-coupled and to minimize the response to temporal stimuli with a random phase spectrum. This ensures that the nonlinear coding is robust to noise and explains why photoreceptors respond linearly to non-informative white noise stimuli and nonlinearly to naturalistic time-series exhibiting local and global phase correlations between their spectral components. It also explains why the models derived using responses to white noise stimuli fail to capture the key nonlinear transformations performed by photoreceptors.
This strategy for processing temporal stimuli does not require dynamic adaptation to stimuli statistics beyond the mean and variance. Specifically, we show that what appears to be an adaptation of the photoreceptor to the statistical structure of the stimuli is in fact explained by the shape of the two-dimensional phase function corresponding to the second-order GFRF of the photoreceptor, which is almost constant along the lines ω i + ω j = constant.
From the point of view of information theory, different encoding of naturalistic and white noise signals can be viewed as a solution to the problem of matching the stimulus (source) with the communication channel in a probabilistic sense in order to achieve an optimal trade-off between two competing goals: minimizing distortion in decoding behaviourally relevant stimuli features and minimizing the information rate, that is, the energetic costs associated with phototransduction [49].
The fact that the second-order phase function computed for a separate photoreceptor model, derived for Caliphora, exhibits similar characteristics (see S8 Fig) to those of Drosophila, suggests that the maximization of sensitivity to phase aligned or coupled frequency components is a fundamental 'design principle' of fly photoreceptors, which may well apply to other sensory neurons. Previous studies of the primary visual cortex [42,50] have postulated the presence of nonlinear mechanisms that are sensitive to phase correlations. Here we demonstrate for the first time that in Drosophila these mechanisms operate at the photoreceptor level. The main benefit of implementing nonlinear encoding at photoreceptor level is that it may facilitate the efficient decoding of key stimuli features such as edges, at higher visual processing levels. Essentially, as we have illustrated earlier, edge maps could easily be extracted from the nonlinear component of the response using a simple threshold decoder. Given that neural circuits are highly optimized, one would expect that in the absence of the nonlinear photoreceptor code the higher visual processing stages would have to implement additional computations, leading to more complex downstream neural circuit architectures. At the same time, encoding temporal edges at photoreceptor level, before the information from six photoreceptors has been pooled by the interneurons, should help improve further the signal-tonoise ratio of the encoded features.
The similarities that exist between responses of primate cones and blowfly photoreceptors [51] suggest that nonlinear transformations performed by cone photoreceptors ultimately achieve the same processing goals, albeit using different molecular mechanisms and signal processing steps.
Previous experimental and theoretical studies of early visual processing in humans indicate the existence of detectors that are highly sensitive to features characterized by high phase congruency [52,53]. We speculate that human photoreceptors implement similar nonlinear processing of the visual stimuli to detect phase congruency, which could help explain why neurons in the primary visual cortex can reliably signal phase congruence and how the phase congruency information is extracted from the visual stimuli. Since moving spatial edges generated by saccadic eye movements or by moving objects generate temporal edges at photoreceptor level, selectively encoding these temporal features and enhancing their salience, should facilitate downstream processing of the spatio-temporal visual stimuli for edge and motion detection [54].
The simple technique we used to separate the second-order response directly from experimental recordings could easily be used to test this hypothesis for mammalian retinal cones.
One could envisage a simple model where eye saccades map localized spatial edges onto temporal edges that are encoded and enhanced by photoreceptors, enabling the downstream neural circuits to use timing in addition to spatial information to detect edges and group them into contours.
As it is not clear how the higher-order spiking neurons could implement efficiently this processing step, it could be that all downstream visual processing relies on the phase congruency information generated by photoreceptors to the extent that the absence of this information may incapacitate downstream feature detectors. If this were true, applying such nonlinear transformations to the visual stimuli, prior to delivering these to retinal or ganglion cells, may improve significantly the performance of artificial retinas [55,56].