Machine learning discovers numerous new computational principles supporting elementary motion detection

Domain: chemistry
Published: 2026-03-03
Editorial outcome: major_revision

Abstract
Motion direction detection is a fundamental visual computation that transforms spatial luminance patterns into directionally tuned outputs. Classical models of direction selectivity rely on temporal asymmetry, where motion detection arises through either delayed excitation or inhibition. Here, I used biologically inspired machine learning applied to retinal and cortical circuits to uncover receptive field architectures capable of direction selectivity. These include mechanisms based on asymmetric synaptic properties, spatial receptive field variations, new roles for pre- and postsynaptic inhibition, and previously unrecognized kinetic implementations. Conceptually, these circuit architectures cluster into eight computational primitives underlying motion detection, four of which are previously undescribed. Many of the solutions rival or outperform classical models in both robustness and precision, and several exhibit enhanced noise tolerance. All mechanisms are biologically plausible and correspond to known physiological and anatomical motifs, offering fresh insights into motion processing and illustrating how machine learning can uncover general principles of neural computation.

Introduction
Visual sensation begins with photoreceptors, but they do not encode motion information. To perform motion computations, the visual system must transform spatially separated luminance signals to an output that is stronger in specific directions 1-4. Elementary motion detectors, composed of shallow networks with directionally untuned inputs, represent the simplest neural substrate that can support the computation of direction selectivity (DS). Over 60 years ago, two foundational studies proposed temporally asymmetric circuitry as the basis for DS in elementary motion detection circuits. The core computation of the Hassenstein & Reichardt (H&R) detector is delay-and-compare, where signals from spatially offset inputs are multiplied after the input activated first in the preferred direction is delayed by a low-pass filter 5. Similarly, Barlow and Levick's (B&L) model compares fast excitatory signals with delayed inhibition from the "null" side -- the position from which the null direction stimulus enters the cell's receptive field (RF) -- to suppress responses to motion in that direction 6. While these classical models have been extended and refined by more recent research, revealing several distinct biological implementations, the core computational motifs of kinetic filtering applied to preferred side excitation or null side inhibition remain central to our understanding of motion detection 2,3,7-28.

The anatomical simplicity of elementary motion detectors belies their potential for remarkable functional diversity. The neurons in the input layer can process sensory stimuli in complex, nonlinear ways, depending on their filter characteristics and the statistical properties of the stimulus 29,30. Accordingly, past research identified several distinct solutions to motion detection in biological circuits mediated by asymmetric filtering of sensory input through dissimilar synaptic dynamics and receptive field properties 8,9,22-24,31-37. However, it remains unclear whether additional circuit-based solutions exist and to what extent known algorithmic implementations of directional selectivity encompass all potential operational principles of elementary motion detectors.

In this study, I employed machine learning to systematically explore how functional architectures of spatial and temporal filters of RF center and surround components are capable of supporting motion detection. Exploration of the vast combinatorial space formed by variability in individual RF components and their interactions revealed multiple distinct circuit configurations capable of performing motion computations with multiple asymmetries in the RF composition capable of leading to DS levels that were often comparable to those observed in the classical models. Notably, despite structural differences among these circuits, their information processing consistently converged on a small set of algorithmic implementations, which I refer to as ‘computational primitives’. These included classical models like H&R and B&L, as well as previously unrecognized circuit motifs. Examination of these algorithmic implementations showed they can drive robust direction selectivity and may sustain enhanced noise tolerances. These findings highlight the power of machine learning to uncover previously unrecognized solutions and provide algorithmic-level insights into neural computations mediated by simple networks that often elude human intuition.

Results
Machine learning framework to study motion computations
The ability to detect motion direction was examined using feedforward, shallow network architectures. In the retina, previous studies have shown that excitatory inputs from bipolar cells to direction selective ganglion cells (DSGCs) are organized into distinct populations arranged sequentially along the cell’s preferred motion axis 38. In this model, stimulus activation was first filtered through a linear spatial RF, where the overlap between the stimulus and the RF shape determined the potential for graded transmitter release. Temporal filtering, together with readily releasable pool (RRP) dynamics, then governed the actual rate of synaptic drive onto the postsynaptic neuron (see Methods for more details ). To evaluate the motion detection capacity of this architecture, I conducted simulations in which four excitatory input populations -- each potentially having distinct receptive field properties -- were spatially distributed along an arbitrarily defined preferred direction axis onto a multicompartmental DSGC model (Fig. 1a). Only excitatory inputs were included in this initial set of simulations.

Consistent with biological observations that most bipolar cell terminals lack directional tuning 39, all presynaptic inputs in these models were set to produce similar responses to stimulation in different directions. Directional selectivity of the DSGC was quantified using a direction selectivity index (DSI), computed as the vector sum of peak subthreshold membrane potential at the soma across five stimulus velocities (Fig. 1a) 40,41. While subthreshold responses generally show weaker directional tuning than spiking output 33,36,42-48, this approach allowed robust detection of DS mechanisms across a range of stimulus intensities without the confounding effects of spiking nonlinearities (Supplementary Fig. 1).

Machine learning was implemented as follows: starting with randomly initialized models, the fittest model (determined by the highest directional performance) seeded the next generation. The parameters describing the characteristics of presynaptic RFs were adjusted using a genetic algorithms-based machine learning paradigm 49,50. In all cases, including the seed values, parameters were set to range within biologically relevant boundaries. In most simulations, constraints were put on the difference in RF description of the presynaptic inputs. Some parameters were allowed to vary independently across presynaptic groups. To explore the parameter space beyond the constraints of the initial seed values, other parameters were 'shared' or mutated jointly for all presynaptic cells. Overall, this approach enabled the identification of optimal configurations for directional tuning among presynaptic populations for various aspects of their RF structure.

To validate the effectiveness of this search strategy, I tested the ability of the algorithm to find the H&R solution for optimal spatiotemporal distribution of units with distinct kinetics. Evolved simulations (n = 100 independent runs) converged on a singular outcome where response dynamics became progressively faster with increasing distance from the preferred side of the DSGC, leading to synchronized peak responses and robust postsynaptic summation in the preferred direction (DSI = 23.8 ± 1.6%, Fig. 1b), as postulated by the H&R model 5,26,38.

Direction selectivity computation with differences in presynaptic receptive field sizes and orientations
The spatial RF properties can influence motion processing and impact signal transformation in H&R model families 31,32,49, but could they mediate DS computations in the absence of different temporal filters? To investigate whether and how inputs that only differ in spatial RF properties could compute DS, I evolved networks with presynaptic groups constrained to have identical kinetics but potentially variable extent of their RFs (Fig. 1c). In this configuration, all presynaptic neurons produced identical responses when stimulated by a full-field stimulus. However, the representation of moving bars was dissimilar, with larger RFs activated earlier, for longer duration and smaller amplitude, reflecting the slower engagement of the RF by the stimulus 49,51,52. Evolved models demonstrated that differences in spatial RF properties can lead to strong directional tuning (DSI = 22.1 ± 0.8%) and that the optimal functional circuit configuration that supported DS was formed with progressively increasing RF diameters in the input population along the motion preference axis (Fig. 1 c, c ).

The spatial arrangement of the presynaptic RFs in the optimal solution offers an insight into how synchronized engagement of the presynaptic populations is achieved. RFs from the different input populations converge near the distal region of the postsynaptic dendritic arbor at the location where the preferred direction stimulus first encounters the cell (Fig. 1c ). This near-alignment facilitates motion computations that begin substantially earlier than in the kinetic (H&R-like) implementation (Fig. 1b-e, Supplementary Fig. 1) 53, potentially compensating for delays in retinal circuitry, a property that is important in circuits where speed of motion detection is paramount 54-56.

Spatial differences between presynaptic cells are not limited to size differences. Indeed, not all RFs are perfectly round, and orientation-selective neurons with elliptical RFs are frequently observed in DS circuits 36,57-61. Implementing models with non-circular RFs revealed that DS computation can take advantage of presynaptic populations that vary only in their orientation. In the optimal configuration, the long RF axis of the presynaptic population transitioned from perpendicular to parallel along the preferred direction of motion. This arrangement of orientation-selective cells closely mirrored the RF size-based solution in terms of directional selectivity levels (DSI = 21.2 ± 0.7%), velocity tuning, and detection time (Fig. 1d-e).

Both the kinetic and spatial solutions described above produced comparable selectivity to bars of varying widths and supported strong directional spiking responses when sodium and potassium channels were added to the postsynaptic soma (Supplementary Fig. 1). From an algorithmic perspective, organizations of the presynaptic spatial RFs promoted an 'anti-H&R' computation where neurons with slower motion responses, mediated by larger RF sizes or orientation of the long axis of an oval RF parallel to the preferred motion direction, innervated the null side of the postsynaptic cell and thus were recruited last in the preferred-direction activation (Fig. 1c-e).

To confirm that DS was driven by asymmetric presynaptic circuitry rather than postsynaptic voltage effects 62-70, I examined directional tuning in models where all input groups had identical properties. In these cases, the DSI of the best-performing models was only 2.4 ± 0.1%, indicating a minimal contribution of postsynaptic processing to the computation (Fig. 1e, ‘symmetric RF’).

Mechanisms of direction selectivity computation with presynaptic surrounds
The investigation of how presynaptic RF architecture influences DS computation would be incomplete without considering the effects of surround suppression -- a canonical property of most visually active cells 71 -- that is known to impact motion computations 14,72,73. In the models presented below, the surround component was implemented using a separate spatiotemporal filter, spatially aligned with the center RF. Surround activation reduced the net synaptic drive: the full RF signal was computed as the center minus surround conductances and then rectified to prevent negative synaptic conductance. To systematically examine the role of presynaptic surrounds in DS, I simulated inputs with varying surround strength, kinetics, and spatial properties, including size and shape (Fig. 2).

Consistent with the larger number of parameters describing the center-surround RF, these surround models generated diverse solutions, which fell into two main categories: circuits having preferred-side inputs with strong and fast surrounds that delayed and narrowed the full RF responses (Fig. 2a, b, and c), or organization of presynaptic cells with stronger and slower surrounds on the null-side portion of the postsynaptic dendritic arbor (Fig. 2a , b , and c ).

Functionally, although asymmetric surrounds could promote DS through H&R and anti-H&R computations, many solutions did not produce motion responses with sufficiently distinct kinetics, which is the cornerstone of these computations (e.g., Fig. 2c ). Instead, across many circuit configurations, surround suppression strategically blocked responses at specific time points to temporally align the peak waveforms of the full RF across all presynaptic inputs (Fig. 2b, c, and c ). Importantly, this temporal alignment occurred only during motion responses (Fig. 2d), as RF activation dynamics were shaped by the differential spatiotemporal engagement of center-surround components with the moving stimulus.

Overall, in models with variable surrounds, kinetic models achieved the highest DSI levels for moving bar stimulation (17.5 ± 1.4%, Fig. 2e) and for drifting gratings (Supplementary Fig. 2). However, combination of variable center and surround components produced stronger directional tuning (DSI = 28 ± 1.3%) in models with variable center and surround diameters compared to those with variable center and surround kinetics (DSI = 23.4 ± 0.9%, Fig. 2e), the latter being close to the performance of models with free spatiotemporal centers (DSI = 24.6 ± 1%). These spatial solutions matched the directional performance of models with unconstrained center-surround properties (DSI = 27.9 ± 1%, Fig. 2e).

Direction selectivity with synaptic weights
Algorithmic implementations of DS computation rely on asymmetric processing across different information channels that feed into the detector 1. Prior work suggested that asymmetry in synaptic weights between the input population and the detector may be sufficient to promote directional capabilities 74. To better understand the factors enabling a circuit with different synaptic weights to identify the sequence of presynaptic activation, I first examined this question in a minimal model that featured an input layer composed of two cells with identical filters that summed linearly at the detector (Fig. 3a). Already in this simple scenario, the summation of two inputs with similar kinetics but different amplitudes was found to be directionally asymmetric. The maximal signal was reached with sharply rising, slowly decaying responses when the weaker input was recruited first (Fig. 3a). In some models, an alternative solution with slower kinetics and reversed weights was reached, but the directional selectivity of this solution was about half of the optimal model (DSI = 6% vs. 11%, Fig. 3a, inset).

Encouraged by these results, I pursued the examination of weight-based solutions in the full DSGC model, probed at multiple velocities and directions (Fig. 3b, c). As its minimal counterpart, this biologically grounded circuit could be trained to generate directional preference, with the optimal solutions being a band-like organization of stronger synaptic inputs perpendicular to the motion preference axis (Fig. 3b). As in minimal models, two solution types were apparent. About half of the models evolved to have two bands of stronger synapses, with the band on the null side of the DSGC having more numerous inputs, mirroring the optimal solution seen in minimal models. In other cases, double strong-to-weak transitions evolved, each following the logic of the alternative minimal model solution (Fig. 3b, right). Directional tuning was observed in all simulations across multiple velocities; however, as expected, it was less robust compared to variable-RF models (DSI = 7.4 ± 0.4%, Fig. 3b, c). Importantly, however, this result indicates that differential spatiotemporal filters are not required for the detection of motion direction, demonstrating a potential for synaptic weight differences in implementing DS computations.

Computational primitives in elementary motion detectors
To identify additional DS algorithms in elementary motion detection circuits, I modeled various configurations of excitatory cells where multiple RF and synaptic parameters were free to vary between the input groups (Supplementary Figs. 3, 4). Critically, all solutions converged on one of the four computational primitives revealed in this study, namely H&R, anti-H&R, amplitude, and temporal alignment (Fig. 4a). The probability of observing the computations in different functional circuits is illustrated as a two dimensional plot, in which the H&R and anti-H&R implementations, which rely on temporal differences in motion responses, form the horizontal axis. The vertical axis indicates the degree of similarity of normalized response waveforms in the preferred direction (temporal alignment) or at the time of stimulus onset (amplitude). Some solutions could be explained with a mix of primitives, for example, by forming H&R-like slower signals that are also delayed in time on the preferred side of the detector, as was often seen in circuits with varying surrounds (Fig. 4b).

This examination of the parameter space revealed that amplitude-based computational primitives yielded the lowest selectivity. This mechanism was optimal only in models where the free parameter was synaptic weight or synaptic facilitation (Supplementary Fig. 3). But what algorithm produces the best directional performance? Removing all constraints on RF description illustrated that models can utilize diverse strategies to compute comparable DS levels (DSI = 27.8 ± 1%, Supplementary Fig. 4). Computational primitives provided a convenient platform for understanding signal processing in circuits with complicated filters. Across simulations, it was clear that the models evolved towards temporally aligning the responses in the preferred direction, either as a lone computational implementation or combined with H&R or anti-H&R kinetic differences, singling it as the best strategy to achieve the best selectivity without involving postsynaptic inhibition (Fig. 4c).

Minimal requirements for DS mechanisms
The RF configurations described above differ substantially from the classical H&R model, which relies primarily on temporal filtering. This raises a question: what RF properties are truly required for DS computations?

In the formulation used here, each presynaptic RF contained two nonlinear components. First, synaptic dynamics incorporated depletion of the RRP, which altered both the gain and the shape of the input–output transformation. Second, because synaptic conductances were constrained to be non-negative, they imposed a rectifying nonlinearity on the visual signal. To determine whether these nonlinearities were essential for DS, I again examined a minimal model containing only two presynaptic neurons and a linearly integrating postsynaptic detector. This simplified architecture retained the same RF definitions as the full multicompartmental model but was far easier to interpret. Using this framework, models in which center-surround RF kinetics, spatial extents, or relative amplitudes were allowed to vary continued to exhibit robust DS (Supplementary Fig. 5).

Replacing the nonlinear synaptic release mechanism with a purely linear formulation and removing signal rectification reduced, but did not eliminate, the DS produced by the spatial- and amplitude-based solutions when stimulated with moving bars or drifting gratings (Supplementary Figs. 5, 6). Thus, nonlinear RF interactions enhance, but are not strictly required for DS computations.

As shown above, center-surround interactions alone can generate DS. Indeed, spatial- and amplitude-varying models with linear synaptic release lost their DS capabilities when the surround component was not implemented (Supplementary Figs. 5, 6). Further analysis revealed that either nonlinear RF transformations or center–surround structure was sufficient to support direction discrimination for moving bars, whereas robust DS responses to drifting gratings required models in which the center and surround components were allowed to differ in temporal kinetics (Supplementary Figs. 5, 6).

Enhanced performance in elementary motion processing models with a postsynaptic inhibitory drive
The models thus far examined circuit architectures where the detector receives innervation from populations of excitatory neurons. To study how postsynaptic inhibition affects directional processing in elementary motion circuits, I introduced an additional group of input cells that processed visual stimuli independently -- with the same organizational rules as described above -- and inhibited the postsynaptic cell. As before, I first tested whether this architecture could replicate known DS computations, in this case, the classic B&L formulation. The B&L model relies on spatially offset RFs, which can be achieved by allowing the presynaptic groups to have different synaptic weights 6. Training models with variable strengths of the inhibitory drives achieved impressive directional performance (DSI = 50.8 ± 0.8%, Fig. 5a). The optimal solutions assigned a stronger inhibitory input on the null side of the dendritic arbor, leading to delayed recruitment of inhibition in the preferred direction, which is the essence of B&L computation (Fig. 5a) 37,46,75,76.

Notably, solutions to other functional circuit configurations did not always converge on B&L-like temporal interactions. Some models with kinetically diverse inhibitory cells, like the one shown in Fig. 5b, or with short-term plasticity of inhibitory synapses (Supplementary Fig. 7), were solved with an 'anti-B&L' arrangement, where early inhibition preceded delayed excitation in the preferred direction. In circuits with spatially varying inhibitory inputs, a different solution with a double-humped shape of the inhibitory profile that introduced a pause-in-inhibition around the time of the peak excitatory conductance was frequently observed (Figs. 5c, Supplementary Figs. 8, 9).

The information provided by the inhibitory drive can influence motion computations in ways that may not be reflected by the basic spatiotemporal filtering mediated by RF centers described above. Asymmetry in inhibitory RF surrounds can support impressive tuning levels (Fig. 5d-f, Supplementary Figs. 9, 10). Furthermore, while DS performance of models with only the center or surround RF diameter was allowed to vary was rather low (Figs. 5g, h), it was significantly improved when both RF components of the inhibitory cells were variable (Fig. 5i, Supplementary Fig. 10).

As with excitatory DS circuits, I found four distinct computational primitives provided by inhibitory drive (Fig. 6). These include the classical B&L, ‘anti-B&L’, and ‘pause-in-inhibition’ mechanisms described above. Similar algorithmic solutions were observed in circuits evolved to produce DS responses to drifting gratings (Supplementary Fig. 11). Unexpectedly, while moving bar stimulation did not solicit solutions with directionally tuned inhibition, this computational primitive was observed as an optimal solution with drifting gratings (Supplementary Fig. 11).

More complicated interactions emerged when excitatory and inhibitory inputs were allowed to evolve. The two communication channels in this scenario supported a rich combinatorial space of possible interactions between various RF components. In general, models with varying dynamics outperformed simulations with varying RF sizes or orientations: kinetic vs. spatial models reached peak DSI levels of 66.8 ± 2.4% vs. 56.2 ± 8.2% (Supplementary Figs. 9, 10). This enhanced performance of the temporal RF models could be attributed to more effective engagement of both excitatory and inhibitory mechanisms of DS, sometimes by combining H&R with B&L organization motifs (Supplementary Figs. 8-10) 14,77,78 or anti-H&R and anti-B&L mechanisms, with the latter often seen as the best performing solution in models unrestrained in all excitatory and inhibitory parameters (DSI = 73.1 ± 2.4%, Fig. 6c).

Robustness of direction selectivity mechanisms across distinct neural architectures.
To evaluate the generalizability of the computational principles identified in DSGC models, I applied the same machine learning-based exploration of DS to models of layer 2/3 cortical pyramidal neurons. The cortical system differs markedly from the retinal circuit in several key ways. First, synaptic input arises from spiking thalamo-cortical and cortico-cortical afferents, unlike the graded synaptic transmission characteristic of retinal circuits. Second, the spatial arrangement of synapses on the postsynaptic neuron does not reflect the retinotopy of their receptive fields 79,80. Third, the asymmetric morphology of pyramidal cells, along with strong electrical filtering in apical dendrites, alters the integration of synaptic inputs and could influence motion computations 62.

Yet despite these structural and biophysical differences, the same DS mechanisms emerged in cortical models (Fig. 7). Together with similar findings in the minimal models (Fig. 3, Supplementary Figs. 3, 5), this result supports the conclusion that computational primitives guiding DS implementation described in this work are robust and largely independent of postsynaptic architecture.

Robustness of direction selectivity mechanisms across variable stimulus presentations
For neuronal computations to be effective, they must remain robust in the face of noise or variability in sensory input. Previous studies have identified both population-level and cell-intrinsic mechanisms that support noise resilience in retinal and cortical circuits involved in direction and orientation selectivity 43,81-86.

To assess the noise tolerance of the different models described above, I evaluated their direction tuning in response to moving bars with variable acceleration. In each simulation run, the bar moved consistently in one direction, but its velocity fluctuated according to the noise level (Fig. 8a). This design ensured that performance reflected the circuit’s ability to detect motion in the preferred direction, independent of confounding effects from motion in the opposite direction, as would occur if the velocity was negative. Thus, any degradation in directional tuning directly indicated reduced resilience to temporal variability within the stimulus.

As expected, increasing velocity variability led to a reduction in the DSI levels across most models (Fig. 8, Supplementary Fig. 12). Notably, two exceptions were seen. First, DS computation implemented with varying sizes of the center RF component performed exceptionally well and maintained high DSI values across all noise levels, in both excitatory and inhibitory setups (Fig. 8b-e). In contrast, models relying on variations in the size of the surround RF component performed poorly; their selectivity deteriorated rapidly with increasing noise, becoming virtually undetectable at higher variability levels (Supplementary Fig. 12).

Methods to identify circuit components participating in DS computations
The conventional approach to identifying presynaptic mechanisms underlying DS typically involves whole-cell recordings or measurements of extracellular spikes or calcium signals. Most studies using this approach have employed stimuli designed to reveal temporal differences in input kinetics or spatial offsets between excitation and inhibition -- strategies motivated by classical models such as the H&R correlator, the energy model, and the B&L framework (Fig. 9a, b) 5-7. Simulations of models trained to implement these computational primitives confirm that these methods are effective at uncovering distinct spatiotemporal activation patterns, namely the presence of oriented space-time RFs, central to DS computation in those scenarios (Fig. 9c).

However, this traditional approach fails to identify DS circuits that rely on alternative mechanisms, such as spatial differences in receptive field structure or center-surround interactions (Fig. 9d-g, Supplementary Fig. 13). As a result, many of the DS-capable circuit configurations uncovered in this study remain inaccessible using conventional recording paradigms. Instead, these circuits require direct examination of the presynaptic population. Experiments targeting this layer -- via electrical recordings, calcium imaging of cell bodies, or measurements of synaptic activity (e.g., pre- or postsynaptic calcium signals or neurotransmitter release) -- could be combined with the same 1D noise paradigm to reveal asymmetries in center-surround RF structure among presynaptic cells (Fig. 9d-g, Supplementary Fig. 13). This strategy, combined with more sophisticated RF mapping from motion responses 49, could be used to uncover the broader range of computational primitives participating in motion detection proposed in this study.

Methods
Multicompartmental simulations were performed using NEURON 8.2 (www.neuron.yale.edu/neuron). DSGC (352 segments) and cortical layer 2/3 (69 segments) morphologies were adapted from 49,83. The initial global passive parameters were as follows: passive conductance = 4e-4 S/cm2, membrane capacitance = 1 µF/cm2, reversal potential = -60 mV and axial resistance = 150 Ωcm. Simulations were executed for 3 – 4 seconds with a time step of 0.1 ms.

The postsynaptic cells received innervation from 100 excitatory and potentially 100 inhibitory inputs, making one synapse each. Presynaptic responses and analysis was done in Igor Pro 9 (www.wavemetrics.com) or Python 3.10. In retinal configuration, presynaptic RFs were distributed randomly in a circle covering 300 µm centered on the soma of the DSGC and connected to the nearest postsynaptic dendrite. Cortical circuits covered a larger area (400 µm) and the connectivity was made to random postsynaptic branches.

To determine presynaptic RF activation, the spatial overlap between two-dimensional Gaussian functions describing the center and surround RF components and the shape of the stimulus (see below) was computed separately for each time step (Δt = 1-10 ms).

Area (µm2) corresponds to the normalized fraction of the stimulated RF area at time t, computed for each component. For full-field static flashes, the entire center and surround were activated when the stimulus was presented. For moving stimuli, area was the pixel-by-pixel sum of the area of the Gaussian function describing the center or surround RF component that spatially overlapped with the stimulus (equation 1):

area = ∑∑stimulus × e{ -[(xθ)/size_x]^2 - [(yθ)/size_y]^2 } / (√ln2 * √ln2)

Where size (µm) is the full width at half maximum (FPWHM) amplitude of the corresponding RF component. xθ and yθ represent centered and potentially rotated x and y pixel coordinates, calculated as (equation 2):

xθ = (x-center_x)*cos(θ)+(y-center_y)*sin(θ)
yθ = -(x-center_x)*sin(θ)+(y-center_y)*cos(θ)

Subsequently, RF center and surround responses at time t were determined independently from the following equations (equation 3):

RF_t = (area_t - RF_t-1) * Δt / rise_time + RF_t-1 * adaptation_t-1
adaptation_t = max (0, adaptation_t-1 - RF_t * Δt / decay_time)

Where adaptation represents depletion of the RRP. In models where RRP dynamics were not simulated, the adaptation term was set to 1 for the entire duration of the simulation. Rise and decay times were bound to be larger than Δt. The widths were positive and were smaller than 200 µm for the center component.

The integration of center and surround RFs was calculated as the difference between the RF components' activation (equation 4):

RF_full,t = ω_syn * max(0, RF_center,t - RF_surround,t * ω_surround)

Where ω_syn is the synaptic weight and ω_surround is the strength of the surround, bounded between 0 and 1. In retinal models, the postsynaptic synaptic conductance was set to RF_full values. Synapses in cortical circuits were activated by simulated presynaptic spikes, whose rates were generated based on RF_full values: at each time point, a spike was generated if a random value drawn from a uniform distribution in the range [0, 0.3] was larger than the RF_full,t.

Short-term plasticity
The amplitude of synaptic release was constant unless noted otherwise.
In simulations with short-term dynamics, short-term depression was implemented using a depression factor syn_dep,amp in the range [0, 1]. With a starting value of 1, representing a full readily releasable pool, synaptic depression was modulated by RF activation at each time point by equation 5:

syn_dep,t = syn_dep,t-1 * syn_fac,amp * (1 - syn_dep,amp * RF_full,t)

Where syn_dep,t is the depression factor at time t, which recovered with a time constant syn_dep,τ (ms) as follows (equation 6):

syn_dep,t = syn_dep,t-1 * (1 - syn_dep,t-1 / syn_fac,τ / Δt)

To calculate the presynaptic release, RF_full,t was multiplied by syn_dep,t.

Short-term facilitation was calculated using the following algorithm: facilitation value syn_fac,t and steady-state facilitation of the releasable pool syn_fac,amp were initially set to zero. At each time point, the value of syn_fac,t was recalculated (equation 7):

syn_fac,t = syn_fac,amp * RF_full,t + syn_fac,t-1

The facilitation value at each time point had a slow decay to zero with a time constant syn_fac,τ (units, ms; equation 7):

syn_fac,t = syn_fac,t - syn_fac,t * Δt / syn_fac,τ

Last, presynaptic release was increased by a factor of 10*syn_fac,t. As syn_fac,t ranged from zero to one; facilitation could increase release amplitude by a factor of ten.

Voltage-gated channels
To enable spiking responses, sodium and potassium conductances were introduced at the postsynaptic soma, following 95, with maximum densities of 10,000 and 2,000 pS/µm², respectively. Spike counts were quantified as the number of voltage crossings at 0 mV.

Training of genetic algorithms
For each generation, 10 NEURON models were executed simultaneously. In most simulations, the stimuli consisted of full-field bars moving at five different speeds: 0.25, 0.5, 1, 2 and 4 mm/s across 12 directions, presented within 800 µm-wide square arena. The height of the bar was 1 mm, its duration was 0.2 seconds. The intensity of the stimuli was set to 1 and the background was set to zero (AU).

In most models, the presynaptic population was divided into 4 groups. A separate RF description was instantiated with random values drawn from a uniform distribution for each of the presynaptic groups. In the synaptic weight model, one group with the same RF was used.

The following limits to the initial values were used:
Response delay (only in models with kinetic changes), 0 – 50 ms. Center width FWHM, 10 – 200 µm. Surround width FWHM > 10 µm. Activation / inactivation / depression / facilitation time constants > 10 ms. Surround strength, 0 – 1. Depression / facilitation factors, 0 – 1.

In the first generation, random values were assigned to seed each network model. Likewise, random initial values were chosen for passive conductance (ranging from 1e-5 to 1e-3 S/cm2), and axial resistance (constrained between 50 and 200 Ωcm).

In each model run, somatic membrane potential was recorded. For each stimulation speed, the DS index was calculated as a vector sum of the peak postsynaptic somatic potentials normalized by the scalar sum of responses (equation 8):

DSI = ∑(R_θ * cos(θ)) / ∑R_θ

Where R_θ is the response magnitude at direction θ, with θ = 0 set as the preferred direction. To reward generation of large postsynaptic responses, the following directional metric was computed (equation 9):

directional metric = DSI * tanh(maxR_θ)

Where maxR_θ was the largest response value across all directions of stimulation.

Subsequently, the models were ranked based on the average directional metric calculated over the five stimulation speeds. The best-performing model seeded the next generation of candidate solutions. Mutations were then applied to 9 models, leaving one intact to prevent search degeneration. For each parameter describing the presynaptic RF and postsynaptic properties, a normal distribution with a mean of 1 and a standard deviation of 10% was used to determine a scaling factor. This scaling factor was multiplied with the parameter value and combined with a random value drawn from a uniform distribution ranging from -0.015 to 0.015. Two mutation types were implemented: (1) Independent parameters, which were allowed to vary between input groups. (2) Shared parameters that mutated only in the first input group and the new values were applied without change to other input groups.

Typically, the models were evolved over 300 generations. To ensure adequate search quality, low-performing models were evolved over 1000 generations. At a minimum, 50 independent runs were performed for each circuit configuration.

Minimal models
Simplified DS models consisted of two presynaptic cells and a postsynaptic detector that summed presynaptic waveforms. Genetic algorithms trained the minimal models based on detector responses, with a reward function similar to the full model. The key difference was that the direction selectivity index was calculated from two directions and for a single velocity.

Analysis
Amplitude-based DS mechanism was determined by analyzing the similarity of motion response waveforms. First, synaptic responses to a bar moving at 1 mm/s were temporally aligned to ensure independence from synaptic position. Next, response overlap was calculated for all time points where at least one of the synaptic waveforms exceeded zero. Overlap was considered absent if less than 75% of the waveforms reached 20% or less of their peak response. Otherwise, the overlap was computed using the average waveform across the synaptic population for that time point. Finally, the mean value of the overlap waveform was used as the metric for the amplitude-based mechanism.

The extent of the alignment-based DS mechanism was computed similarly, but using actual motion responses in the preferred direction and the overlap required 25% of the population waveforms to reach at least 20% or more of their peak responses.

The presence of H&R and anti-H&R was determined as follows: for each synaptic input, the duration of the motion response was calculated as the FWHM amplitude. A linear fit to the durations of synaptic responses determined if the responses were widening or narrowing along the preferred motion axis, corresponding to anti-H&R and H&R mechanisms, respectively.

Finally, the metric for the mechanism was calculated using the narrowest and widest responses as follows (equation 10):

metric_(anti)H&R = (1 - narrowest / widest)^2

The position along the amplitude-alignment and (anti-H&R)-H&R axes was taken as the difference between the corresponding metrics.

To determine the inhibitory mechanism, two time points that divided the area under the curve of the excitatory response in the preferred direction into equal thirds were found. Anti-B&L and B&L mechanisms were computed from the fraction of the inhibitory drive in the same direction within the first and the last phase of the excitatory response, respectively.

The participation of the inhibitory amplitude-based algorithm was computed by equation 11:

metric_amp,inhib = max (0, (amp_Null - amp_Pref) / (amp_Null + amp_Pref))

Where amp_Null and amp_Pref represent the peak amplitude of the inhibitory drive in the null and preferred directions, respectively.

The 'pause' mechanism was computed when both anti-B&L and B&L mechanisms were below 20% (equation 12):

metric_pause = (amp_Null - amp_Pref*) / (amp_Null + amp_Pref*)

Where amp_Null was computed as before and amp_Pref* represents the peak amplitude of the inhibitory drive within the middle time window of the excitatory response.

Position along the cardinal axes was computed similarly to the excitatory mechanisms. A computational primitive was considered to be present if the position distance was more than 20% from the origin along the corresponding axis.

Variable acceleration and 1D noise stimulation
To evaluate the robustness of evolved models under variable stimulus conditions, I introduced motion noise by updating the bar’s velocity every 50 ms. Each new speed was sampled from the positive half of a normal distribution. The width of the distribution was set by the imposed noise level, up to a standard deviation of 5 mm/s. These variable-velocity trajectories were then presented in the preferred and null directions to calculate single-trial DSI values. DSI statistics for each noise level were obtained by averaging across 50–100 trials.

For 1D noise stimulation, I used 20 µm-wide vertical bars that spanned the entire visual field. Bar contrast was randomly assigned as either 0% or 100% with equal probability. Contrast values were updated every 20–200 ms. To analyze receptive field characteristics, I examined responses at either the model detector or presynaptic cells, focusing exclusively on the ON (rising) response phase. Signal fluctuations following bar onset were averaged in 20 ms time bins for each spatial location. These measurements produced space-time response maps (Fig. 9), from which the temporal RF component was derived by analyzing the shape of the response over the central position, and the spatial RF was extracted along the different spatial positions at the time of the maximal response.

Drifting gratings stimulation
In some simulations, the models were trained to optimize responses to drifting sinusoidal gratings rather than moving bars. The gratings, presented at 100% contrast on a gray background, moved across the cell’s receptive field at a speed of 1 mm/s. Five spatial frequencies were tested: 100, 200, 400, 800, and 1600 µm/cycle. DSI was calculated using the steady-state peak-to-peak voltage difference measured over the final 30% of each trial instead of the peak response amplitude.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "In this paper, Poleg-Polsky uses a genetic algorithm to find model parameters for circuits that are direction-selective, with connectivity constraints imposed by measured connectivity and neural morphology. The claim is that some of these motifs represent entirely new algorithms for direction-selectivity, beyond what has canonically been considered in the field: the dynamic differences associated Hassenstein-Reichardt correlators, Barlow-Levick veto mechanisms, and motion energy models. If this claim is correct, then this is an important and interesting paper.\nI found this paper clearly written, and the models all looked like they performed as described. However, I had some fundamental questions about choices made by the author and how they affect the results. In some cases, as I describe below, I have doubts about whether some of these models should be considered direction-selective.\nMajor issues\n1) The author uses maximum responses to assess DSI throughout. However, I’m unaware of any existing or proposed biophysical mechanism that would measure the maximum response of a neuron over time in order to make such a comparison. If the stimulus is linearly filtered, then a nonlinearity is required to generate a direction-selective mean signal. Taking the maximum is one option, but an awkward one to implement. The motion energy model supposes one could square (or rectify and square) a signal, then average it, in order to compare to signals – that is far more plausible. The difference is not just academic; I’m including a short snippet of Matlab code at the end of this review that shows how such a difference in nonlinearity can invert the direction-selectivity of a system. I think this work needs to carefully assess the impact of the chosen nonlinearity on its results.\n2) I found the equations excessively difficult to read, in part because they used words instead of variables, but also because there were unit issues that I could not get past. (Equation and line numbers throughout would be nice to have.) In equation 1, for instance, area_t seems to have a numerator in the exponent in which time has a unitless number (length/(speed*t)) subtracted from it. Units should definitely match here. Similarly RF_t in equation 2 shows that RF_t is equal to a set of terms including one that is RF_{t-1}/(rise time), which again, has units that do not work out. I believe these are intended to be discrete versions of differential equations, but I don’t think they are correctly written, and make it difficult to evaluate the operations being performed.\n3) This may be the most important issue. When I was reading, I spent a lot of time thinking about the example in Figure 1c, which shows how identical temporal processing but increasing receptive field sizes results in directional responses to the moving bar (using the max of the trace). In the discussion, the author says that a minimal model with only linear filtering and the max() nonlinearity reproduces this result and the others in Figure 1, Figure 2, and S2. My intuition for this RF-size change model is that when the bar comes from the large side, its entry into the component RFs is staggered in time, but when it comes from the small size, it hits small and large RF edges simultaneously, creating a higher max value when all the inputs are summed. If we rewrite this is a single spatiotemporal linear summation, it will in its simplest form just have, for instance, a higher weight on one side than the other and bilobed temporal filters throughout (so that it responds primarily on entry, as in the figures in 1c). This proposed filter has the same property described here: larger max responses to a bar moving one way than the other.\na. However, the linear RF for the model with changing receptive field widths is spatiotemporally *separable* by virtue of the identical dynamics at all points in space. This means it can be written as h(x)f(t). Adelson & Bergen require spatiotemporal inseparability for a reason: if the spatiotemporal filter is separable, one can convolve the spatial component first, then the temporal component second, and then show that for drifting sinusoids, the results for one direction vs. the other is just a phase shift in output sinusoids of identical amplitude. That is, it is provably the case that such a linear filter cannot yield directional responses to a sinusoid. This is why there exists such a strong emphasize on oriented space-time filters to generate direction-selective responses.\nb. However, I believe such an LN model does yield the proposed responses to a moving bar. If so, should we think of that as a direction-selective model, if it responds equally to drifting sinusoids in both directions? Drifting sinusoids are a standard tool in visual neuroscience, so why should one prefer a bar to a sinusoid (or vice versa) for evaluating direction-selectivity? How should we think of these direction-selective models if they’re stimulus dependent, even for super simple stimuli? I view this as a fundamental issue for this paper – why optimize using moving bars if the results are not direction-selective to moving sinusoids?\nc. (It may also be that the full model including all the synaptic dynamics and facilitation — not the simple LN model approximation referenced in the discussion — is direction-selective to drifting sinusoids. In this case, it is incumbent on the author to sort out which aspects of the model are required: for instance, it could be the rectification at the model synapse is important. But the claims seem as though they would be quite different from the current ones being made.)\nd. This broad formulation becomes a question for many of the other proposed mechanisms that appear to be able to be approximated by a separable space-time filter.\nMore minor comments\n1) DSI color scale in figure 4 isn’t that useful. Maybe move it to the range involved? 0-0.3?\n2) I think there are two prior papers doing similar work using the fly connectome that might be appropriate to cite here: Mano et al. 2021 and Lappalainen et al. 2024. With reference to bars vs. sinusoids, its useful to note that these both use natural scene inputs to optimize their models.\n3) I had a difficult time understanding the meaning of anti-HR and anti-BL.\n4) Discussion sentence: “For this reason, the models developed here are simpler than classic frameworks such as the full Hassenstein-Reichardt correlator or the energy model (Adelson & Bergen, 1985; Hassenstein, 1956), both of which rely on more elaborate nonlinear computations and specialized spatiotemporal filtering”. I don’t think this sentence is even remotely true. This work relies on an unknown mechanism finding the maximum of a timetrace, relies on quite complicated synaptic facilitation and dynamics, and integrates over entire neurons using the NEURON simulator. The motion energy model relies on linear filtering and squaring, while the HR model relies on linear filtering and multiplying. In fact, all the extra nonlinearities in this paper make it rather hard to assess in some cases what is really causing direction-selectivity in order to compare it with these simple models (see major comment above).\n5) The paper says that the optimized models were “Analyzed systematically”, but this seems to mean “understood by the process of staring”. It’s fine to analyze this way, but is there a way to plot the model parameters or cluster them or show that there’s really only 4 or 6 dimensions?\n6) In the discussion, the author notes he did not observe DS inhibition in these models. But are the models rich enough to generate DS inhibition? That would require a nonlinearity upstream of the inhibitory input to the cell. It was hard to tell whether this was in principle possible. It is also a bit surprising, considering the well-known direction selectivity of the SACs upstream of DS-RGCs that confer some if not most of their DS properties.\n7) Since this an entirely computational project, the code should be made freely available -- I did not see any links to a repository.\nMatlab code snippet:\n% test of when max resp vs. other measure could lead to opposite\n% results for directionality, using a linear model to make the response and different\n% nonlinearities to evaluate strength\nf0 = zeros(1,40);\nf1 = f0; % filter at point 1\nf2 = f0; % filter at point 2\nf1(11:21) = 1; % mostly low\nf1(11) = 2; % high at beginning\nf2(6:16) = 1; % mostly low\nf2(16) = 2; % high at end\nfigure;\nsubplot(2,1,1);\nplot(f1); ylabel('f1');\nsubplot(2,1,2);\nplot(f2); ylabel('f2');\nxlabel('time step');\nr_plus = f1 + circshift(f2,[0 5]); % this simulates two delta function inputs displaced in time and space\nr_minus = circshift(f1,[0 5]) + f2; % and the opposite direction\nfigure;\nsubplot(2,1,1);\nplot(r_plus); ylabel('direction plus');\nsubplot(2,1,2);\nplot(r_minus); ylabel('direction minus');\nxlabel('time step');\ndisp('max of the two:');\ndisp(['r_plus,r_minus = ' num2str(max(r_plus)) ', ' num2str(max(r_minus))]);\ndisp('sum square of two:');\ndisp(['r_plus,r_minus = ' num2str(sum(r_plus.^2)) ', ' num2str(sum(r_minus.^2))]);", "reviewer_2": "In this manuscript, the author uses a machine learning approach to uncover algorithmic building blocks for direction selectivity in models of visual neurons. They find that, when optimized for direction selectivity, these models tend to fall into one of several computational “primitives,” some of which align with existing models of direction selectivity like the Hassenstein-Reichardt correlator or the Barlow/Levick model, which gives some credibility to the overall approach. Interestingly, this genetic algorithm approach also reveals some new mechanisms for direction selectivity that rely on things like differences in upstream receptive field sizes and orientations. Another key finding here is that downstream/postsynaptic integration schemes don’t much alter the basic upstream building blocks that the genetic algorithm learns, indicating that these mechanisms might be shared across different circuits like retina, cortex, and perhaps other visual systems. Overall, this is an interesting manuscript that has important implications for our understanding of direction selectivity. The genetic algorithm approach is very interesting and somewhat novel (or at least unusual), and it seems well suited for this question, although I do wish there were more details about the evolution of the models (see below). I have two primary concerns about the manuscript that I think need to be addressed. First, the details of the genetic algorithm and the evolution of the models was sparse, which made it hard to evaluate the results, in particular the claim about there being distinct computational “primitives.” Second, the narrow range of stimuli used makes it hard to evaluate how general these mechanisms are.\nMajor issues:\n1. How are solution families / computational primitives defined?\nThroughout the paper, it is not clear how distinct the model solutions are, and how often each type of solution was arrived at. How did you draw boundaries between solution families? The major concern here is that when looking at this huge diversity of models that emerge, you are imposing structure based on the expectations of what the primitives should look like, rather than the structure really being present.\n- Can you show clustering for more of these solution families, as you do in Fig. 3?\n- How did you classify a model in Fig. 4 as a “pure” model or a “hybrid” model if presumably a single simulation results in a model that belongs to one or another cluster?\n- What fraction of the evolved networks came to each solution class? Does changing the hyperparameters of the genetic algorithm (mutation rate, etc) change this?\n2. Stimulus dependence\nIt seems like all of these simulations were run with moving bars oriented perpendicular to the axis of motion. How dependent are the primitives on the identity of the stimulus? What if you used a thicker/large moving feature? Or even widefield motion? I don’t think you need to re-run every simulation but some exploration of different stimuli would bolster the claim that the model is learning computational primitives and not just strange solutions that only work for one stimulus (the gap-in-inhibition solution in Fig. 5 comes to mind).\n3. Biological feasibility\nThe paper claims that these computational primitives are biologically feasible in pretty strong terms. E.g. the abstract “All mechanisms are biologically plausible and correspond to known physiological and anatomical motifs.” I suppose plausibility is subjective but throughout the results there are models that have a distinctly non-biological feel to them, like things with non-monotonic changes in RF size or orientation across the NP axis (Fig. 2, 5), or a single presynaptic cell with a huge RF (Fig. 6). Can this be addressed or softened?\n4. Null / preferred nomenclature\nThe labeling of the “null side” and “preferred side” seems backwards in places. The text says the null side is the “position from which the null direction stimulus enters the RF” which is my understanding of the existing convention. But Fig. 1a has the null direction stimulus entering on what is labeled the “pref side.” And again in the description of Fig. 5 the text reads “… stronger inhibitory input on the preferred side …” but again isn’t that the null side? That would agree with the general idea that the BL model relies on “null side” inhibition – a description included in, for example, the Fried 2002 paper cited there and I’m sure other papers as well. Can you clear up this discrepancy since it really confused me while reading the paper.\nMinor\n- Can you add a legend item for the dotted line in 1d (top)\n- I’m confused about the description of the model in Fig. 3 – text says “band-like organization of stronger synaptic inputs along the motion preference axis” but isn’t the band perpendicular to the NP axis? Panel b looks like there are stronger inputs on both the N and P side of the RF – how is this similar to the 2 input model and how does it generate DS responses?\n- Call out to Fig. 1g – there is no Fig. 1g. I assume this should only refer to Supp Fig. 1?", "reviewer_3": "This is an important and timely manuscript that systematizes different computational approaches to motion selectivity. The manuscript is well written and technically correct. It was interesting to read how motion selectivity can be achieved not only through spatial differences in temporal filtering as int he HR detector but also through spatial differences in receptive field size, ratios, and surround properties. One aspect that I can suggest for improvement is that it was not clear how well clustering into different primitives worked, whether clusters were well separated or it was just a continuum between different axes (temporal kinetics, amplitude, receptive field size and asymmetry). In other words, if there is strong difference in orientation differences does it preclude observing differences in dynamics. Or to pose the question in yet another way: Figure 4B shows the distribution of different possible effects individually but what about correlations between them (combinations of primitives).", "reviewer_4": "I appreciate the effort the author has made in exploring computational principles underlying elementary motion detection.\nThe manuscript is ambitious in scope, presenting a wide range of simulated circuit architectures and proposing that these converge on a small number of “computational primitives” for motion detection. The use of minimal, tractable feedforward models as a way to isolate and examine core features of direction selectivity is a commendable strategy, and I agree with your stated rationale that such simplification can reveal useful insights.\nHowever, despite the technical breadth of the work, I find that the manuscript does not convincingly achieve its conceptual goals, nor does it provide sufficient clarity about the key questions it is addressing. My concerns center around the overall framing and purpose of the study, rather than any specific technical errors in the simulations. Below, I outline several major issues that I believe need to be addressed for the study to be considered suitable for publication.\n1. Lack of a Clear Scientific Question or Rationale\nThe manuscript surveys a large number of possible circuit configurations that can implement direction selectivity (DS), but it remains unclear what the main scientific question is. Are you testing specific hypotheses about biological mechanisms? Are you proposing new models for experimental validation? Or are you primarily interested in mapping the design space of motion-detecting circuits? This ambiguity makes it difficult to assess the contribution of the work.\nYou state in the introduction that the study aims to explore how basic receptive field features and spatial arrangements give rise to directional tuning. However, the simulations appear largely exploratory and descriptive, lacking a well-defined problem or hypothesis. Without a clear rationale, it is difficult to judge the value or relevance of the results beyond their existence.\n2. The Study Feels Preliminary Rather Than Conclusive\nAlthough the manuscript presents numerous models and simulation results, the overall impression is one of an exploratory analysis or early-stage investigation. Many of the findings—such as the rediscovery of Hassenstein-Reichardt (H&R) and Barlow-Levick (B&L)–like computations—are already well known and expected. The identification of new combinations or variants of these does not in itself constitute a major conceptual advance unless their biological plausibility or explanatory power is convincingly demonstrated.\nThe models are presented in a systematic but somewhat catalog-like fashion, and the paper often reads like a compilation of simulation experiments rather than a focused scientific narrative. This is not inherently a flaw, but I believe that, for a journal like Nature Communications, a stronger framing and clearer implications are expected.\n3. Ambiguity of Machine Learning’s Role\nYou describe using machine learning to explore a \"vast combinatorial space\" of circuit solutions. However, the exact nature of the machine learning method, its justification, and its biological relevance are not clearly explained. While evolutionary algorithms can certainly be used to optimize network parameters, the study does not make clear what has been learned from this process, other than the fact that many solutions exist.\nMoreover, the manner in which models are selected and evaluated—based on performance metrics like directional selectivity index (DSI)—follows an ML benchmarking logic that is not necessarily meaningful in the context of biological interpretation. There is little discussion of how these models relate to known circuitry or why certain parameter regimes might be favored by real neural systems. This weakens the biological impact of the findings.\n4. Biological Plausibility Is Not Adequately Addressed\nAlthough you note that all mechanisms are “biologically plausible,” there is little evidence or discussion to support this claim. Many of the circuit motifs presented could just as easily be dismissed as mathematical possibilities with no relevance to actual neural systems. Without experimental data or stronger theoretical justification, it is impossible to tell which of these mechanisms the brain might actually use.\nFurthermore, because the design space is enormous, the models explored here necessarily occupy a minuscule region of the total parameter space. This raises concerns about the representativeness of the results. If no constraints (biological, theoretical, or functional) are applied to restrict the search space, then the significance of any particular solution is unclear.\n5. Lack of Hypothesis-Driven Structure\nA key issue throughout the manuscript is the absence of hypotheses that can be tested or evaluated. This is particularly problematic given that many of the findings are not surprising: for instance, that asymmetric RFs or synaptic weights can produce DS. The absence of falsifiability or predictive power makes the work more descriptive than explanatory.\nAs you note in the discussion:\n\"I found a strikingly large number of distinct functional organizations of a simple feedforward network capable of extracting motion direction from sensory input.\"\nThis sentence encapsulates the core problem: while the variety of possible implementations is acknowledged, the manuscript does not provide a framework for determining which of them are meaningful. The absence of such a framework reduces the scientific value of the result.\n6. Presentation Style and Scope\nThe extensive number of figures (1–10), while individually clear, contributes to a sense of fragmentation. There is a lack of narrative cohesion tying the results together or building toward a central conclusion. The manuscript feels more like an organized summary of simulations than a focused argument.\nThe use of phrases like “previously unrecognized circuit motifs” or “new computational primitives” may overstate the novelty of what are, in essence, variations on known mechanisms. Without empirical grounding or theoretical analysis, the use of such terms may risk being misleading.\nWhile the technical execution of the simulations is sound, and the general idea of using machine learning to explore motion detection mechanisms is potentially valuable, the current manuscript lacks the conceptual clarity and scientific focus needed to support its claims. I encourage the author to reconsider the framing of the work, tighten the narrative around specific questions, and critically assess the biological meaning of the results.\nI hope this feedback is helpful in guiding the next steps for this promising but currently preliminary line of research.\nMinor Comments\nThe manuscript lacks page and line numbers, which makes it difficult to reference or discuss specific passages in detail. Including these elements is important for facilitating constructive feedback and precise editorial communication.\nThere is inconsistency in figure panel labeling. For example, panels are labeled as \"Ci, Cii\" in the figures but referred to simply as \"1C\" in the main text. This unconventional and inconsistent notation can be confusing for readers and reviewers alike. Please consider standardizing the figure references to align with common scientific writing conventions."}, "author_response": "I appreciate the constructive feedback provided by the reviewers. Below you will find my\nresponses to your concerns and suggestions marked in blue in this document and in the\nmanuscript text.\nI apologize for not including the link to simulation code in the main text, it appeared in the\nsubmission metadata: https://github.com/PolegPolskyLab/DS-mechanisms. The information is\nnow included in the manuscript under ‘code avaliability’\nReviewer #1 (Remarks to the Author):\nIn this paper, Poleg-Polsky uses a genetic algorithm to find model parameters for circuits that\nare direction-selective, with connectivity constraints imposed by measured connectivity and\nneural morphology. The claim is that some of these motifs represent entirely new algorithms for\ndirection-selectivity, beyond what has canonically been considered in the field: the\ndynamic differences associated Hassenstein-Reichardt correlators, Barlow-Levick veto mechanisms, and\nmotion energy models. If this claim is correct, then this is an important and interesting paper.\nI found this paper clearly written, and the models all looked like they performed as described.\nHowever, I had some fundamental questions about choices made by the author and how they\naffect the results. In some cases, as I describe below, I have doubts about whether some of\nthese models should be considered direction-selective.\nMajor issues\n1) The author uses maximum responses to assess DSI throughout. However, I'm unaware of any\nexisting or proposed biophysical mechanism that would measure the maximum response of a\nneuron over time in order to make such a comparison. If the stimulus is linearly filtered, then a\nnonlinearity is required to generate a direction-selective mean signal. Taking the maximum is\none option, but an awkward one to implement. The motion energy model supposes one could\nsquare (or rectify and square) a signal, then average it, in order to compare to signals – that is\nfar more plausible. The difference is not just academic; I'm including a short snippet of Matlab\ncode at the end of this review that shows how such a difference in nonlinearity can invert the\ndirection-selectivity of a system. I think this work needs to carefully assess the impact of the\nchosen nonlinearity on its results.\nThe reviewer is correct that peak membrane potential differences between trials do not always\ntranslate directly into proportional changes in spike output, as active properties in some\nneurons can rarely lead to complex nonlinear behaviors (see for example (Wienbar & Schwartz,\n2022)). However, in both cortical and retinal direction-selective cells, stronger depolarization in\nthe preferred direction is the norm, and this makes the voltage-based analysis biologically\nrelevant.\nThe layout of the simulation closely follows classical experimental design and well-established\nanalysis techniques, where the extent of direction selectivity is measured from responses to\ndifferent trials in which stimuli move in different directions and the peak response or number of\nspikes are used to calculate the vector sum of the response, as I and others have done in the\nDSGCs. The following representative citations (too many to cite all!) from the cortical literature\n(Freeman, 2021; Lien & Scanziani, 2018; Wilson et al., 2018) and the fly DS system (Badwan et\nal., 2019; Gruntman et al., 2021; Maisak et al., 2013) show similar, peak-based analysis.\nPlease also note that not all DS cells spike, with one well-known example being starburst\namacrine cells that generate directionally tuned calcium signals in their dendrites, but the\nanalysis is quite similar (Euler et al., 2002; Morrie & Feller, 2018; Poleg-Polsky et al., 2018).\nTo further confirm the correspondence between subthreshold voltage signals and spiking\noutput, I have now included simulations in which models were first trained on membrane\npotentials of the detector (the analog nature of the 'subthreshold' voltage signal has smooth\nderivatives, greatly simplifying model training) and then equipped post hoc with somatic\nvoltage-gated sodium and potassium channels. As shown in the new Supplementary Figure 1,\nthis approach yields robust, directionally tuned spiking responses that closely mirror the\nunderlying voltage selectivity.\n2) I found the equations excessively difficult to read, in part because they used words instead of\nvariables, but also because there were unit issues that I could not get past. (Equation and line\nnumbers throughout would be nice to have.) In equation 1, for instance, area_t seems to have a\nnumerator in the exponent in which time has a unitless number (length/(speed*t)) subtracted\nfrom it. Units should definitely match here. Similarly RF_t in equation 2 shows that RF_t is equal\nto a set of terms including one that is RF_(Badwan et al., 2019)/(rise time), which again, has\nunits that do not work out. I believe these are intended to be discrete versions of differential\nequations, but I don't think they are correctly written, and make it difficult to evaluate the\noperations being performed.\nThank you for noticing. Units were added to all equations. I broke equation 1 into two\ncomponents and now describe how rotation of elongated RFs was implemented. Old equation 2\nis now fixed (it was missing a time step value)\n3) This may be the most important issue. When I was reading, I spent a lot of time thinking\nabout the example in Figure 1c, which shows how identical temporal processing but increasing\nreceptive field sizes results in directional responses to the moving bar (using the max of the\ntrace). In the discussion, the author says that a minimal model with only linear filtering and\nthe max() nonlinearity reproduces this result and the others in Figure 1, Figure 2, and S2. My\nintuition for this RF-size change model is that when the bar comes from the large side, its entry\ninto the component RFs is staggered in time, but when it comes from the small size, it hits small\nand large RF edges simultaneously, creating a higher max value when all the inputs are\nsummed. If we rewrite this is a single spatiotemporal linear summation, it will in its simplest\nform just have, for instance, a higher weight on one side than the other and bilobed temporal\nfilters throughout (so that it responds primarily on entry, as in the figures in 1c). This proposed\nfilter has the same property described here: larger max responses to a bar moving one way than\nthe other.\na. However, the linear RF for the model with changing receptive field widths is spatiotemporally\n*separable* by virtue of the identical dynamics at all points in space. This means it can be\nwritten as h(x)f(t). Adelson & Bergen require spatiotemporal inseparability for a reason: if the\nspatiotemporal filter is separable, one can convolve the spatial component first, then the\ntemporal component second, and then show that for drifting sinusoids, the results for one\ndirection vs. the other is just a phase shift in output sinusoids of identical amplitude. That is, it is\nprovably the case that such a linear filter cannot yield directional responses to a sinusoid. This is\nwhy there exists such a strong emphasize on oriented space-time filters to generate direction-\nselective responses.\nI completely agree with the reviewer that exploring whether proposed mechanisms can respond\nto stimuli beyond a moving bar is important. As the new Supplementary figure 2 (recreated\nhere) shows, both simple 2-input and the full models can create strong direction-selective\nresponses to drifting sinusoidal gratings. This is true for models that vary in the spatial or\ntemporal characteristics, and for models in which different manipulations of the surround were\nallowed, similar to figure 2. The revised manuscript now shows that all examined RF\ncomponents are capable of responding preferentially to a specific direction of a drifting grating\nwhile utilizing the algorithmic solutions found with moving bar stimulation (new Supplementary\nfigures 3 and 10).\nSupplementary Figure 2. Responses of kinetic and spatial solutions to drifting grating\nstimulation. a, Top, space-time plot showing pixel intensity of a drifting grating stimulus\nactivating two inputs positioned 200 µm apart. The model was trained to produce stronger\nactivation (measured as the peak-to-trough amplitude during the steady-state phase) for\nmotion from bottom to top. Schematic on left shows the temporal and spatial RF properties of\nthe evolved model; horizontal line widths indicate RF size. Middle, stimulus intensity at the RF\ncenters (gray) and response kinetics of the two cells (green). Bottom, linear summation of the\ntwo detectors (black). Vertical bars indicate response amplitude. Note the stronger signal in the\npreferred direction. Gray, linear combination of the stimulus intensity traces shows similar\namplitude. b, Responses to five spatial frequencies in the full model with variable RF kinetics. c-\nd, As in a-b, but for models with variable RF diameters. In both cases, circuit architecture of the\noptimal solutions trained on drifting gratings were comparable to those trained on moving bars.\nReviewer's intuition is almost correct regarding the solution the models found for the spatially-\ndiverse RFs. I the attached reviewer's figure 1, I illustrated the principle of such computation.\nFigure R1. Toy examples of direction selectivity computation in models with spatially distinct\nreceptive fields.\na, Schematic of two spatially offset cells with different RF sizes (red cell: large RF, blue cell: small\nRF; RF outlines shown as dotted curves). Soma positions are indicated by small circles. The RFs\nare positioned so that their edges meet at a single point. A narrow bar stimulus moves across the\nRFs, starting either from the side farthest from the meeting point (null direction, left, stimulus\nmotion is from left to right) or from the opposite side (preferred direction, right, stimulus motion\nis from right to left). The spatial arrangement is depicted from the stimulus's point of view, so\ntemporal progression is left-to-right. Inset: same arrangement shown from an external observer's\nperspective, showing the actual direction of the stimulus (from right to left).\nb, Model responses when each cell encodes stimulus appearance within its RF. In this toy\nexample, the cell responds only when the stimulus reaches the cell. In the null direction,\nresponses are temporally staggered, as the bar first enters the red cell's RF; in the preferred\ndirection, responses are synchronous, enabling stronger summation. Vertical dotted lines mark\nstimulus entry into each RF.\nc, DS computation in a scenario more representative of the manuscript's models, where each\ncell's input is proportional to the degree of overlap between the stimulus and its RF (dotted\ncurves), followed by temporal filtering in both cells (solid curves). The temporal filter is identical\nin both cases, but because the overlap of the stimulus with the large RF sustains is longer, its\nfiltered response is also prolonged.\nd, Summed signal in a downstream linear detector depicting stronger activation in the preferred\ndirection.\nIn principle, a model that responds to the entrance of the stimulus to the RF, can lead to DS\ndetection that is speed invariant; in the arrangement shown in Fig. R1 a, the preferred-direction\nstimulus is always going to enter the RFs of both cells at the same time, drive responses and\nstrong summation. Such a solution can be implemented with spiking cells, but I doubt if it would\nbe able to operate over different contrasts, and I expect noise tolerance to be quite low.\nA more realistic signal processing is shown in Fig. R1 b-c, where the cells' activation at any given\ntime is a function of the overlap between the stimulus and the RF. Because this activation is fully\nsymmetric for both directions (Fig. R1c, the arrangement of the dotted curves is a mirror version\nof each other), a temporal filter is required (Fig. R1c, solid curves). In this example, the filter is\nimplemented with an adaptive synapse that mimics the depletion of the readily releasable pool\n(RRP). This two-step arrangement is conceptually similar to processing in bipolar cells, which\nsample (mostly) independent photoreceptor drives at their dendrites, which are summed and\ntransferred as a membrane potential signal to the axon, where voltage-gated channels and\nsynaptic dynamics perform the temporal filtering step.\nBelow is the code for recreating the model in Fig. R1 c-d (Igor syntax):\n// preferred direction (x, y, last dim represents the two cells):\n•make/o/n=(400,400,2) m_RF= (p-200+60*r)^2+(q-200)^2<(100-60*r)^2\n// null direction\n•make/o/n=(400,400,2) m_RF= (p-200-60*r)^2+(q-200)^2<(100-60*r)^2\n// shared code\n//convolve the stimulus with the RF (assume width=1 point)\n•matrixop/o w_sum=sumrows(m_RF)\n//normalize the stimulus (not strictly required for the model to work)\n•matrixop/o w_max=maxcols(w_sum)\n•w_sum/=w_max[0][0][r]\n//this wave will hold the RRP availability\n•make/o/n=(400,2) w_RF_adaptation=0\n//RRP is decreased as a function of RF activation\n•w_RF_adaptation[1,]=w_RF_adaptation[p-1][q]+w_sum[p][0][q]\n//response profile, computed from the spatial RF activation and available RRP\n•make/o/n=(400,2) w_RF_time=w_sum[p][0][q]*(1-w_RF_adaptation[p][q]/w_RF_adaptation[399][q])\n//combination of the two cells' responses in downstream detector\n•make/o/n=(400) w_RF_full= w_RF_time[p][0]+w_RF_time[p][1]\nI took this detour to illuminate the difference between staggered responses and the solutions\nachieved with the genetic algorithm and presented in Figure 1. Importantly, the temporal\nenvelopes of the signals that are transmitted to the detector during object motion are different\nbetween the presynaptic populations. Thus, although presynaptic cells in the spatial models\n(varying RF size and varying RF orientation) have similarly shaped responses to static stimuli,\nthey are able to impose differential temporal filtering on moving stimuli, including drifting\ngratings. Thus, while the reviewer is correct that purely linear models would fail to respond to\ndrifting gratings, the spatial models presented in the manuscript have sources of nonlinearity\nthat allow them to have DS performance that matches the DS levels seen in kinetic models.\nThe reviewer raises an important difference between the kinetic solutions (Hassenstein &\nReichardt and Adelson & Bergen models) and the spatial solutions discussed above. Kinetic\nsolutions require oriented space-time filters; spatial solutions do not. I have mentioned this\npoint in the manuscript in the results describing DS computations by the surround (Fig. 2d) and\nthis was the motivation to show how these circuits can be identified experimentally (Fig. 10).\nOverall, this observation strengthens the claim that these are different solution families and\nsome intuitions/requirements embedded in one solution family don't apply to the other.\nb. However, I believe such an LN model does yield the proposed responses to a moving bar. If\nso, should we think of that as a direction-selective model, if it responds equally to drifting\nsinusoids in both directions? Drifting sinusoids are a standard tool in visual neuroscience, so\nwhy should one prefer a bar to a sinusoid (or vice versa) for evaluating direction-selectivity?\nHow should we think of these direction-selective models if they're stimulus dependent, even for\nsuper simple stimuli? I view this as a fundamental issue for this paper – why optimize using\nmoving bars if the results are not direction-selective to moving sinusoids?\nc. (It may also be that the full model including all the synaptic dynamics and facilitation — not\nthe simple LN model approximation referenced in the discussion — is direction-selective to\ndrifting sinusoids. In this case, it is incumbent on the author to sort out which aspects of the\nmodel are required: for instance, it could be the rectification at the model synapse is important.\nBut the claims seem as though they would be quite different from the current ones being\nmade.)\nd. This broad formulation becomes a question for many of the other proposed mechanisms that\nappear to be able to be approximated by a separable space-time filter.\nAgain, thank you for the suggestion. As I now show, both the simple and the full models are\ndirectionally selective to drifting gratings. This provides further supporting evidence for the\nclaim that the different new solutions discovered here (spatial, surround contribution, etc) are\ncomparable in their ability to provide DS information to established kinetic solutions.\nMore minor comments\n1) DSI color scale in figure 4 isn't that useful. Maybe move it to the range involved? 0-0.3?\nFixed\n2) I think there are two prior papers doing similar work using the fly connectome that might be\nappropriate to cite here: Mano et al. 2021 and Lappalainen et al. 2024. With reference to bars\nvs. sinusoids, its useful to note that these both use natural scene inputs to optimize their\nmodels.\nThank you for the suggestion, references added!\n3) I had a difficult time understanding the meaning of anti-HR and anti-BL.\nMy goal in this work is to analyze DS circuits at both the algorithmic and implementation levels.\nThe literature contains many examples of different circuit implementations of established\nmodels, such as the H&R correlator. It is well known that multiple circuit configurations can\nrealize the same computation. One of the questions I addressed was whether fundamentally\nnew algorithmic solutions to the DS problem exist.\nTo explore this, it is first necessary to clearly define what constitutes a particular algorithm and\nwhat lies outside of it. In the H&R model, the core principle is the presence of distinct temporal\nfilters arranged in a specific configuration: the slower (delayed) filter faces the preferred\ndirection so that it is driven first. The same logic applies when detecting more complex stimuli,\nsuch as drifting gratings. Reversing the filter order in the H&R algorithm reverses the preferred\ndirection.\nHowever, in many spatially based solutions, the opposite occurs: the positions of the faster and\nslower filters are swapped, yet the detector still responds to the original preferred direction (Fig.\n1). I refer to this computation as \"anti-H&R.\"\nSimilarly, the essence of the Barlow–Levick algorithm is delayed inhibition in the preferred\ndirection. In the \"anti-B&L\" solution, inhibition precedes excitation. In both the B&L and anti-\nB&L schemes, excitation and inhibition overlap in the null direction.\n4) Discussion sentence: \"For this reason, the models developed here are simpler than classic\nframeworks such as the full Hassenstein-Reichardt correlator or the energy model (Adelson &\nBergen, 1985; Hassenstein, 1956), both of which rely on more elaborate nonlinear\ncomputations and specialized spatiotemporal filtering\". I don't think this sentence is even\nremotely true. This work relies on an unknown mechanism finding the maximum of a timetrace,\nrelies on quite complicated synaptic facilitation and dynamics, and integrates over entire\nneurons using the NEURON simulator. The motion energy model relies on linear filtering and\nsquaring, while the HR model relies on linear filtering and multiplying. In fact, all the extra\nnonlinearities in this paper make it rather hard to assess in some cases what is really causing\ndirection-selectivity in order to compare it with these simple models (see major comment\nabove).\nThis is why I have included the simple 2 cell model with linear summation! All the key results\nfrom the full NEURON simulation were reproduced with just two inputs performing relatively\nsimple, nearly linear, processing.\nBy comparison, the full Hassenstein–Reichardt model and the Adelson–Bergen (A&B) model are\nconsiderably more complex than many of the excitatory models presented in this manuscript.\nThe full H&R model involves interactions between two DS detectors, while the A&B model\nincorporates separate ON and OFF channels.\n5) The paper says that the optimized models were \"Analyzed systematically\", but this seems to\nmean \"understood by the process of staring\". It's fine to analyze this way, but is there a way to\nplot the model parameters or cluster them or show that there's really only 4 or 6 dimensions?\nThis is only partially true. I am a strong proponent of \"staring\" at the data (by which I mean\nexploring multiple types of visualizations to gain different perspectives on the underlying\ncomponents). To this end, I dissected all models in a way that displays all major parameters.\nOne such plot is reproduced below.\nFigure R2. Illustration of components of the 'flexible surround size' model family.\nEach row represents an optimal solution from the trained models. Columns, from left to right:\nMean DSI measured across 5 stimulation velocities (tuning was assessed in 12 directions). Radar\nplot depicting the computational primitives recruited by each solution, color-coded by DSI.\nExcitatory conductance in the null and preferred directions, with strength represented by color\nintensity. The vertical axis represents spatial position over time. ω is the relative strength of the\nsurround for the four input populations. Activation and inactivation time constants are shown\non a logarithmic scale. \"Width\" refers to the half-width at half-max amplitude of the\ncorresponding RF component. The dendrogram shows the distance between solutions,\ncalculated using hierarchical clustering. The models highlighted in Fig. 2c are indicated by black\nsymbols. Note that the only parameter that varies between input populations is the width of\nthe surround; other parameters are shared, but, of course, differ between specific\nimplementations.\nDuring this initial period of detailed inspection, I began to identify recurring computational\nmotifs. In some solutions, input kinetics aligned with H&R predictions. Others showed similarly\nshaped response profiles, a configuration distinct from the H&R scheme, but these profiles were\ntime-shifted to occur simultaneously in the preferred direction. As I show in Fig. 2d, this shift\nwas absent for stationary stimuli. Still other solutions displayed different temporal shapes\naltogether, yet in positions inconsistent with predictions from kinetic models as discussed\nabove.\nAfter cataloging these patterns, I derived equations to distinguish between the different\nalgorithms (see Methods for full details) and assigned names to the resulting computational\nprimitives.\n6) In the discussion, the author notes he did not observe DS inhibition in these models. But are\nthe models rich enough to generate DS inhibition? That would require a nonlinearity upstream\nof the inhibitory input to the cell. It was hard to tell whether this was in principle possible. It is\nalso a bit surprising, considering the well-known direction selectivity of the SACs upstream of\nDS-RGCs that confer some if not most of their DS properties.\nAs noted in the Discussion, the answer is both \"yes\" and \"no.\" In our simulations, the\npresynaptic layer was constrained to have no DS responses. This was intentional: to (1) limit the\nparameter space (see discussion below), (2) adhere to the definition of an elementary motion-\nprocessing circuit with non-DS inputs, and (3) prevent the model from exploiting a trivial\n\"shortcut\" solution.\nNevertheless, the models could generate inhibitory DS responses in the postsynaptic cell using\nthe same mechanisms that produced excitatory DS responses. However, even when inhibitory\nkinetics were allowed to vary, the models did not converge on H&R–like temporal profiles for\ninhibition. The reason was simple: other strategies consistently achieved stronger DS. This result\nchanged when the models were optimized on the drifting grating task. As I show in the new\nSupplementary Figure 10, some models evolved to produce stronger inhibition in the null\ndirection. The Discussion now contains text that compares directionally tuned inhibition in bar\nand drifting grating stimulation .\n7) Since this an entirely computational project, the code should be made freely available -- I did\nnot see any links to a repository.\nMatlab code snippet:\n% test of when max resp vs. other measure could lead to opposite\n% results for directionality, using a linear model to make the response and different\n% nonlinearities to evaluate strength\nf0 = zeros(1,40);\nf1 = f0; % filter at point 1\nf2 = f0; % filter at point 2\nf1(11:21) = 1; % mostly low\nf1(11) = 2; % high at beginning\nf2(6:16) = 1; % mostly low\nf2(16) = 2; % high at end\nfigure;\nsubplot(2,1,1);\nplot(f1); ylabel('f1');\nsubplot(2,1,2);\nplot(f2); ylabel('f2');\nxlabel('time step');\nr_plus = f1 + circshift(f2,[0 5]); % this simulates two delta function inputs displaced in time and\nspace\nr_minus = circshift(f1,[0 5]) + f2; % and the opposite direction\nfigure;\nsubplot(2,1,1);\nplot(r_plus); ylabel('direction plus');\nsubplot(2,1,2);\nplot(r_minus); ylabel('direction minus');\nxlabel('time step');\ndisp('max of the two:');\ndisp(['r_plus,r_minus = ' num2str(max(r_plus)) ', ' num2str(max(r_minus))]);\ndisp('sum square of two:');\ndisp(['r_plus,r_minus = ' num2str(sum(r_plus.^2)) ', ' num2str(sum(r_minus.^2))]);\nReviewer #1 (Remarks on code availability):\nI did not find any code referenced in the methods or available with the manuscript.\nReviewer #2 (Remarks to the Author):\nIn this manuscript, the author uses a machine learning approach to uncover algorithmic building\nblocks for direction selectivity in models of visual neurons. They find that, when optimized for\ndirection selectivity, these models tend to fall into one of several computational \"primitives,\"\nsome of which align with existing models of direction selectivity like the Hassenstein-Reichardt\ncorrelator or the Barlow/Levick model, which gives some credibility to the overall approach.\nInterestingly, this genetic algorithm approach also reveals some new mechanisms for direction\nselectivity that rely on things like differences in upstream receptive field sizes and orientations.\nAnother key finding here is that downstream/postsynaptic integration schemes don't much alter\nthe basic upstream building blocks that the genetic algorithm learns, indicating that these\nmechanisms might be shared across different circuits like retina, cortex, and perhaps other\nvisual systems. Overall, this is an interesting manuscript that has important implications for our\nunderstanding of direction selectivity. The genetic algorithm approach is very interesting and\nsomewhat novel (or at least unusual), and it seems well suited for this question, although I do\nwish there were more details about the evolution of the models (see below). I have two primary\nconcerns about the manuscript that I think need to be addressed. First, the details of the\ngenetic algorithm and the evolution of the models was sparse, which made it hard to evaluate\nthe results, in particular the claim about there being distinct computational \"primitives.\"\nSecond, the narrow range of stimuli used makes it hard to evaluate how general these\nmechanisms are.\nMajor issues:\n1. How are solution families / computational primitives defined?\nThroughout the paper, it is not clear how distinct the model solutions are, and how often each\ntype of solution was arrived at. How did you draw boundaries between solution families? The\nmajor concern here is that when looking at this huge diversity of models that emerge, you are\nimposing structure based on the expectations of what the primitives should look like, rather\nthan the structure really being present.\n- Can you show clustering for more of these solution families, as you do in Fig. 3?\n- How did you classify a model in Fig. 4 as a \"pure\" model or a \"hybrid\" model if presumably a\nsingle simulation results in a model that belongs to one or another cluster?\n- What fraction of the evolved networks came to each solution class? Does changing the\nhyperparameters of the genetic algorithm (mutation rate, etc) change this?\nPrompted by this comment, I now present the distribution of all individual solutions within the\ncomputational primitives space. This addition substantially increases the information available\nto the reader, revealing both the frequency of each primitive and their combinations.\nThe computational primitives form a continuum, quantified by equations provided in the\nMethods. A primitive was considered present if its positional displacement along the\ncorresponding axis exceeded 20% of the origin. These quantitative definitions formalize the\nintuitions I developed by \"staring\" at the data (see response to Reviewer 1, comment #5).\nImportantly, every model exhibited at least one primitive (if this were not the case, the list of\nprimitives would be incomplete). While additional algorithmic solutions to motion detection in\nthese architectures may exist, they do not form stable local minima in the solution space.\nIndeed, some primitives emerged only under highly constrained parameter sets -- for example,\nthe all-excitatory amplitude solution achieved only a modest DSI (Fig. 3). Models with flexible\nspatial or temporal RF properties generally did not optimize toward this primitive, as superior\nperformance was attainable with H&R, anti-H&R, or alignment-based computations (Fig. 4).\nSome models expressed more than one primitive. For instance, an all-excitation solution might\nshow waveform width differences along the motion axis, but smaller than those in other\nimplementations. According to my metrics, such a case would plot off the anti-H&R and H&R\naxes, leaning toward the \"Alignment\" mechanism pole (Fig. 4).\nFor clarity and more rigorous quantification, I divided the primitive space into eight octants:\nsolutions in cardinal sectors were labeled \"pure,\" while those on the diagonals were classified as\n\"mixed.\"\nAs suggested by the reviewer, I now also show the clustering of all solution families in the\ncorresponding supplementary figures.\nOverall, the genetic-algorithm hyperparameters had little effect on the final outcome. As one\nwould expect, a very low mutation rate could prevent convergence within a reasonable\ntimeframe, while an excessively high rate would be counterproductive by stochastically pushing\nthe model away from stable gradient exploration. Across the tested range, mutation rates of 5–\n10% performed equally well, while 20% often failed. Populations of 10 models were as effective\nas 50, and most simulations converged to a reasonable solution in fewer than 100 generations,\nbut were always continued to 300 generations to obtain the full potential of the solution.\n2. Stimulus dependence\nIt seems like all of these simulations were run with moving bars oriented perpendicular to the\naxis of motion. How dependent are the primitives on the identity of the stimulus? What if you\nused a thicker/large moving feature? Or even widefield motion? I don't think you need to re-run\nevery simulation but some exploration of different stimuli would bolster the claim that the\nmodel is learning computational primitives and not just strange solutions that only work for one\nstimulus (the gap-in-inhibition solution in Fig. 5 comes to mind).\nI agree with the reviewer that this is an important benchmark for DS mechanisms. I have\nincluded exploration of responses to different bar widths and drifting sinusoidal grating (see\ndetailed responses to reviewer 1). These simulations revealed a similar range of algorithmic\nsolutions as observed with moving bars (see new Supplementary Figures 2, 3, 10). I believe this\ngreatly strengthens the claim that these are, in fact, general-purpose DS solutions.\n3. Biological feasibility\nThe paper claims that these computational primitives are biologically feasible in pretty strong\nterms. E.g. the abstract \"All mechanisms are biologically plausible and correspond to known\nphysiological and anatomical motifs.\" I suppose plausibility is subjective but throughout the\nresults there are models that have a distinctly non-biological feel to them, like things with\nnon-monotonic changes in RF size or orientation across the NP axis (Fig. 2, 5), or a single presynaptic\ncell with a huge RF (Fig. 6). Can this be addressed or softened?\nThe main reason for this effect can be attributed to model training. I did not attempt to\ngenerate monotonic function, nor were these required to replicate known solutions such as the\nH&R implementation with kinetically – distinct inputs (see Fig. R3)\nFigure R3. Different implementations of H&R detector with flexible kinetics solution families.\nAs in Figure R2, each row represents an optimal solution from the trained models, allowed to\nvary in input kinetics and lag time, indicating the offset between stimulus arrival to the cell and\ninitiation of the response. Please note that some parameters change monotonically in some\ncases, but often have a prominent peak or trough indicative of slow / fast inputs, respectively,\nsurrounded by members with less different dynamics.\nVisual neuroscientists are well aware that presynaptic populations often include cells with\nwidely differing kinetic properties, as illustrated above.\nLikewise, neurons in the early visual system exhibit markedly diverse receptive field shapes and\nsurround organizations. It is therefore reasonable to expect that circuit wiring could exploit\nthese features to construct a directional detector. The minimal models demonstrate that\nasymmetry in these parameters, aligned with the preferred motion axis, can effectively promote\ndirection selectivity. Importantly, some simulations converged on solutions where asymmetric\ncells were positioned away from the stimulus edge, underscoring the generality of the\nmechanism. This shows that the effect does not depend on highly tuned or contrived\ndistributions -- a point I view as a strength, highlighting the biological plausibility of the solution\nrather than suggesting an artificial construct.\n4. Null / preferred nomenclature\nThe labeling of the \"null side\" and \"preferred side\" seems backwards in places. The text says the\nnull side is the \"position from which the null direction stimulus enters the RF\" which is my\nunderstanding of the existing convention. But Fig. 1a has the null direction stimulus entering on\nwhat is labeled the \"pref side.\" And again in the description of Fig. 5 the text reads \"…stronger\ninhibitory input on the preferred side …\" but again isn't that the null side? That would agree\nwith the general idea that the BL model relies on \"null side\" inhibition – a description included\nin, for example, the Fried 2002 paper cited there and I'm sure other papers as well. Can you\nclear up this discrepancy since it really confused me while reading the paper.\nApologies for the confusion; the figures and the text were fixed with the existing convention.\nMinor\n- Can you add a legend item for the dotted line in 1d (top)\nAdded\n- I'm confused about the description of the model in Fig. 3 – text says \"band-like organization of\nstronger synaptic inputs along the motion preference axis\" but isn't the band perpendicular to\nthe NP axis? Panel b looks like there are stronger inputs on both the N and P side of the RF –\nhow is this similar to the 2 input model and how does it generate DS responses?\nFixed ‘perpendicular’. I have added the following text to the results section to explain and\ncompare the band-like patterns:\n“As in minimal models, two solution types were apparent. About half of the models evolved to\nhave two bands of stronger synapses, with the band on the null side of the DSGC having more\nnumerous inputs, mirroring the optimal solution seen in minimal models. In other cases, double\nstrong-to-weak transitions evolved, each following the logic of the alternative minimal model\nsolution (Fig. 3b, right).”\n- Call out to Fig. 1g – there is no Fig. 1g. I assume this should only refer to Supp Fig. 1?\nFixed\nReviewer #2 (Remarks on code availability):\nI don't see any code in the supplementary materials. I cannot evaluate the code.\nReviewer #3 (Remarks to the Author):\nThis is an important and timely manuscript that systematizes different computational\napproaches to motion selectivity. The manuscript is well written and technically correct. It was\ninteresting to read how motion selectivity can be achieved not only through spatial differences\nin temporal filtering as int he HR detector but also through spatial differences in receptive field\nsize, ratios, and surround properties. One aspect that I can suggest for improvement is that it\nwas not clear how well clustering into different primitives worked, whether clusters were well\nseparated or it was just a continuum between different axes (temporal kinetics, amplitude,\nreceptive field size and asymmetry). In other words, if there is strong difference in orientation\ndifferences does it preclude observing differences in dynamics. Or to pose the question in yet\nanother way: Figure 4B shows the distribution of different possible effects individually but what\nabout correlations between them (combinations of primitives).\nI appreciate the reviewer's suggestion and changed the computational primitives plots\naccordingly (reproduced below).\nNew Figure 4\nNew Figure 7\nIndeed, in some instances, there is evidence of clusters, which tend to correspond to different\ncircuit implementations mentioned in the main text. Excitation-only solutions (Fig. 4) were more\nlikely to be a combination of two primitives (as indicated by models along the diagonals),\nwhereas inhibitory outcomes were almost always restricted to a single mechanism.\nReviewer #4 (Remarks to the Author):\nI appreciate the effort the author has made in exploring computational principles underlying\nelementary motion detection. The manuscript is ambitious in scope, presenting a wide range of\nsimulated circuit architectures and proposing that these converge on a small number of\n\"computational primitives\" for motion detection. The use of minimal, tractable feedforward\nmodels as a way to isolate and examine core features of direction selectivity is a commendable\nstrategy, and I agree with your stated rationale that such simplification can reveal useful\ninsights.\nHowever, despite the technical breadth of the work, I find that the manuscript does not\nconvincingly achieve its conceptual goals, nor does it provide sufficient clarity about the key\nquestions it is addressing. My concerns center around the overall framing and purpose of the\nstudy, rather than any specific technical errors in the simulations. Below, I outline several major\nissues that I believe need to be addressed for the study to be considered suitable for\npublication.\n1. Lack of a Clear Scientific Question or Rationale\nThe manuscript surveys a large number of possible circuit configurations that can implement\ndirection selectivity (DS), but it remains unclear what the main scientific question is. Are you\ntesting specific hypotheses about biological mechanisms? Are you proposing new models for\nexperimental validation? Or are you primarily interested in mapping the design space of\nmotion-detecting circuits? This ambiguity makes it difficult to assess the contribution of the\nwork.\nYou state in the introduction that the study aims to explore how basic receptive field features\nand spatial arrangements give rise to directional tuning. However, the simulations appear\nlargely exploratory and descriptive, lacking a well-defined problem or hypothesis. Without a\nclear rationale, it is difficult to judge the value or relevance of the results beyond their\nexistence.\nThe main goal of this work is to identify the foundational principles that shape receptive field\nstructures underlying direction selectivity. In the revised introduction, I now explicitly frame this\nguiding question. Understanding how receptive field properties influence motion processing is\nnot an esoteric issue, but rather a central and still unresolved question in visual neuroscience.\nHundreds, if not thousands, of studies have been dedicated to this topic. I believe the present\nwork, building on this extensive foundation, contributes new and valuable insights.\nThe reviewer asks why the manuscript cannot be more narrowly focused. The reason is\nprecisely because the question is so foundational. A narrow treatment would miss the broader\nimplications. Addressing this problem requires not only theoretical models but also frameworks\nthat can be tested experimentally, a systematic mapping of the design space of motion\ndetection circuits, and consideration of whether these principles extend to applied domains\nsuch as neuromorphic engineering. Limiting the scope would risk reducing a fundamental\ninquiry to a case study, undermining the generality and impact of the work.\nIn general, theoretical studies are valuable because they can define broad organizing principles\nthat unify mechanisms across systems and guide future experiments. For this reason, I believe a\nbroad scope is not a weakness but a strength of the manuscript.\n2. The Study Feels Preliminary Rather Than Conclusive\nAlthough the manuscript presents numerous models and simulation results, the overall\nimpression is one of an exploratory analysis or early-stage investigation. Many of the findings—\nsuch as the rediscovery of Hassenstein-Reichardt (H&R) and Barlow-Levick (B&L)–like\ncomputations—are already well known and expected. The identification of new combinations or\nvariants of these does not in itself constitute a major conceptual advance unless their biological\nplausibility or explanatory power is convincingly demonstrated.\nThe models are presented in a systematic but somewhat catalog-like fashion, and the paper\noften reads like a compilation of simulation experiments rather than a focused scientific\nnarrative. This is not inherently a flaw, but I believe that, for a journal like Nature\nCommunications, a stronger framing and clearer implications are expected.\nThe \"rediscovery\" of the H&R and B&L mechanisms is not redundant but absolutely essential for\nthe value of this work. There are two key reasons for this.\nFirst, failure to identify these classical models would represent a major flaw. If the approach\ncannot recover well-established DS solutions, there would be no reason to trust it in uncovering\nnovel ones. Demonstrating that the framework reliably rediscovers the canonical mechanisms\nprovides a critical validation step that establishes confidence in its utility.\nSecond, and more importantly, the H&R and B&L models provide benchmarks for excitatory and\ninhibitory strategies of DS. In this project, I intentionally used both simplified motion detectors\n(linear integrations of two inputs) and elaborate multicompartmental DSGC models with\nrealistic presynaptic input structure. In each case, the H&R and B&L solutions serve as reference\npoints that illustrate how established DS mechanisms emerge in the two architectures.\nBy systematically altering single parameters in these models, I then showed that the same\nmodeling framework can give rise to qualitatively different DS implementations. These new\nsolutions are not trivial mixtures or variants of existing models. Rather, they reveal distinct\nalgorithmic strategies for motion detection, ones that fall outside the canonical excitatory and\ninhibitory categories (see more discussion and examples below). The comparison to H&R and\nB&L therefore, does not dilute the novelty of the findings-it underscores it, by demonstrating\nthat the new mechanisms emerge alongside but not within the established frameworks.\n3. Ambiguity of Machine Learning's Role\nYou describe using machine learning to explore a \"vast combinatorial space\" of circuit solutions.\nHowever, the exact nature of the machine learning method, its justification, and its biological\nrelevance are not clearly explained. While evolutionary algorithms can certainly be used to\noptimize network parameters, the study does not make clear what has been learned from this\nprocess, other than the fact that many solutions exist.\nMoreover, the manner in which models are selected and evaluated—based on performance\nmetrics like directional selectivity index (DSI)—follows an ML benchmarking logic that is not\nnecessarily meaningful in the context of biological interpretation. There is little discussion of\nhow these models relate to known circuitry or why certain parameter regimes might be favored\nby real neural systems. This weakens the biological impact of the findings.\nI appreciate the reviewer's concern regarding the role of machine learning in this study. Let me\nclarify both the methodological rationale and the biological relevance.\nFirst, the choice of an evolutionary algorithm is deliberate: this approach is well-suited for\nsystematically exploring the combinatorial explosion of possible circuit configurations while\nremaining agnostic about specific mechanistic assumptions. Far from being a generic\noptimization exercise, this strategy allows one to uncover a range of solutions that traditional\nintuition- or hypothesis-driven modeling might miss. As I stated in response to previous\nconcern, the rediscovery of canonical Hassenstein-Reichardt and Barlow-Levick mechanisms\nvalidates the method and demonstrates that it does not simply produce arbitrary solutions but\ninstead converges on biologically meaningful architectures when they exist in the design space.\nSecond, the study is not claiming that \"many solutions exist\" as a trivial observation. Rather, it\nshows that distinct algorithmic strategies for motion detection emerge from small perturbations\nin model parameters, and that these solutions can diverge qualitatively from known models.\nThis has direct biological significance: it suggests that direction selectivity can be implemented\nthrough multiple, functionally robust mechanisms, potentially explaining the diversity of DS\ncircuits observed across species and cell types (see more below).\nThird, the use of DSI here is not abstract benchmarking. With over 25 years of recording and\nanalyzing single-cell responses in-vivo and in-vitro, I can easily identify realistic synaptic\nresponses. In this work, I dedicated significant effort to finding the optimal optimization scheme\nthat would be biologically relevant. I have resolved an ambiguity in the text to indicate that the\nML did not optimize for a simple DSI metric, as the strategy often leads to very small responses\n(it is easier to have better DSI with small inputs due to the impact of the driving force on large\nEPSPs). The exact equation appears in the methods, and it contains two terms - DSI and\npostsynaptic response amplitude to force the models to produce strong postsynaptic\ndepolarizations. Overall, the fitness function rewards models that have strong and directionally\ntuned signals, which is biologically meaningful. Only after the training was complete, DSI was\nused to report on trained performance.\nWhy DSI? It is a standard and widely used measure of direction selectivity in experimental\nneuroscience. By using this metric, the models are evaluated on the same basis that real\nneurons are characterized.\nFinally, the machine learning component is not an end in itself but a tool: 1) to validate the\nframework against classical solutions, 2) to discover qualitatively new DS mechanisms, and 3) to\nmap parameter regimes onto biologically interpretable circuit features. Far from weakening the\nbiological impact, this approach strengthens it by revealing principles that are generalizable,\nexperimentally testable, and potentially explanatory of cross-species diversity in DS circuits.\n4. Biological Plausibility Is Not Adequately Addressed\nAlthough you note that all mechanisms are \"biologically plausible,\" there is little evidence or\ndiscussion to support this claim. Many of the circuit motifs presented could just as easily be\ndismissed as mathematical possibilities with no relevance to actual neural systems. Without\nexperimental data or stronger theoretical justification, it is impossible to tell which of these\nmechanisms the brain might actually use.\nFurthermore, because the design space is enormous, the models explored here necessarily\noccupy a minuscule region of the total parameter space. This raises concerns about the\nrepresentativeness of the results. If no constraints (biological, theoretical, or functional) are\napplied to restrict the search space, then the significance of any particular solution is unclear.\nI respectfully disagree with the reviewer's assessment that biological plausibility is inadequately\naddressed. The study does not present arbitrary \"mathematical possibilities,\" but instead is\nexplicitly grounded in biologically relevant constraints at multiple levels.\nFirst, the models are constructed from known circuit components of the retina / cortex and\ndirection-selective cells. Inputs are drawn from realistic presynaptic populations, and the\nmulticompartmental DSGC and L2/3 models incorporate experimentally established biophysical\nproperties. The mechanisms that emerge are therefore not abstract mathematical constructs\nbut solutions built on the same substrates available to real neurons.\nSecond, the evolutionary algorithm is not an unconstrained search through an infinite\nparameter space. The design space is intentionally bounded by biological considerations, such\nas synaptic weight ranges, short-term kinetics and receptive field structure, so that all solutions\nremain within physiologically reasonable regimes. This is precisely why the method successfully\nrediscovers the canonical mechanisms: the constraints are sufficient to recover established\nbiology.\nThird, it is important to distinguish between circuit implementations (of which there can be\nmany) and the underlying algorithms. These represent different levels of analysis, much in line\nwith Marr's classic framework. At the computational level, the task is to detect stimulus motion.\nAt the algorithmic level, the question is which computational principles can solve this task. At\nthe implementation level, the focus is on how specific circuit \"hardware\" realizes those\nalgorithms.\nGoing into this project, I anticipated that multiple distinct circuits could rearrange input filters in\ndifferent ways to implement established algorithmic solutions. Demonstrating this is itself\nvaluable, as it provides concrete hypotheses for how diverse neural circuits might arrive at\nsimilar functional outcomes and offers a foundation for future experimental testing.\nThe more significant finding, however, is that the known set of algorithmic solutions to direction\nselectivity is incomplete. This work reveals four new computational primitives that are distinct\nfrom each other and from the canonical H&R and B&L models. These results suggest that\nmotion detection is not underpinned by a small, closed set of strategies, but rather by a broader\nalgorithmic space that biological systems may draw upon in diverse ways. This suggests that DS\ncan be implemented in multiple robust ways, which aligns with the empirical observation that\ndifferent species and even different retinal circuits use diverse strategies for motion\ncomputation. The diversity of solutions is not a weakness but a biologically meaningful insight\ninto the degeneracy and flexibility of neural circuit design.\nFinally, the representativeness concern is directly addressed by the scope of the simulations.\nIn the text, I avoided quantifying the parameter space directly because I think it is an unnecessary\ndistraction. To give a sense of a single aspect of this space, below I show a schematic of different\nRF components that were examined in this study:\nFigure R4. Exponential growth of model complexity with RF components.\nLeft: Illustration of some of the receptive field components considered in this work: both\nexcitatory and inhibitory presynaptic cells can have center and surround components. Each can\ndiffer between presynaptic populations in four parameters: synaptic weight, response kinetics,\nspatial diameter, and degree of circularity. Configurations commonly analyzed in the literature\nare shown in gray. Right: Number of possible circuit architectures assuming each parameter can\neither be shared across or vary between presynaptic cells. This estimate was derived from a\nbinary analysis, yielding: N = 2parametrs . For illustration, if only synaptic weight is\narchitectures\nconsidered, two architectures exist: uniform vs. variable weights. With two parameters, four\ncombinations are possible: shared/shared, variable/shared, shared/variable, and\nvariable/variable.\nNote that in the manuscript, I also examined synaptic short-term plasticity, which adds another\n24 possibilities. Clearly, testing all 220 possible circuit configurations is infeasible within a single\nstudy. Instead, I focused on targeted perturbations of individual model components and\norganized the results in a roughly linear progression and examined questions that were not\nadequately addressed before, including:\n1. How can spatial center RF characteristics impact motion computations in the absence of\nasymmetric kinetics?\n2. Can RF surround mediate DS?\n3. After RF properties are computed in the input population, can the amplitude and\ndynamics of signal transformation at the synapse to the detector affect DS?\n4. How do responses to questions 1-3 change with an added inhibitory presynaptic\npopulation?\n5. Are the newly discovered DS mechanisms relevant when the input is corrupted by noise?\n6. How can they be detected experimentally?\nThus, while the combinatorial design space is indeed vast, the study systematically samples\nacross it and compares outcomes against both canonical models and experimentally relevant\nbenchmarks (e.g., DSI). This approach does not claim to exhaustively map every possibility but\nrather to identify the range of plausible solutions consistent with known biology and to highlight\nnew mechanisms worthy of future experimental testing.\nIn short, the models are biologically grounded, the search space is meaningfully constrained,\nand the results provide both validation against classical mechanisms and discovery of novel,\ntestable principles. Far from being arbitrary, the findings directly inform how neural systems\nmay exploit circuit flexibility to achieve direction selectivity and provide a roadmap to uncover\nthem experimentally.\n5. Lack of Hypothesis-Driven Structure\nA key issue throughout the manuscript is the absence of hypotheses that can be tested or\nevaluated. This is particularly problematic given that many of the findings are not surprising: for\ninstance, that asymmetric RFs or synaptic weights can produce DS. The absence of falsifiability\nor predictive power makes the work more descriptive than explanatory.\nAs you note in the discussion:\n\"I found a strikingly large number of distinct functional organizations of a simple feedforward\nnetwork capable of extracting motion direction from sensory input.\"\nThis sentence encapsulates the core problem: while the variety of possible implementations is\nacknowledged, the manuscript does not provide a framework for determining which of them\nare meaningful. The absence of such a framework reduces the scientific value of the result.\nI again respectfully disagree with the characterization that this work lacks hypothesis-driven\nstructure. The central hypothesis motivating the study is that direction selectivity can arise from\na broader set of circuit implementations and algorithmic solutions than those currently\nrecognized. This is a testable, falsifiable hypothesis: if the framework had recovered only the\nclassical models, then the conclusion would have been that no fundamentally new solutions\nexist or can't be uncovered using this approach. Instead, the discovery of additional\ncomputational primitives directly supports the hypothesis and expands the known solution\nspace.\nThe reviewer is correct that certain findings, such as the fact that synaptic weights can generate\nDS, were reported before. These serve an essential role as positive controls.\nThe reviewer is concerned that the manuscript reports that asymmetric RFs can produce DS. But\nRF asymmetry is key! The existence of asymmetric RFs was not only known but is an absolute\nrequirement. In fact, the motivating question for this work can be reformulated as “what are\nthe asymmetric RF structures that can contribute to DS?” I find that essentially any asymmetry\ncan be exploited by the model to generate DS responses. Many of such solutions are not\nbiologically relevant - and I did not show them here. These include solutions with highly\ndepolarized membrane potentials, models with extremely prolonged – seconds-long -\nresponses, etc. But many more appear to be fully consistent with known physiology, and some\nare simply fascinating. To provide an example:\nMany cortical cells are orientation selective. This classic property was discovered by Hubel and\nWiesel more than 60 years ago, who showed how orientation selectivity arises from elongated\nRFs, which themselves can be created by a simple combination of circular center-surround\ninput.\nWe know that many cortical cells can also produce DS responses. The current work proposes a\nsimple organization scheme to convert inputs whose only asymmetry is in the orientation of\ntheir RFs to DS output at the detector. This result is conceptually as simple as the creation of\noriented RFs by integration of presynaptic inputs from cells with displaced circular RFs, but as\nfar as I know, this link between orientation and DS was not explored before. In the text, I\nmentioned recent work from the Awatramani group who found that the inputs to DSGCs can\nhave oriented RFs (Hanson et al., 2023). Could these RFs contribute to DS? We do not know\nbecause the theoretical underpinning is missing and the key experiment – how are these\noriented RFs organized as a function of the null-preferred axis was not performed.\nFor these reasons, far from being \"merely descriptive,\" the work provides a structured\nframework for identifying, categorizing, and comparing algorithmic solutions to motion\ndetection. The contribution is not only in enumerating many possible implementations, but in\ndistilling them into algorithmic classes that reveal new organizing principles. This directly\ngenerates predictions that can be evaluated experimentally.\nIn short, the manuscript is both hypothesis-driven and predictive:\nHypothesis: DS is supported by implementations and algorithmic strategies beyond the\ncanonical models.\nTest: Apply an unbiased, biologically constrained search to see whether new solutions emerge.\nOutcome: Rediscovery of classical models (validation) and identification of new computational\nprimitives and demonstration how these could be realized in DS circuits.\nPrediction: These circuits/primitives point to experimentally testable mechanisms and broaden\nthe conceptual framework for DS.\nThus, rather than reducing the scientific value, the generality and hypothesis-generating\nstructure of this work substantially increase its impact by providing both novel predictions and a\nprincipled foundation for future experiments.\n6. Presentation Style and Scope\nThe extensive number of figures (1–10), while individually clear, contributes to a sense of\nfragmentation. There is a lack of narrative cohesion tying the results together or building toward\na central conclusion. The manuscript feels more like an organized summary of simulations than\na focused argument.\nThe use of phrases like \"previously unrecognized circuit motifs\" or \"new computational\nprimitives\" may overstate the novelty of what are, in essence, variations on known mechanisms.\nWithout empirical grounding or theoretical analysis, the use of such terms may risk being\nmisleading.\nI understand the reviewer's concern that the different solutions identified in this work might\nappear too similar, essentially variations on a common feedforward framework. It is true that all\nthe models examined are feedforward circuits: they sample visual information through\nspatiotemporal receptive fields at the input layer and transmit the integrated signals to a\npostsynaptic detector via time-varying synaptic inputs. The key differences, however, lie in the\nparameters that define receptive field organization and, in some models, the kinetics of synaptic\nrelease. As elaborated above, these parameter choices give rise to qualitatively distinct\nmechanisms of direction selectivity, several of which are novel and not captured by classical\nformulations.\nWith that said, I appreciate the reviewer for raising this point, as it highlights the importance of\nclarifying these distinctions for readers who may not be experts in this area. In response, I have\nadded a new supplementary figure (reproduced below) that directly contrasts temporal and\nspatial DS solutions, illustrating their distinct computational principles and demonstrating how\neach can be implemented in neuromorphic hardware. I hope this addition makes the\ndifferences between the key implementations clear and strengthens the overall presentation of\nthe work.\nSupplementary Figure 13. Possible implementations of artificial motion detectors.\na, Schematic of a simple neuromorphic motion detector based on Hassenstein-Reichardt\ncorrelator logic. Two photodetectors (p1 and p2) feed into an analog integrator (adder). DS\narises from differential temporal filtering: p2 is subjected to a delay or low-pass filter, producing\nstronger summation when motion engages the detectors in the upward direction. b, Similar\nschematic for a spatially based motion detector. Here, direction selectivity emerges from\ndifferences in photodetector size: p1 has a larger receptive surface area, while no temporal filter\nis applied. As demonstrated in this work, such spatial asymmetry can generate robust upward\nmotion responses with selectivity comparable to the kinetic implementation.\nWhile the technical execution of the simulations is sound, and the general idea of using\nmachine learning to explore motion detection mechanisms is potentially valuable, the current\nmanuscript lacks the conceptual clarity and scientific focus needed to support its claims. I\nencourage the author to reconsider the framing of the work, tighten the narrative around\nspecific questions, and critically assess the biological meaning of the results.\nI hope this feedback is helpful in guiding the next steps for this promising but currently\npreliminary line of research.\nMinor Comments\nThe manuscript lacks page and line numbers, which makes it difficult to reference or discuss\nspecific passages in detail. Including these elements is important for facilitating constructive\nfeedback and precise editorial communication.\nFixed\nThere is inconsistency in figure panel labeling. For example, panels are labeled as \"Ci, Cii\" in the\nfigures but referred to simply as \"1C\" in the main text. This unconventional and inconsistent\nnotation can be confusing for readers and reviewers alike. Please consider standardizing the\nfigure references to align with common scientific writing conventions.\nFixed\nReferences:\nBadwan, B. A., Creamer, M. S., Zavatone-Veth, J. A., & Clark, D. A. (2019). Dynamic\nnonlinearities enable direction opponency in Drosophila elementary motion\ndetectors. Nat Neurosci, 22(8), 1318-1326. https://doi.org/10.1038/s41593-019-\n0443-y\nEuler, T., Detwiler, P. B., & Denk, W. (2002). Directionally selective calcium signals in\ndendrites of starburst amacrine cells. Nature, 418(6900), 845-852.\nhttp://www.ncbi.nlm.nih.gov/entrez/query.fcgi?cmd=Retrieve&db=PubMed&dopt=\nCitation&list_uids=12192402\nFreeman, A. W. (2021). A Model for the Origin of Motion Direction Selectivity in Visual\nCortex. J Neurosci, 41(1), 89-102. https://doi.org/10.1523/JNEUROSCI.1362-\n20.2020\nGruntman, E., Reimers, P., Romani, S., & Reiser, M. B. (2021). Non-preferred contrast\nresponses in the Drosophila motion pathways reveal a receptive field structure\nthat explains a common visual illusion. Curr Biol.\nhttps://doi.org/10.1016/j.cub.2021.09.072\nHanson, L., Ravi-Chander, P., Berson, D., & Awatramani, G. B. (2023). Hierarchical\nretinal computations rely on hybrid chemical-electrical signaling. Cell Rep, 42(2),\n112030. https://doi.org/10.1016/j.celrep.2023.112030\nLien, A. D., & Scanziani, M. (2018). Cortical direction selectivity emerges at\nconvergence of thalamic synapses. Nature, 558(7708), 80-86.\nhttps://doi.org/10.1038/s41586-018-0148-5\nMaisak, M. S., Haag, J., Ammer, G., Serbe, E., Meier, M., Leonhardt, A., Schilling, T.,\nBahl, A., Rubin, G. M., Nern, A., Dickson, B. J., Reiff, D. F., Hopp, E., & Borst, A.\n(2013). A directional tuning map of Drosophila elementary motion detectors.\nNature, 500(7461), 212-216. https://doi.org/10.1038/nature12320\nMorrie, R. D., & Feller, M. B. (2018). A Dense Starburst Plexus Is Critical for Generating\nDirection Selectivity. Curr Biol, 28(8), 1204-1212 e1205.\nhttps://doi.org/10.1016/j.cub.2018.03.001\nPoleg-Polsky, A., Ding, H., & Diamond, J. S. (2018). Functional Compartmentalization\nwithin Starburst Amacrine Cell Dendrites in the Retina. Cell Rep, 22(11), 2898-\n2908. https://doi.org/10.1016/j.celrep.2018.02.064\nWienbar, S., & Schwartz, G. W. (2022). Differences in spike generation instead of\nsynaptic inputs determine the feature selectivity of two retinal cell types. Neuron.\nhttps://doi.org/10.1016/j.neuron.2022.04.012\nWilson, D. E., Scholl, B., & Fitzpatrick, D. (2018). Differential tuning of excitation and\ninhibition shapes direction selectivity in ferret visual cortex. Nature, 560(7716),\n97-101. https://doi.org/10.1038/s41586-018-0354-1"}, {"round": 2, "reviewer_comments": {"reviewer_1": "I thank the author for a clear set of answers to the reviewer questions, including mine. Many concerns of mine were addressed.\nHowever, I’m still troubled by not understanding how the author’s models become direction selective to drifting sinusoids when the inputs have nominally identical kinetic parameters, but different spatial parameters (in addition to the spatial offset). This is main point 3 and minor point 4 in the initial review and rebuttal. The models in question are the simple models Figure S2c. Here, the author agrrees that “purely linear models” would not be able to be direction-selective, but then adds that the spatial models presented here as the simple models have sources of nonlinearity to drive direction-selectivity. (I’m glad we’re on the same page about this fact about linear filtering.) Further, the author asserts that these models do not need to be spatiotemporally oriented and somehow generate direction selectivity in another way.\nIt seems reasonable to me that there could be nonlinearities in these models that generate the direction-selectivity where purely linear filtering (equal to the linear spatial and temporal filtering parameters in the model) would not. However, I don’t think it is sufficient to just say that there is some unidentified nonlinearity at play here. Since this spatially varying unit RF size is portrayed as counter to 40-odd years of dogma on oriented receptive fields, and since the author has complete knowledge of these models, I think it is incumbent on the author to narrow down exactly which nonlinearities are permitting this and how they effect processing to generate the direction-selectivity. (It looks as though the adaptation equation ReLU is the only nonlinearity in the photodetector equations upstream of the linear summation, if I’m understanding the simple model correctly.)\nConcretely, I could easily imagine that this photodetector nonlinearity results in changes in kinetics when the receptive field is larger. If so, then this effect could rather be thought of as a parameterization problem—the model parameters seem to imply that the kinetics are the same for the two channels, but in reality they are not. Such a result would be consistent with the standard oriented RF interpretation; it would only mean that the orientation of the RF is obscured in the model’s parameterization into linear and nonlinear components. If this is what’s going on, then I do not think that this spatial RF size change constitutes a form of direction-selectivity distinct from the prior kinetic ones. Instead, it would be another manifestation of the same kinetic differences, but using recruited nonlinearities to arrive at the kinetic difference. There may be some other way that this model becomes DS; if so, and if the kinetics of the inputs to these stimuli truly are the same, then I do think the results could be framed as distinct from prior models.\nOverall, I find the text in the manuscript a bit less insightful about these issues than the rebuttal. The front-end nonlinearities appear to be central to what’s going on and this should be emphasized, and the nonlinearities explained with care, including their potential effects on kinetics when coupled to changes in spatial receptive field extents. Figure S13 is representative of the issue. The explanation makes only references to nominally linear processing in p1 and p2 here, referring to kinetics and spatial size. However, panel b cannot be direction-selective to sinusoids if all the processing is linear. Where are the nonlinearities? Which ones matter? Why do they get recruited more or differently for large than for small receptive field sizes? How do they generate the direction-selective responses in model S13b? I think these questions are critical in all the models where spatial receptive field differences are said to create DS signals without explicit parameter differences in the kinetics (inhibitory input strength, size, etc.). Are these changes in spatial receptive fields actually creating kinetic differences when the model is presented with sinusoids? If so, it’s a completely different interpretation than saying that circuits may be DS by employing larger receptive fields on one side than the other.\nOne last point on this: the author seems to have refit the models to make them DS when presented with sinusoids. With the moving bar, an entirely linear system can be DS with the max-response metric. It might be informative to compare the two solutions to find out what nonlinearities were invoked in the sinusoidal case.", "reviewer_2": "My major points 1, 2, and 4 were sufficiently addressed.\nMy point 3, about biological feasibility was dismissed with no changes to the manuscript that I can discern. I don't think I'm being unreasonable, and another reviewer also raised a concern about the claims of biological feasibility. To be clear, it isn't that I believe that diverse spatial and temporal tuning properties in upstream neurons is non-biological. But rather that in the present manuscript the DS that is observed in some models appears to rely on very specific, non-monotonic, or otherwise \"strange\" arrangement of those inputs. The author's response suggests that the particular arrangement doesn't matter, and that there are families of solutions that have different arrangements that work perfectly well. If this is the reason for the apparent non-biological flavor of some of the models presented in the main figures then I believe this should be at least demonstrated and explained to the reader, given the paper's strong claims about biological feasibility.", "reviewer_3": "The manuscript has improved. Especially useful was the addition of comparison between moving bars and gratings.\nHowever, in the second reading of the manuscript, the work does strike as delivering new conceptual advance. My take-away message from the manuscript is that directional selectivity can be achieved in any system that possesses asymmetry in in receptive field shape and/or dynamics. If that is correct, perhaps it is useful to state this up front. I also agree with other reviewers concerns, especially those from Reviewer 4.", "reviewer_4": "I appreciate the authors’ careful and detailed responses to the initial review. Below I offer a few additional comments on several key issues.\nFirst, my earlier concern was not about the breadth of the topic per se. It is important not to conflate a “specific” research goal with a “narrow” one. A study can address a large, foundational question and still define a clear, specific objective; conversely, even a very narrow topic can be presented in a vague or diffuse manner.\nSecond, it remains unclear how the reported results can be meaningfully linked to actual brain mechanisms. What concrete experimental approaches could be used to explore such links, and what evidence supports the claim that the reported results are biologically plausible? Stating that the model is “constructed from known circuit components” does not by itself guarantee relevance to real neural systems. In a design that is not centered on testing a sharply defined hypothesis but instead on simulating a moderately realistic network and searching for interesting outcomes, it is difficult to know how specific modeling choices or constraints influence the results. Moreover, the model conditions used here — shallow architecture, random initialization, etc. — differs substantially from biological circuits, which makes strong claims about biological reality hard to accept.\nRelatedly, the manuscript asserts that the work “reveals four new computational primitives distinct from the canonical H&R and B&L models.” Such predictions would indeed be exciting if the model could further lead to experimental observations in the brain. However, producing novel patterns from a partially biologically inspired network is always possible, and without stronger justification it is not clear that other investigators should be motivated to invest in experimental validation.\nThe authors summarize their approach as “both hypothesis-driven and predictive,” with the stated hypothesis that “DS is supported by implementations and algorithmic strategies beyond the canonical models.” As written, this hypothesis is so general that it is effectively unfalsifiable; it would be surprising if any reasonable dataset failed to support it. This makes the study appear less truly hypothesis-driven and more exploratory or post hoc in nature. The outcome — rediscovery of classical models and identification of new computational primitives — does not convincingly advance the field, as validation of classical models has been addressed extensively in prior work. Likewise, the prediction of “experimentally testable mechanisms” leaves open the key question of whether the newly proposed primitives are genuinely expected to exist and be observable in the brain, and if not, why not.\nIn summary, while the modeling framework is interesting, the manuscript would benefit greatly from a clearer articulation of specific, falsifiable hypotheses, stronger justification for the biological plausibility of the proposed mechanisms, and a more rigorous discussion of how the predictions could be tested experimentally."}, "author_response": "I agree with the reviewer that clarifying how non-kinetic solutions generate motion enhancement for drifting sinusoidal gratings is essential. In response, I substantially revised the presentation of these solutions and now include a more explicit, mechanistic comparison of the different computational strategies that can produce direction selectivity in this regime.\nDefinitions used in this response.\nFor clarity, I refer to operations that preserve the sinusoidal form of the response (such as shifting, multiplication, or convolution) as linear. These operations keep each presynaptic response a sine wave (or a sum of sine waves). In contrast, processes that alter the waveform shape are termed nonlinear. For simplicity, I assume that the postsynaptic detector performs a simple linear summation of its presynaptic inputs, and I quantify DS as the difference in peak-to-peak amplitude of this summed postsynaptic response. For example, a fluctuation from -70 mV to -50 mV carries more directional signal than a fluctuation from -61 mV to -59 mV; various downstream mechanisms could further amplify this difference, but such steps are beyond the scope of this discussion.\nNonlinearities in the RF formulation\nThe model includes two sources of nonlinearity:\n• Readily releasable pool (RRP) dynamics - synaptic release depletes the RRP, which in turn reduces subsequent release. This depletion-recovery loop is nonlinear.\n• Center–surround rectification (ReLU-based operation that combines center and surround components). This nonlinearity is absent in models that use the center component alone.\nIn contrast, the photoreceptor drive onto each RF is integrated linearly: the RF template (a 2D Gaussian) is multiplied by the stimulus at each time step. As a result, drifting sinusoidal gratings always generate sinusoidal signals following spatial RF integration; smaller RFs generate larger-amplitude sinusoids, whereas larger RFs generate smaller-amplitude ones, but phase and frequency remain fixed.\nThis behavior contrasts with moving bars, where changes in RF size create differences in response kinetics. With a bar, small RFs experience a brief, transient activation while large RFs integrate the stimulus over a longer time window.\nSolutions for drifting gratings\nDo drifting gratings require nonlinear computations to achieve DS?\nNot necessarily. Classical temporal models can generate DS purely through phase shifts. For example, presynaptic cells with fixed delays produce sine waves of identical frequency but offset phase. Linear summation of two such inputs with different delays is direction selective, a well-known solution that can be implemented compactly (e.g., code accompanying (Borst, 2018)).\nIn my model, RF size or amplitude alone cannot directly shift the temporal phase. However, linear solutions remain possible via center-surround interactions. When center and surround RFs are allowed to differ in their temporal activation timescales, their linear combination yields direction-dependent phase relationships, even without invoking the nonlinear components. These mechanisms are presented in the new Supplementary Figure 6 and described in the revised text. As seen in that figure, linear center-surround models achieved DSIs of ~15–25%, compared to ~50% when nonlinearities (RRP depletion and rectification) are present. In comparison, models in which center and surround share identical kinetics perform worse, even when RRP nonlinearity is retained (DSI ≤10%).\nTo eliminate ambiguity about the linear interactions, I extended the Borst model to include explicit spatial RFs and surrounds. The full implementation (<200 lines of code; roughly half is plotting) is available at:\nhttps://github.com/PolegPolskyLab/DS-mechanisms/tree/Dec2025\nminimal_example.py\nRegarding Figure S13 (now S14): this figure is intentionally schematic. Its purpose is to illustrate the conceptual distinction between the sensor (photodetector/RF) and the connection to the detector (wire/synaptic release). It is not meant to depict biophysical detail but rather to highlight the key point that, although the signals arriving at the detector must differ in their temporal envelopes to generate DS, the specific implementation at the sensor level can vary widely.\nThe figure serves as a simplified, didactic contrast (an artificial “instrument analogy”) to communicate this high-level concept. Importantly, the revised manuscript now includes a much more rigorous and explicit mechanistic analysis of the linear and nonlinear components that generate DS (Figs. S5–S6). These new analyses directly address the reviewer’s concerns by demonstrating precisely how different RF and synaptic components contribute to temporal asymmetries.\nGiven this expanded, detailed treatment of the mechanisms, I believe Fig. S14 continues to serve its intended clarifying role and can remain unchanged.\nConclusion\nBoth linear and nonlinear architectures can generate DS for moving bars and drifting sinusoidal gratings.\n• Kinetic models shift the response phase directly and therefore do not require nonlinear elements, although, as shown in Figs. S5-S6, nonlinearities substantially improve DSI.\n• Models constrained to vary only RF amplitude or spatial extent cannot directly manipulate phase, but can nonetheless generate DS through linear center-surround interactions when the two components have distinct temporal properties.\n• Biologically grounded nonlinearities (RRP dynamics, center-surround rectification) significantly enhance DS but are not strictly required for the computation.\nReviewer #2 (Remarks to the Author):\nMy major points 1, 2, and 4 were sufficiently addressed.\nMy point 3, about biological feasibility was dismissed with no changes to the manuscript that I can discern. I don't think I'm being unreasonable, and another reviewer also raised a concern about the claims of biological feasibility. To be clear, it isn't that I believe that diverse spatial and temporal tuning properties in upstream neurons is non-biological. But rather that in the present manuscript the DS that is observed in some models appears to rely on very specific, non-monotonic, or otherwise \"strange\" arrangement of those inputs. The author's response suggests that the particular arrangement doesn't matter, and that there are families of solutions that have different arrangements that work perfectly well. If this is the reason for the apparent non-biological flavor of some of the models presented in the main figures then I believe this should be at least demonstrated and explained to the reader, given the paper's strong claims about biological feasibility.\nI agree with the reviewer’s point. Accordingly, I modified most of the example implementations so that the simulations now produce monotonic functions. Achieving monotonicity is straightforward in models where only one or a few parameters are allowed to vary. However, when the parameter space becomes large, models frequently converge on non-monotonic functions. This behavior is expected, for the reasons outlined below.\nAll RGC models in the manuscript are innervated by four distinct groups of excitatory cells and, in some cases, an additional four groups of inhibitory presynaptic cells. When no explicit constraints are imposed, each presynaptic group can evolve independently.\nIn such a high-dimensional parameter space, a given algorithmic solution can be realized in multiple, equally valid structural implementations. For example, in several models the optimal solution was effectively to silence a subset of inputs. The models can accomplish this in several interchangeable ways, including: (1) generating strong, prolonged surround supression that shunts center activation; (2) evolving center and surround RFs with matching kinetics and spatial parameters, thereby canceling each other; or (3) reducing the synaptic weights from those presynaptic cells to zero.\nBecause these presynaptic populations operate independently, any of these solutions may arise randomly yet produce the same functional outcome. Furthermore, when synaptic weights for a presynaptic group are driven to zero, the remaining RF parameters of that group become arbitrary, as they no longer influence the postsynaptic cell’s response. Consequently, in sufficiently complex models, examining a single parameter dimension (e.g., RF kinetics alone) may not yield a monotonic dependence. Instead, the appropriate approach is to evaluate the model's complete functional output, for instance, the postsynaptic cell’s final response. When viewed at this level, the underlying computational logic becomes clear and follows the computational primitives described in the main text.\nOverall, among the ~60 circuit implementations presented in the manuscript, only about five exhibit non-monotonic parameter dependencies, all of which arise from the complex high-dimensional scenarios described above.\nReviewer #2 (Remarks on code availability):\nI could not get main.py to run. I suggest making a requirements.txt file for installation, or making the respository an installable python package, and adding a readme with installation and basic entry point instructions.\nBoth requirements.txt and readme.txt should be present. NEURON can be finicky to install. I suggest using Google Colab, where the installation of NEURON always works.\n#install neuron for python\n!pip install neuron -q\n# clone the github\n!git clone --branch Sep2025 --recurse-submodules https://github.com/PolegPolskyLab/DS-mechanisms.git\n#install requirements\n%pip install -U pip\n%pip install -r requirements.txt\n# compile the neuron mod files\n%cd DS-mechanisms/mod\n!nrnivmodl\n# copy the x84_64 folder to the main directory\n!mv -v x86_64 ..\n%cd ..\n# this should now work!\n!python main.py\nIn addition, I created a simplified version of the model based on (Borst, 2018), showing how DS can be created in different RF configurations in response to drifiting gratings stimulation. It can be accessed from:\nhttps://github.com/PolegPolskyLab/DS-mechanisms/tree/Dec2025\nminimal_example.py\nReviewer #3 (Remarks to the Author):\nThe manuscript has improved. Especially useful was the addition of comparison between moving bars and gratings. However, in the second reading of the manuscript, the work does strike as delivering new conceptual advance. My take-away message from the manuscript is that directional selectivity can be achieved in any system that possesses asymmetry in in receptive field shape and/or dynamics. If that is correct, perhaps it is useful to state this up front. I also agree with other reviewers concerns, especially those from Reviewer 4.\nThis is accurate and I modified the text to reflect this insight better\nReviewer #4 (Remarks to the Author):\nI appreciate the authors’ careful and detailed responses to the initial review. Below I offer a few additional comments on several key issues.\nFirst, my earlier concern was not about the breadth of the topic per se. It is important not to conflate a “specific” research goal with a “narrow” one. A study can address a large, foundational question and still define a clear, specific objective; conversely, even a very narrow topic can be presented in a vague or diffuse manner.\nSecond, it remains unclear how the reported results can be meaningfully linked to actual brain mechanisms. What concrete experimental approaches could be used to explore such links, and what evidence supports the claim that the reported results are biologically plausible? Stating that the model is “constructed from known circuit components” does not by itself guarantee relevance to real neural systems. In a design that is not centered on testing a sharply defined hypothesis but instead on simulating a moderately realistic network and searching for interesting outcomes, it is difficult to know how specific modeling choices or constraints influence the results. Moreover, the model conditions used here — shallow architecture, random initialization, etc. — differs substantially from biological circuits, which makes strong claims about biological reality hard to accept.\nRelatedly, the manuscript asserts that the work “reveals four new computational primitives distinct from the canonical H&R and B&L models.” Such predictions would indeed be exciting if the model could further lead to experimental observations in the brain. However, producing novel patterns from a partially biologically inspired network is always possible, and without stronger justification it is not clear that other investigators should be motivated to invest in experimental validation.\nThe authors summarize their approach as “both hypothesis-driven and predictive,” with the stated hypothesis that “DS is supported by implementations and algorithmic strategies beyond the canonical models.” As written, this hypothesis is so general that it is effectively unfalsifiable; it would be surprising if any reasonable dataset failed to support it. This makes the study appear less truly hypothesis-driven and more exploratory or post hoc in nature. The outcome — rediscovery of classical models and identification of new computational primitives — does not convincingly advance the field, as validation of classical models has been addressed extensively in prior work.\nWhat advances the field (among others) is the identification of new computational primitives.\nLikewise, the prediction of “experimentally testable mechanisms” leaves open the key question of whether the newly proposed primitives are genuinely expected to exist and be observable in the brain, and if not, why not.\nIn summary, while the modeling framework is interesting, the manuscript would benefit greatly from a clearer articulation of specific, falsifiable hypotheses, stronger justification for the biological plausibility of the proposed mechanisms, and a more rigorous discussion of how the predictions could be tested experimentally.\nDS models:\nIt is possible to make a model of the DS system that is complex enough, which would have multiple solutions that are not particularly interesting because they lack biological relevance. For example, deep neural networks can reach multiple solutions because of the huge parameter space one can cover, but the majority of these solutions will have nothing to do with biology and instead reflect abstract mathematical properties of the network.\nDeep neural networks can approximate any function, this generality is precisely why they are not well-suited for the present question. Demonstrating that a sufficiently deep network can compute DS would add little insight: such networks can implement virtually any mapping. Instead, the models I propose here are concise and highly constrained throughout the main text of the manuscript. I do explore “unconstrained” shallow models – but these are only mentioned briefly and presented in supplementary materials for a reason: there is not that much information one can glean from these complicated solutions.\nConfusingly, going to deep networks is exactly what the reviewer is suggesting. If the reviewer means to explore a fully biologically realistic deep model that represents processing in the retina, LGN, V1, and beyond, this is an enormously complex undertaking, well beyond the scope of the current work and in practice only feasible with resources on the scale of the Blue Brain project. More importantly, as I showed in the first rebuttal, even in a simplified two-layer system, the solution space is already extremely rich; adding additional layers would expand this space exponentially and obscure, rather than clarify, the underlying computational principles.\nI chose a shallow model architecture for two main reasons. First, it aligns with the classical frameworks of Hassenstein & Reichardt and Barlow & Levick, which operate at the same level of abstraction. Contrary to the reviewer’s suggestion, a shallow architecture is an appropriate approximation for the retinal circuit composed of bipolar cells synapsing onto direction-selective ganglion cells, and it captures the essential computational transformations known to occur at this stage. The revised main text now includes this information.\nSecond, in contrast to the reviewer’s assertion, the models I am presenting are experimentally tractable. My organization of the manuscript follows logically from the RF composition and lends itself to simple experimental validation that requires only measurement of the RF structure. Another major advantage of simplifying, rather than complicating the analysis, is the concrete, testable model predictions (see more about these critical points below).\nThe reviewer complains about randomized parameter initialization. I think they are missing a key point. All the parameters are still constrained within the biologically relevant range – it is only where within the range they are instantiated that the randomization is applied (this is now explicitly mentioned in the revised main text).\nNaturally, I spend much effort to describe the solutions reached in the simulations. That desn’t mean that such solutions are common. The opposite is true. DS in the first generation, composed of randomly initialized models, is invariantly low. Effective solutions occupy a tiny fraction of the parameter space; this is true for classical models and for the new computational primitives and circuit implementations described here.\nExperimental evidence:\nThe central conceptual finding of this project is that basic RF components beyond kinetics can robustly support DS computation. The manuscript provides specific, testable predictions. For example, the models predict that DS should be enhanced when presynaptic RF sizes increase progressively along the sequence of activation in the preferred direction. This prediction emerged consistently across all models constrained to vary only in RF center size, even when initialized with random values - strong evidence that this configuration represents a uniquely optimal solution under those constraints. Random initialization was essential for discovering this convergence.\nTo address the reviewer’s concern about biological plausibility, I now cite (Matsumoto et al., 2019) in the revised Discussion. In that study, the authors examined the temporal dynamics of bipolar cell inputs to direction-selective ganglion cells and reported clear differences in temporal kinetics consistent with a Hassenstein-Reichardt-type computation. As part of their basic characterization of visual response properties, they also measured the spatial extent of the receptive fields. They briefly mention that RF size increases along the preferred direction of activation. This observation was not pursued further in their work, as the prevailing view in the field holds that RF size does not contribute to DS tuning.\nFor clarity, I reproduce their Fig. 2F below. As shown, there is both a sharpening of the temporal RF and a widening of the spatial RF in the preferred direction (upward in their convention, matching the convention used in my work). Note that the yellow signal most likely originates from glutamatergic amacrine cells, and mediates a small fraction of the excitatory drive to the cell.\nI compare these findings directly with predictions from my models. In panel (a), the left and middle plots correspond to models in which only kinetic parameters or only spatial RF size can vary (as shown in Fig. 1). The right plot shows a model in which both parameters can vary simultaneously. In panel (b), I present the resulting DSI values across stimulus velocities. As predicted, the combined model performs better than models restricted to varying only kinetics or only RF size.\nThe alignment is striking: the experimental results in (Matsumoto et al., 2019) mirror the qualitative and quantitative predictions of my models.\nWith this information and given the long history of receptive-field mapping in visual neuroscience, I am unsure why the reviewer finds the proposed “concrete experimental approaches” unclear. Techniques to measure RF size, structure, and dynamics are standard and widely used, and I outline how well-established RF-mapping approaches, used in vision science for over 70 years, can be employed to test these simple predictions experimentally (Fig. 9).\nTestable predictions and hypotheses:\nThe familiar maxim “all models are wrong, but some are useful” is particularly relevant here. Every theoretical model necessarily simplifies biological reality. In this work, I present models at two different levels of abstraction and across two different systems.\nAs detailed above, I have addressed the reviewer’s concerns regarding experimental validation. However, contrary to the reviewer’s suggestion, I believe the simplified models are often more valuable than attempts to fully replicate biological complexity, an undertaking that may not even be achievable in practice. Simplified models allow for clarity, systematic dissection, and direct insight into the underlying computational principles.\nFor instance, I predict that an anti-Barlow & Levick (anti-B&L) algorithm exists in the brain. The computation itself is conceptually straightforward: inhibition that precedes excitation in the preferred direction. Observing this computation in highly simplified models is not a weakness; rather, it demonstrates the generality and robustness of the underlying algorithm.\nThis perspective is consistent with the history of the canonical B&L model. When first described in the mammalian retina (Barlow & Levick, 1965), the model operated at a comparably abstract level. Only later did experimental work reveal that starburst amacrine cells exhibit the spatially offset RFs and slower dynamics predicted by the model. Importantly, the same conceptual algorithm applies across diverse systems, including mammalian cortex and the Drosophila visual system, despite the absence of starburst cells in those circuits. The implementation differs, but the computation is conserved (see instance (Borst & Helmstaedter, 2015)). Thus, models that are not tied to detailed biological realism can still be highly informative and broadly applicable.\nMy goal in this work is to discover novel algorithms and implementations; naturally, there is not yet direct experimental evidence for all of the models presented. That does not diminish their scientific value. To suggest that predictions should only be studied once experimental evidence already exists would be to eliminate the role of theory in advancing neuroscience. The reviewer writes that “such predictions would indeed be exciting if the model could further lead to experimental observations in the brain.” But this is precisely the point: the model generates experimentally testable predictions.\nBelow, I list several of these concrete, falsifiable predictions along with the experimental approaches (using established and widely available methods) that could be used to test them:\nPrediction Example experimental approach\nTo maximize DS, RF size should increase Model validated experimentally by\nin the preferred direction recording glutamatergic signals to\ndirection selective ganglion cell.\nSynaptic weights by themselves can Model validated experimentally by\npromote DS in two distinct configurations. postsynaptic patch clamp recordings\nDifferences in RF surround size, kinetics, Surround measurements can be\nand amplitude can promote DS in two performed directly from presynaptic cells\ndifferent implementations each, and these or using analysis of glutamate signals\nimplementations require different onto the postsynaptic cell.\ncombinations of RF properties\nA shallow network with excitatory inputs Calcium or glutamate imaging to record\nonly can combine the responses of the the spatiotemporal patterns of the\npresynaptic inputs in four distinct ways, or excitatory input.\nin other words, four algorithms exist for\nexcitation-only DS (one was known\npreviously).\nIn a system with an inhibitory drive, the Surround measurements can be\namplitude, kinetics and size of the center performed directly from presynaptic cells\nand surround RF components of the (using calcium imaging) or by analyzing\ninhibitory drive can be combined in GABAergic signals onto the postsynaptic\nspecific ways (which I explore at length in cell.\nthe manuscript) to generate DS\nAt least two new algorithmic solutions to These algorithms can be tested by\nhow excitation and inhibition can interact comparing the temporal dynamics of the\nto generate DS. excitatory and inhibitory drives in\npreferred and null directions using\npostsynaptic voltage clamp or calcium\nimaging of presynaptic (E and I)\npopulations.\nFurther, I respectfully disagree with the assertion that an explicitly framed falsifiable\nhypothesis is essential to appreciate the contribution of this work. Many influential\ntheoretical studies in neuroscience were not organized around a single, narrowly\nformulated hypothesis, yet they clearly generated testable (and ultimately falsifiable)\npredictions.\nA classical example is Rall’s work on dendritic integration (e.g., (Rall, 1967)), in which\nmathematical models were used to analyze signaling dynamics in cells with dendrites.\nRall derived numerous predictions about attenuation and filtering in passive dendrites,\nmany of which were later confirmed experimentally. It is not obvious what a concise\n“falsifiable hypothesis” would have been for this entire body of work, beyond something\ntrivial along the lines of “dendrites influence somatic responses,” which adds little\nconceptual value and does not capture the depth of the theoretical advances.\nSimilarly, in the next generation of computational studies inspired by Rall, Mainen &\nSejnowski (Mainen & Sejnowski, 1996) applied cable theory to models with active\nconductances and showed that differences in dendritic morphology can produce\nmarkedly distinct firing patterns, even when the underlying ion channels are similar. The\nkey insight is that morphology and compartmental integration shape nonlinear\ninteractions with voltage-gated channels. Here again, while their simulations produced\nmany testable predictions, the overarching “hypothesis” might be summarized as\n“morphology affects function,” which, as the reviewer notes in another context, is so\nbroad that it is essentially unfalsifiable.\nIn all of these examples, the central scientific question is not “whether” but “how.” Long\nbefore Rall, neuroscientists knew that dendrites affect signal spread. Mainen &\nSejnowski were not recognized because they were the first to propose that dendritic\nstructure can influence action potential generation, but because they elucidated specific\nmechanisms by which it does so.\nIn the same spirit, my manuscript focuses on how direction selectivity can arise from\nparticular receptive-field structures and circuit motifs. I provide concrete examples of\nhow newly identified algorithms can be implemented in biological circuits, offer higher-\nlevel conceptual frameworks for understanding them, and specify the detailed\nrequirements for these computations. All of this is expressed in precise mathematical\nand computational terms that yield clear, testable predictions, some of which already\nalign with existing data. In this sense, the work follows the tradition of theoretical studies\nthat advance mechanistic understanding even when a single, neatly packaged\nfalsifiable hypothesis is not the primary organizing principle.\nMore broadly, theoretical work plays a crucial role in identifying candidate mechanisms\nthat may otherwise remain unexplored. Publishing these predictions provides the field\nwith concrete hypotheses, quantitative expectations, and clear experimental readouts\nthat can motivate and guide empirical studies. Without such a theoretical framework,\nmany of these potential mechanisms would simply not be considered or tested. In this\nsense, the work is intended to stimulate and support future experimental efforts, not\nreplace them.\nIn summary, the manuscript provides multiple testable predictions, at both the\nimplementation and algorithmic levels, regarding DS computation. Available\nexperimental data already support several of these predictions. Furthermore, I outline\nhow well-established, widely used RF-mapping methods can be used to test the\nadditional predictions generated by this work.\nReferences:\nBarlow, H. B., & Levick, W. R. (1965). The mechanism of directionally selective units in\nrabbit's retina. J Physiol, 178(3), 477-504.\nhttp://www.ncbi.nlm.nih.gov/entrez/query.fcgi?cmd=Retrieve&db=PubMed&dopt=\nCitation&list_uids=5827909\nBorst, A. (2018). A biophysical mechanism for preferred direction enhancement in fly\nmotion vision. PLoS Comput Biol, 14(6), e1006240.\nhttps://doi.org/10.1371/journal.pcbi.1006240\nBorst, A., & Helmstaedter, M. (2015). Common circuit design in fly and mammalian\nmotion vision. Nat Neurosci, 18(8), 1067-1076. https://doi.org/10.1038/nn.4050\nMainen, Z. F., & Sejnowski, T. J. (1996). Influence of dendritic structure on firing pattern\nin model neocortical neurons. Nature, 382(6589), 363-366.\nhttps://doi.org/10.1038/382363a0\nMatsumoto, A., Briggman, K. L., & Yonehara, K. (2019). Spatiotemporally Asymmetric\nExcitation Supports Mammalian Retinal Motion Sensitivity. Curr Biol.\nhttps://doi.org/10.1016/j.cub.2019.08.048\nRall, W. (1967). Distinguishing theoretical synaptic potentials computed for different\nsoma-dendritic distributions of synaptic input. J Neurophysiol, 30(5), 1138-1168.\nhttp://www.ncbi.nlm.nih.gov/pubmed/6055351\nhttps://www.physiology.org/doi/pdf/10.1152/jn.1967.30.5.1138\nI thank all the reviewers for their time and effort. I believe the critiques greatly improved\nthe work."}, {"round": 3, "reviewer_comments": {"reviewer_1": "The author has done a fine and careful job addressing my comments. I especially appreciated Figs S4 and S5, which lay out clear where direction selectivity arises for sinusoids. I have only minor comments.\nMinor comments\n1) Overall, I wonder if there are a few places where this paper might benefit by adopting the terminology of separable/inseparable linear filtering and oriented space-time filter – the terminology from Adelson and Bergen that has been so helpful in that part of the field. I’m thinking especially of places like around line 186, where my reading is that you are saying (for the sinusoids) that you can create oriented space-time RFs by combining multiple inseparable space-time filters. (Though it’s more complicated because of the reLU on each c-s synaptic input, so the simple linear analysis of the motion energy model doesn’t necessarily hold.) For the simplified linear models, which are equivalent to the ME linear filtering step, this terminology might be especially apt.\n2) Line 150: “anti-HR” is being defined as the neurons with longer responses on the null-direction side of the RF. But kinetics are identical, so the longer responses are stimulus dependent. This only looks “anti-HR” with bars as the stimulus. With sinusoids, it would not. Shouldn’t the listed properties be properties of the circuit, not the stimulus-circuit pairing?\n3) Last, I remain a tad incredulous at the claim near line 404 that these numerical models are “simpler” than classical frameworks like the HRC and ME. The reasoning given is that the numerical models are simpler because each model configuration required independent testing. I cannot see how one follows from the other. I think there are many things you could say about the differences between these numerical models and the kind of analytically tractable models represented by the HRC and ME. The HRC and ME models are not biologically realistic, they are in fact overly simplified, so that the nonlinearities and phases and anti-symmetric subtractions are prescribed, rather than allowed to be fit parameters. For instance, the original HRC model, proposed in the 1956 paper, has two parameters: the distance between sensors and the time delay. It has a very restricted architecture to make up for the lack of parameters and to focus explicitly on the pairwise correlations that are the simplest hallmark of motion. However, I don’t think one can argue that it’s not simple. Overall, I think it would be worthwhile to spell out precise differences and how they contributed to the new findings, rather than arguing that these numerical models (with nonlinear synaptic release dynamics and integrated neural morphology!) are somehow just ‘simpler’, which I think is hard to define, hard to justify, and an unnecessarily controversial claim to make.\nMinor minor comments\n1) Line 17 in abstract: “previously undescribed” instead of “newly discovered”?\n2) Line 37: “computational motifs” rather than principles?\n3) Line 79: “nonlinear RRP dynamics”? It might be worthwhile to emphasize linear and nonlinear operations throughout…\n4) Line 115 and beyond: RF size gradient with identical positive only kinetics gets bar DS, but NOT sinusoidal DS. Perhaps some telegraphing of this later caveat is order?\n5) Line 315: for HR+BL formulation, probably want to also cite Leong … Clandinin 2016 JNeuro.\n6) Line 386: caveat also: DS is for specific stimuli; paper now gives examples for detectors that work well for bars but not at all for sinusoids. This seems like a point to highlight. What stimuli should we all be using to measure DS? More of them?\n7) Line 378: I didn’t see a ‘Second, …’ for enumerating clearly the second key insight.\n8) Line 452: Author suggests that sinusoids are more complex than simpler bars. However, these are just two representations of space (regular and Fourier) and I’m not sure there’s a good reason to claim that one is simpler or more complex than the other. If anything, most of the basic motion detection behaviors seem to be for detecting flow fields and stabilizing eyes and bodies – extended stimuli seem better suited for investigating/optimizing such detectors.\n9) Line 456: Gruntman did not observe DS inhibition, since that does not appear to be a component of the fly DS circuit, so shouldn’t be cited in this list.\n10) Supp 14: In b, one could mention that this works as drawn for some stimuli (like bars) but not others, like sinusoids, unless you have spatiotemporally inseparable filters for p1 and p2? Consistent with Supp Fig 5.", "reviewer_2": "All of my concerns have been addressed. Thanks for the thoughtful responses.", "reviewer_3": "The author have adequately addressed the previous comments. The revised discussion clearly describes how the computational predictions made here can be tested in future studies, and how the existing experimental evidence on receptive field asymmetries provides early support for the predictions made here.", "reviewer_4": "The efforts made by the author to address the concerns raised by the reviewers are appreciated. While the revised manuscript and rebuttal may give the impression of engaging with these issues, many of the added explanations appear to stem from a misunderstanding of the underlying questions and are therefore largely irrelevant to the concerns raised.\nAs emphasized by other reviewers, the core issue remains the need to clearly establish the biological relevance of the findings — specifically, how the model’s predictions relate to actual biological principles. While the model may produce simulation results that align with experimental data, this does not necessarily imply that it accurately reflects biological mechanisms, nor does it automatically confer scientific significance. Such significance typically arises from the formulation and validation of clearly falsifiable hypotheses, and this aspect appears to remain underdeveloped in the current manuscript."}, "author_response": "I made the chages to address all the minor concerns raised by the reviewer"}]

Decision letter