Introduction

Starburst amacrine cells play a crucial and well-established role in the computation of direction-selectivity (DS) in the mammalian retina (Demb, 2007; Mauss et al., 2017; Wei, 2018). SACs can be broadly categorized into two subtypes: ON and OFF cells, both of which receive excitatory inputs from diverse populations of bipolar cells. In mice, the majority of input synapses onto SACs are concentrated in the proximal two-thirds of the radially symmetric dendritic tree, while the distal dendrites house the output synapses (Ding et al., 2016; Vlasits et al., 2016). Individual SAC branches exhibit robust calcium responses to objects moving in the outward direction (away from the soma); consequently, different parts of the dendritic tree demonstrate distinct directional preferences (Chen et al., 2016; Ding et al., 2016; Euler et al., 2002; Koren et al., 2017; Morrie & Feller, 2018; Poleg-Polsky & Diamond, 2016; Poleg-Polsky et al., 2018). SACs are essential for generating DS in direction-selective ganglion cells (DSGCs); they provide a combination of cholinergic excitation and GABAergic inhibition that specifically target co-stratifying DSGCs sharing the directional preference of the releasing site (Briggman et al., 2011; Fried et al., 2002; Kostadinov & Sanes, 2015; Lee et al., 2010; Morrie & Feller, 2015; Pei et al., 2015; Sethuramanujam et al., 2021; Soto et al., 2019; Wei et al., 2011; Yonehara et al., 2011; Yonehara et al., 2016).

Numerous mechanisms have been proposed to contribute to the establishment of direction selectivity in SACs. For example, network mechanisms, such as feedback inhibition onto neighboring SACs or SAC– derived cholinergic excitation of presynaptic BCs, were shown to enhance DS (Ding et al., 2016; Hellmer et al., 2021; Munch & Werblin, 2006; Poleg-Polsky et al., 2018). However, it is worth noting that SACs maintain their directional tuning in the absence of cholinergic and GABAergic transmission (Chen et al., 2016), indicating that other mechanisms must contribute to directional tuning. Other proposals emphasized the involvement of cell-autonomous processes that encompass dendritic voltage gradients sustained by tonic excitation from BCs (Gavrikov et al., 2003; Hausselt et al., 2007) and fine-tuned morphology that imposes membrane filtering of postsynaptic signals (Tukker et al., 2004; Vlasits et al., 2016). While these mechanisms can contribute to the development of direction selectivity to some degree, they fall short in predicting the experimentally observed levels of directional tuning. In addition, it remains unclear whether these requirements for tonically active presynaptic release described in the rabbit are fulfilled in the mouse circuit, and cells with abnormal shapes still exhibit near-normal levels of direction selectivity (Morrie & Feller, 2018).

Recently, a reformulation of the classical Hessenstein-Reichardt correlator model has been proposed, which emphasizes the kinetic properties of glutamatergic inputs sampled by SAC dendrites as contributing factors to DS (Greene et al., 2016; Kim et al., 2014). Known as the "space-time wiring" model, this formulation is based on connectomic reconstructions of bipolar–SAC contacts. Both ON- and OFF-SAC dendrites tend to stratify closer to the middle of the inner plexiform layer (IPL) as they branch out from the soma. This unique stratification pattern enables SAC dendrites to receive inputs from different types of BCs in their proximal and distal regions (Ding et al., 2016; Greene et al., 2016; Kim et al., 2014). Importantly, functional studies have uncovered a positive correlation between the transiency of BC axonal release and the distance from the boundaries of the IPL (Baden et al., 2013; Chen et al., 2014; Euler et al., 2014; Franke et al., 2017). Specifically, BC types that stratify closer to SAC somas (OFF-SACs: BC1-2; ON-SACs: BC7) exhibit prolonged release, whereas more distal types (OFF-SACs: BC3-4; ON-SACs: BC5) tend to respond with a transient pattern. According to the space-time wiring hypothesis, these regional differences facilitate a more effective summation of excitatory inputs in the outward direction.

Two main methodologies were taken to examine the existence of the space-time wiring hypothesis in the bipolar-SAC circuit. First, with the caveat that the long, thin dendritic processes in SACs impose significant challenges for an effective space clamp (Poleg-Polsky & Diamond, 2011; Stincic et al., 2016), it is possible to measure the distance-dependence of BC dynamics by analyzing somatic excitatory postsynaptic currents. Although one study reported substantial differences in the shape of the recorded currents, primarily in the OFF-SACs (Fransen & Borghuis, 2017), a second group that recorded from ON-SACs failed to replicate this finding (Stincic et al., 2016).

A second method utilized genetically encoded fluorescent glutamate sensors, such as iGluSnFR, to directly study the dynamics of BC release (Borghuis et al., 2013; Marvin et al., 2013). An examination of glutamate release dynamics confirmed the diverse signaling patterns within the lamina of the inner plexiform layer accessible to SAC dendrites (Franke et al., 2017; Gaynes et al., 2022; Strauss et al., 2022). Furthermore, excitatory drive to ON-DSGCs, which share many glutamatergic inputs with ON-SACs, was found to consist of units with distinct release dynamics, theoretically supporting the type of directional computation proposed by the space-time wiring hypothesis (Matsumoto et al., 2019). Finally, a project published last year described the asymmetric distribution of kinetically distinct iGluSnFR signals in ON-SACs (Srivastava et al., 2022).

However, as the investigation of BC signaling has predominantly focused on release dynamics to static stimulus presentations, the question arises: Can responses to flashed stimuli accurately represent the release waveforms produced by moving objects?

The receptive field engagement in response to a continuously moving stimulus follows a sequential pattern, with the surround being activated first, followed by the center. As a result, the motion response profile becomes a complex function, influenced by factors like RF size, temporal characteristics, and stimulation velocity (Strauss et al., 2022). This integration mode stands in stark contrast to the simultaneous activation of RF components observed during flashes of static objects. Indeed, our research has revealed a significant distinction in BC release between stationary and moving probes (Gaynes et al., 2022). This discrepancy emphasizes the limited capacity of flash responses to predict the fundamental characteristics of BC responses to motion accurately and highlights the importance of examining directional capabilities in SACs with physiologically relevant stimuli.

An additional crucial factor to consider is the dependence of the presynaptic response waveform with respect to the direction of activation. The space-time wiring model posits that BCs generate similar glutamatergic drive in both outward and inward directions. Yet, this assumption may be challenged when moving objects appear within the bipolar RF. Signal processing in such scenarios fundamentally differs from stimuli that are flashed or continuously moving, and these differences can vary significantly among different types of BCs (Gaynes et al., 2022; Strauss et al., 2022). The intricacies of signal transformation within the BC network give rise to important considerations regarding the application of advanced visual protocols utilized in previous studies to investigate the DS properties of SAC dendrites. These protocols were specifically designed to target a small patch of retina beneath a single SAC dendrite or induce circularly-symmetric shrinking or expanding motion to elicit consistent responses across the dendritic tree (Euler et al., 2002; Hausselt et al., 2007; Koren et al., 2017; Oesch & Taylor, 2010). Nonlinear integration in BCs triggered by these visual stimuli may unintentionally introduce a directional bias to BC activation - contradicting the fundamental assumptions of the space-time wiring hypothesis (Gaynes et al., 2022).

Due to the challenges associated with examining the relevance of the space-time wiring model in the context of single dendritic computations, we opted instead to focus on understanding the directional performance of SACs when the entire dendritic tree is stimulated. Our approach involved the integration of machine-learning algorithms with multicompartmental models of the bipolar–SAC network. In addition, we employed glutamate imaging techniques to record the waveforms of excitatory synaptic drive experienced by ON- and OFF-SACs during object motion. We reconstructed the underlying RF structure that supports different bipolar release shapes and subsequently combined them into a detailed model. By doing so, we sought to gain insights into how diverse RF characteristics of BCs contribute to the directional performance observed in SAC, while also comparing them to the maximal theoretical capabilities of space-time models achieved with synthetic presynaptic RFs. Our findings indicate that the receptive field diversity observed in bipolar cells generates response profiles that have the potential to contribute to SAC DS but to levels significantly below the maximum theoretical capabilities observed with synthetic RFs. Overall, by combining machine learning techniques, multicompartmental modeling, and glutamate imaging, our study provides insights into the complex relationship between BCs and SACs in the context of direction selectivity.

Results

SAC morphology imposes substantial signal attenuation but does not preclude inter-branch interactions

Calcium signals occurring in the terminal dendrites of starburst amacrine cells exhibit remarkable direction selectivity, both when visual stimulation is focused on the recorded dendrite and when stimuli engage the entire dendritic tree (Chen et al., 2016; Ding et al., 2016; Euler et al., 2002; Koren et al., 2017; Morrie & Feller, 2018; Poleg-Polsky & Diamond, 2016; Poleg-Polsky et al., 2018).

Because the precise nature of SAC function when a stimulus encompasses the entire cell remains unclear, our first goal was to determine whether whole-cell activation should be understood as a combination of multiple independent processing units or as the sequential recruitment of interacting dendrites. The distinctive dendritic morphology of SACs, characterized by long and narrow-diameter processes, was proposed to lead to the compartmentalization of individual or sister branches and significant signal attenuation to other dendrites (Masland, 2005; Poleg-Polsky et al., 2018). Yet theoretical studies suggested the presence of non-trivial dendritic crosstalk (Ding et al., 2016; Tukker et al., 2004).

To investigate whether the dendritic geometry and passive parameters support the hypothesis of independent computations in isolated dendritic compartments, we conducted a multicompartmental NEURON simulation to assess signal attenuation in reconstructed SACs. For our experiments, we selected a specific recording location on a terminal dendrite (Figure 1A, indicated in red) and administered 100 ms-long 10 pA current steps at various dendritic locations (Figure 1A). Consistent with prior studies, we observed a significant increase in signal attenuation as a function of the distance from the recording site (Stincic et al., 2016; Tukker et al., 2004; Vlasits et al., 2016). On average, approximately 40% of the injected signal was detectable at the recording location when the injection site was within the same branch, whereas stimulating other dendrites yielded negligible responses (Figure 1A, B).

Figure 1:

By utilizing the built-in NEURON ’impedance’ class to investigate signal propagation in relation to morphology, we observed a significant >20-fold attenuation on all branches distanced from the recording location by the soma. Notably, even sister dendrites sharing a common second-order bifurcation displayed a considerable >5-fold signal reduction (Figure 1A).

Subsequently, we examined whether dendritic isolation persists when multiple co-active inputs are present, which more accurately represents the activation pattern elicited by global motion. In mouse SACs, BCs innervate the proximal two-thirds of the dendritic tree (Ding et al., 2016; Vlasits et al., 2016). To replicate this arrangement, we distributed 200 synaptic inputs within a circle of 110 µm radius centered on the soma (Figure 1C). Each synapse exhibited a conductance of 0.1 nS and was stimulated by a 100 ms-long pulse.

Initially, we activated all the synapses that innervated dendrites sharing a common second-order branch. The peak voltage response following focal stimulation was measured at 10.9±3.3 mV (mean±SD, n=100, Figure 1C). Next, we performed a similar manipulation, now stimulating excitatory synapses on the remaining dendritic tree (Figure 1C). Notably, the depolarization recorded from the same branch during the stimulation of other dendrites was higher, reaching 13.2±1.9 mV (Figure 1C).

Hence, despite the substantial decrease in inter-dendrite signal propagation resulting from the cell’s morphology and passive characteristics, the model predicts significant voltage interaction throughout the cell. This observation can be attributed to the disparity in the number of active inputs between the two conditions. Specifically, the combined length of dendrites sharing a common second-order parent constitutes less than 10% of the entire dendritic tree’s length. This proportion also applies to the number of synapses involved in focal or global stimuli. Thus, although the efficacy of individual focal synapses is considerably more pronounced (Figure 1A, B), the larger number of inputs on other dendrites compensates for the signal attenuation caused by an individual synapse – leading to a comparable response to focal and global activation patterns (Figure 1C).

Development of a machine-learning model for enhancing direction selectivity in BC-SAC interactions

Our findings highlight the existence of long-distance interactions between synaptic inputs that extend beyond the boundaries of a single dendritic compartment. Unless explicitly mentioned otherwise, all subsequent simulations and experiments described below employed full-field moving bars. This choice is physiologically relevant as SACs exhibit a high degree of direction selectivity in response to such stimuli (Chen et al., 2016; Ding et al., 2016; Euler et al., 2002; Koren et al., 2017; Poleg-Polsky & Diamond, 2016; Poleg-Polsky et al., 2018). Additionally, full-field stimulation circumvents the nonlinear processing of edges and emerging objects in BCs (Gaynes et al., 2022; Strauss et al., 2022), which could potentially influence the interpretation of responses to stimuli appearing within the SAC receptive field (see discussion).

Previous studies suggested that differences in release profiles between proximal and distal BCs could contribute to directional tuning in isolated SAC dendrites, especially if proximal inputs exhibit more sustained and delayed dynamics (Fransen & Borghuis, 2017; Greene et al., 2016; Kim et al., 2014; Srivastava et al., 2022). However, do similar activation profiles promote direction selectivity during whole-cell integration? To investigate the impact of spatial differences in input kinetics on directional tuning, we extended the multicompartmental SAC model to include BCs, modeled as point neurons with center-surround receptive fields, each proving a single excitatory input to the postsynaptic SAC. Following the space-time model, we considered two distinct populations of BCs innervating the proximal and distal postsynaptic regions (Figure 2A, B). BCs within the same population shared identical RF characteristics, but the timing of their responses varied to account for the spatiotemporal progression of the visual stimulus over the simulated circuit.

Figure 2:

To investigate which bipolar dynamics promote postsynaptic direction selectivity, we employed a machine-learning approach based on evolutionary algorithms (EA, Figure 2C) (Ezra-Tsur et al., 2021). This methodology involved a population of models competing to achieve the training objective. In the initial generation or iteration of the algorithm, the parameters describing RF properties of the proximal and distal BC populations were randomly selected. We stimulated the network with horizontally moving bars and recorded calcium signals on terminal SAC dendrites near the horizontal axis. Direction selectivity index (DSI) was computed from the peak calcium values (Figure 2B, C). Initially, DSI levels were very low (mean (±SD) DSI = 7±3%, Figure S1), with variations among the models due to random initialization. The next generation of candidate solutions was formed by selecting the highest-scoring models to replace the below-average performers. To avoid overfitting, the models were trained with a range of stimulus velocities. Next, the RF properties of the BC population, as well as the global postsynaptic properties (axial resistance, leak, and voltage-gated calcium channel conductance), were subject to random mutations, and this cycle was repeated for 100 generations.

Figure 2D illustrates an exemplary solution achieved by the evolutionary algorithm, which displays transient distal inputs and slower, delayed proximal BCs. Whole-cell integration of these inputs in the SAC resulted in a significant directional preference (DSI across all velocities was 37±3%, Figure 2E).

Notably, the directional tuning profile was found to be dependent on the specified training goal of the model. In our standard configuration, we considered the mean DSI and the amplitude of calcium levels across the recorded dendrites as performance metrics. Including calcium amplitude was crucial in ensuring the robustness of the signal as it prevented the convergence on solutions that exhibited negligible absolute differences in calcium levels. Models solely rewarded based on DSI levels often evolved to minimize voltage and calcium levels, specifically for slow stimulation speeds (Figure S2). This outcome stands in contrast to physiological somatic signals and dendritic calcium responses, which exhibit consistent peak levels across the velocity spectrum (Figure S2) (Ding et al., 2016; Koren et al., 2017). Furthermore, the performance of this alternative approach, as measured by DSI, only marginally surpassed the standard configuration (Figure S2).

Presynaptic RF parameters influencing optimal directional performance in SACs

To assess the relative contribution of BC kinetics to the resulting direction selectivity, we modified the model such that all presynaptic cells shared the same receptive field formulation. It is worth noting that even with the requirement of identical properties in the bipolar population, the algorithm could explore a vast parameter space by adjusting both the RF properties of presynaptic cells and the postsynaptic integration parameters. However, despite this flexibility, the resulting simulations achieved only a fraction of the DS observed with spatially varying receptive fields (Figure 2E). This outcome, obtained by integrating realistic BC responses with morphologically detailed SACs, suggests that the space-time wiring model could have a major influence on directional computations in biological circuits.

Using a similar approach, we further investigated the contribution of individual RF parameters to SAC DS. In the modified model described above, all RF parameters in the proximal and distal presynaptic cells were clamped to the same values. Below, we constructed similarly structured models but allowed a single parameter to differ between the BC populations. Previous studies demonstrated that space-time wiring yields the highest DS levels when presynaptic units exhibit different activation lags or response transiency (Fransen & Borghuis, 2017; Kim et al., 2014; Matsumoto et al., 2019; Srivastava et al., 2022). Consistent with this, models in which proximal and distal BCs differed only in their response delay or decay times showed DSI levels close to those observed in the full model (29±2% and 35±3%, respectively, Figure 3).

Figure 3:

Notably, temporal differences in BC responses were not solely influenced by these parameters; we also found that the spatial extent of the RF could introduce an activation lag, which was inversely proportional to RF size: neurons with wider RFs can respond to a moving object sooner compared to their immediate neighbors with smaller RFs. In line with this observation, the optimal solutions in the full model evolved towards larger RF centers in the distal bipolar cells (42±15 µm vs. 32±9 µm, Figure 3C). The difference in receptive field size was more pronounced when center width was the only varying RF parameter (97±18 µm vs. 22±3 µm), resulting in a measurable lag (approximately 135 ms for stimulus velocity = 0.5 mm/s, Figure 3). The impact of spatial RF properties on directional tuning was not as pronounced as the temporal RF characteristics but managed to elevate the DS by about two-fold over the space-invariant configuration (DSI = 19±4%, p<10 vs. identical BCs, Figure 3C).

The space-time wiring model does not require dendritic isolation

Previous studies investigating the postsynaptic mechanisms underlying direction selectivity in SACs have highlighted the importance of electrotonic isolation of terminal dendrites in promoting directional signals (Ding et al., 2016; Koren et al., 2017; Poleg-Polsky et al., 2018). However, since the space-time wiring mechanism primarily operates on presynaptic circuits, it may be less reliant on postsynaptic compartmentalization. To test this hypothesis, we manipulated inter-dendritic signal propagation by modifying the axial and membrane resistance properties of the model in segments around the soma (Figure 4).

Figure 4:

As anticipated, increasing the axial resistance led to decreased somatic depolarization during simulated visual responses (Figure 4A). The disparity in voltage attenuation between the control configuration and the model with elevated axial resistance was particularly pronounced when stimulating a single SAC branch, confirming the near-perfect isolation of dendrites (Figure 4B). Conversely, reducing the barrier to signal propagation had the opposite effect (Figure 4A, B). However, despite these variations in perisomatic electrical properties, all models generated by the evolutionary algorithm achieved comparable DS levels (Figure 4C).

Collectively, these findings suggest that the effectiveness of the space-time wiring model relies primarily on BC kinetics, which are influenced by the spatiotemporal properties of their receptive field organization. In contrast, the specific details of postsynaptic integration appear to play a lesser role in this DS mechanism.

Diversity of glutamatergic inputs in murine SACs

In the previous sections, we examined the theoretical implementation of the space-time wiring hypothesis through detailed simulations of the bipolar–SAC circuit. To investigate the similarity between the receptive field properties of mouse BCs innervating SACs and the optimal receptive fields predicted by our models, we conducted glutamate imaging using a two-photon microscope on an ex-vivo retina preparation. To specifically target ON- and OFF-SACs, we virally expressed floxed iGluSnFR in ChAT-Cre mice (Figure 5). 3-6 weeks after vector introduction, we performed ex-vivo experiments using a wholemount retina preparation. Visual stimuli consisted of full-field flashes and bars moving alternatively either from left to right or in the opposite direction at five different velocities (Figure 5B) and oriented bars used to map RFs using the filtered back projection approach (Figure 5C).

Figure 5:

To analyze the diversity of BC signaling in response to motion, similarly responding pixels were grouped into regions of interest (ROIs) (Figure 5A, B). The time of response onset depends on intrinsic RF properties shared among BC subtypes and the position of the RF relative to the moving bar stimulus. To disentangle these components, we aligned the responses based on the mean 50% rise time observed during stimulation in both leftward and rightward directions. This alignment negated the dependency on the relative position of the RF, allowing us to focus solely on differences in response kinetics to motion.

Next, we analyzed common response motifs across recording regions. We identified functional release clusters using hierarchical clustering of the aligned and normalized dF/F waveforms evoked by a single motion velocity (0.5 mm/s). This procedure was performed separately for ON and OFF ROIs (n=334 and 135 regions, respectively, Figure 6A, C). We determined the optimal number of functional clusters from the curve describing the within-cluster variance (Figure S3). The analysis identified seven groups for the ON-SAC population (Figure 6A, S3) and six for the OFF-SACs (Figure 6C, S3) as providing the optimal separation (Franke et al., 2017; Gaynes et al., 2022; Matsumoto et al., 2021; Matsumoto et al., 2019; Rasmussen et al., 2020; Strauss et al., 2022). Given that our models highlighted the decay time as a critical indicator of potential contributions to DS, we sorted the functional clusters based on their decay times, with C1 representing the fastest decay (Figure 6B, D). As expected, cluster dynamics correlated with the stratification pattern in the IPL (Figure S4).

Figure 6:

Determination of RF characteristics from motion responses

An intriguing finding emerged from the analysis of the signals forming the glutamate release clusters. While most responses exhibited peak activity around the same time, ON-C2 and ON-C4 (and to a lesser extent ON-C1) displayed significantly shorter activation delays, and OFF-C5 tended to respond later to moving bars (Figure 6A, C). These distinct dynamics were observed consistently across different stimulation directions and trials (Figure S5).

The diversity in response timing within the release clusters could arise from either the temporal characteristics of the presynaptic cells or the spatial extent of their RFs, as discussed earlier. To investigate whether the temporal or the spatial RF properties contribute to the lag and to gain deeper insights into the RFs underlying the glutamatergic signals to SACs, we decided to train center-surround RF models to match the shapes of the recorded waveforms. We adapted model formulation described above for simulating BC dynamics, with the addition of a filter to account for iGluSnFR binding and unbinding dynamics, thereby simulating the sensor-mediated filtering of glutamate signals drive (Armbruster et al., 2020; Hain & Moser, 2023; Srivastava et al., 2022). Separate evolutionary algorithm models were trained for each functional cluster, utilizing five stimuli with different velocities to drive the model. The objective of the training was to replicate the shape of the recorded fluorescent signals during the presentation of moving bars at corresponding velocities (Figure 7A).

Figure 7:

Using this approach, we were able to faithfully reproduce the empirical glutamate dynamics observed in the training dataset (Figure 7A, S6). As an initial validation of the algorithm’s predictive power, we compared glutamate signals to the response of the RF models measured during the presentation of full-field flashes (Figure 7B). Encouragingly, even though these stimuli were not part of the training set, the simulated output reliably replicated experimental waveforms, despite substantially different dynamics in many clusters (Figure 7B, S7).

Furthermore, by leveraging the more interpretable stimulus-response dynamics of the full-field flash paradigm, we found a strong agreement between response delay (lag) and the rise time of the flash response, which closely correlated with the rise time kinetics of the center RF component (Figure 7C). Notably, ON-C4, which exhibited an earlier activation during motion compared to other clusters, had the longest lag as predicted by the model (181±25 ms) and confirmed by the flash response (158±37 ms); about 78 ms longer than the lag recorded for the quickest to respond cluster (80±21 ms, ON-C1). The most plausible explanation for cells with a prolonged activation lag responding earlier to a moving stimulus is a larger RF size. Consistently, the RF model estimated a half-width of 379±13 µm for ON-C4, which was validated through RF mapping (286±55 µm, Figure 7C, S8). A somewhat smaller but still substantial center RF size was predicted and experimentally measured for the second early-response cluster (ON-C2), whereas the remaining functional clusters were found to have narrow RFs, consistent with the classical RF description for BCs (Figure 7C, S6) (Euler et al., 2014; Franke et al., 2017; Kuo et al., 2016; Schwartz et al., 2012; Strauss et al., 2022; Turner et al., 2018; Wienbar & Schwartz, 2018).

Realistic BC dynamics suggest a modest effect of the space-time wiring on directional tuning in ON- and OFF-SACs

After establishing the diversity of excitatory dynamics during visual motion, we proceeded to examine the extent to which their combination could enhance direction selectivity in SACs. To achieve this, we replaced the synthetic bipolar description in the bipolar-SAC evolutionary algorithm model with fits to experimentally determined glutamatergic waveforms. These fits provided a deconvolved representation of the excitatory drive, which is more likely to approximate the actual presynaptic signals reaching the SACs. We note that the difference between the recorded and deconvolved waveforms was minimal, and utilizing the experimentally recorded shapes did not significantly impact the subsequent results (data not shown).

We first employed a simplified model where the experimentally determined clusters drove either the proximal or distal halves of the excitatory inputs to the SAC model. Figure 8 illustrates the peak DS index achieved with all possible proximal and distal functional cluster placement combinations. As anticipated, the optimal DS was attained when proximal clusters with sustained waveforms were combined with transient distal ones. For ON-SACs, the highest DSI = 0.22 ± 0.01 was observed with inputs from ON-C1 and ON-C6. In the case of OFF-SACs, the inputs performed slightly better, with the highest DSI levels (0.28 ± 0.01) associated with OFF-C1 and OFF-C5. However, as these DSI values did not exhibit the levels of directional tuning achieved with synthetic BC receptive fields (Figure 2), we asked whether the model, restricted to contain only two input populations at a time, is not rich enough to take full advantage of the heterogeneity of the experimentally measured signals. If this is the case, allowing SACs to integrate over more presynaptic clusters should improve directional outcomes.

Figure 8:

Inspired by a previous model (Fransen & Borghuis, 2017), we distributed presynaptic filters over concentric annuli, preserving input symmetry along the soma-dendritic tip axis as postulated by the space-time wiring hypothesis. Specifically, we grouped the inputs into ten µm-wide bins based on their distance from the postsynaptic soma. Each bin was assigned a single functional cluster RF description drawn at random. The search algorithm was designed to determine the spatial distribution of clusters that produced the most robust directional tuning across a range of stimulus velocities (Figure 9).

Figure 9:

The proposed solutions converged on similar outcomes, irrespective of the initial distribution. We observed that optimal DS was achieved in both ON and OFF circuits when two clusters dominated the proximal and distal innervation regions (Figure 9A, D). Consequently, the levels of DS obtained in these simulations did not significantly differ from the findings of the simple proximal-distal parameter exploration presented earlier (ON-SACs: DSI = 0.22 ± 0.008, OFF-SACs: DSI = 0.28 ± 0.008, Figure 9C, F). Based on these results, we can conclude that while excitatory signals, even when optimally combined, contribute to DS to some extent, they can only account for a fraction of the DS level observed empirically.

Discussion

In this study, through a combination of experimental approaches and computational modeling, we have explored the theory and the evidence for bipolar organization supporting the space-time wiring model. The original space-time wiring hypothesis was conceptualized for individual SAC dendrites, and the difference in BC response shape was considered through the lens of the dynamics of responses to static stimuli (Kim et al., 2014). Although these circumstances were thought to simplify the analysis of dendritic computations, it is now clear that signal preprocessing in the bipolar population introduces unforeseen challenges to this approach: Our prior work demonstrated a low correlation between the dynamics of BC responses to motion and static stimuli (Gaynes et al., 2022). Here, we replicate this result on a new dataset (Figures 7, S6), reinforcing the notion that signal kinetics, crucial to the space-time wiring model, are not a fixed property of the cell but are instead dependent on the stimulus. For example, we demonstrate that presynaptic signals can exhibit prolonged flash-lags yet generate early responses to motion, a phenomenon explained by their spatial RF structure (Figure 7). In addition, it is vital to consider the implications of using probes that activate only a small portion of the SAC dendritic tree in isolation or all dendrites through the same stimulation pattern. Such probe designs are likely to introduce nonlinear signal processing in BCs and, importantly, result in asymmetric bipolar responses for different stimulation directions (Gaynes et al., 2022; Strauss et al., 2022). This poses a concern to the space-time model as it assumes consistent presynaptic waveforms independent of the stimulus direction (Kim et al., 2014).

To overcome these issues, it becomes crucial to devise probes that effectively bypass both limitations. One approach is to design stimuli that either appear or begin to move outside the receptive field of the BCs innervating the postsynaptic cell. By doing so, we can ensure that the probes do not trigger asymmetric responses and maintain the consistency of presynaptic waveforms, as required by the space-time model. A simple stimulus that fulfills this criterion is a continuously moving object that appears and disappears at a considerable distance (several hundred microns) from the cells of interest (Gaynes et al., 2022; Wu et al., 2023). By employing this stimulus design, we can effectively mitigate the nonlinear signal processing in BCs and preserve the foundational assumption of consistent presynaptic waveforms in the space-time model. This approach allows for more accurate investigations into the dynamics of BC responses and their contribution to the directional tuning in SACs and DSCGs.

While the use of full-field stimuli restricts the ability to investigate the responses of individual branches with high granularity, it provides an opportunity to explore a broader question: the potential impact of interactions among multiple dendrites on direction selectivity. Dendritic compartmentalization in SACs is promoted by morphological factors (Tukker et al., 2004; Wu et al., 2023) and by active processes, such as the presence of potassium channels (Ozaita et al., 2004). Additionally, in the mouse, SAC-SAC inhibition predominantly occurs in the perisomatic region, where it modulates communication between branches and maintains the output synapses of SACs within the optimal range for effective DS processing (Ding et al., 2016; Poleg-Polsky et al., 2018). Similar enhancement of DS signaling is mediated by metabotropic glutamate receptors. By limiting calcium influx, they effectively reduce the spread of electrical signals within the dendritic arbor, enhancing the electrotonic isolation of individual processing units (Koren et al., 2017). Interestingly, our results demonstrate that in contrast to these postsynaptic mechanisms, compartmentalization has little effect on directional tuning produced by the space-time wiring model (Figure 4). This unexpected outcome highlights the distinction between pre and postsynaptic computations within the bipolar-SAC circuit.

Our findings establish an upper limit on the contribution of the space-time wiring mechanism to SAC direction selectivity. While we observed a slightly higher direction selectivity in the OFF population, consistent with a previous study (Fransen & Borghuis, 2017), our investigation reveals that the waveforms produced by presynaptic cells in response to moving objects do not fully exploit the potential for directional tuning that can be achieved with synthetic BCs, as depicted in Figure 9. These results strongly indicate that additional postsynaptic mechanisms are necessary to account for the observed direction selectivity capabilities in a manner consistent with empirical evidence.

In this study, we employed evolutionary algorithms, a powerful machine-learning method, to investigate the mechanisms of DS computations in SACs. Evolutionary algorithms have broad applicability and have been successfully utilized in examining the mechanisms of directional computations of SACs before (Ezra-Tsur et al., 2021). While they place minimal constraints on the model type or examined features, it is important to note that evolutionary algorithms are most effective when dealing with models with a relatively constrained feature space. While they may be susceptible to getting trapped in local minima, our largest model, which consisted of 23 free parameters (9 x 2 presynaptic + 5 postsynaptic, as described in the methods), consistently evolved towards similar solutions regardless of the starting parameters. In cases where evolutionary algorithms do not reliably converge on an optimal solution, alternative machine-learning approaches can be employed. For instance, a recent study utilized a neural network system that combined Hassenstein-Reichardt correlators with biologically inspired receptive fields to describe solutions for computations in collision avoidance neurons of flies (Zhou et al., 2022).

In this work, we present a novel application of model fitting to map receptive fields based on responses to moving stimuli. Previous studies have used stimuli moving in different directions to estimate the spatial location of the RF (Gaynes et al., 2022; Wienbar & Schwartz, 2018). We expand upon this approach by incorporating multiple probe velocities. This extension allows efficient RF determination from a brief visual stimulation (approximately 2-3 minutes). However, it is essential to acknowledge that our fits to the data were based on the assumption of a simple center-surround RF structure and may not fully capture the complexity of more diverse RF types or subcellular release specialization present in some BCs. In particular, our study did not specifically investigate the measurement of directional release from individual axonal terminals, as previously demonstrated in type 2 and 7 bipolar cells that innervate OFF- and ON-SACs (Matsumoto et al., 2021). In addition, it is presently unknown whether the contacts made by DS boutons on SACs are aligned with the soma-dendritic tip axis to enhance DS in the outward direction (Srivastava et al., 2022). These aspects represent avenues for future research and could provide further insights into the mechanisms underlying direction selectivity in retinal circuits.

We were intrigued by the discovery of presynaptic cells exhibiting broad receptive fields that spanned hundreds of microns. Our findings bear a resemblance to a recent study that described elongated receptive fields in type 5A bipolar cells (Hanson et al., 2023), challenge the classical understanding of bipolar RF structure, which suggests that the extent of the center component is limited to the size of their dendritic fields, typically below 50 µm in the mouse (Euler et al., 2014; Franke et al., 2017; Gaynes et al., 2022; Kuo et al., 2016; Schwartz et al., 2012; Strauss et al., 2022; Turner et al., 2018; Wienbar & Schwartz, 2018), but see (Asari & Meister, 2012) for an example of large bipolar RFs in the salamander. One plausible mechanistic explanation for the observed extensive RFs is the presence of gap-junction coupling among neighboring BCs (Kuo et al., 2016), which could be part of the adaptive changes in the bipolar–SAC circuits to repeated stimulation (Ankri et al., 2020; Vlasits et al., 2014).

Overall, by combining experimental methods and computational modeling, our study contributes to a deeper comprehension of the bipolar organization supporting the space-time wiring model. We address the limitations of solely focusing on static stimuli and shed light on the dynamic nature of signal processing in bipolar cells. This knowledge is crucial for elucidating the intricate mechanisms underlying visual computations in the DS circuit in the retina.

Methods

NEURON simulation

Multicompartmental simulations were performed using NEURON 8.2 (www.neuron.yale.edu/neuron) on four reconstructed SAC morphologies (https://neuromorpho.org/neuron_info.jsp?neuron_name=185exported, https://neuromorpho.org/neuron_info.jsp?neuron_name=cell8_wt_traces, https://neuromorpho.org/neuron_info.jsp?neuron_name=cell2_wt_traces, https://neuromorpho.org/neuron_info.jsp?neuron_name=cell1_wt_traces). The number of segments varied from 495 to 879. The diameter of branches further than 30 µm from the soma was set to 200 nm. The initial global passive parameters were as follows: passive conductance = 4e S/cm, membrane capacitance = 1 µF/cm, reversal potential = -60 mV and axial resistance = 150 Ωcm. N-type calcium conductance model, adapted without modifications from (Benison et al., 2001), was distributed throughout the entire dendritic tree. In order to calculate the internal calcium levels, a calcium diffusion mechanism was incorporated into all dendrites. This mechanism had a time constant of 50 ms, and the resting calcium levels were set to 100 nM.

The SAC received innervation from 200 BC inputs, which were randomly distributed across proximal dendrites (within a distance of <110 µm from the soma). Each synapse was assumed to be associated with a single BC, and the RF center of the BC was aligned with the postsynaptic position. The BC RFs consisted of Gaussian-shaped components for the center and surround regions.

Stimuli were presented within a 1 mm-wide arena, matching the experimentally presented stimuli (see below). To determine RF activation, the spatial overlap between the Gaussian functions describing the center and surround RF components and the shape of the stimulus was computed separately for each time step (Δt = 1 ms).

Area_t 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_t was the sum of the area of the Gaussian function describing the center or surround RF component that spatially overlapped with the stimulus:

Where width is the FWHM of the corresponding RF component. Subsequently, RF center and surround responses at time t were determined from the following equations:

The integration of center and surround RFs was calculated in a conductance-based model of synaptic integration, determined from the excitatory and inhibitory driving forces:

Where R_in is the input resistance. Because the magnitude of RF activation depended on its size, R_in was set to be proportional to the spatial extent of the center.

RF activation in each BC was temporally adjusted based on stimulus velocity and the spatial position:

The following constraints were imposed: RF amplitude, surround strength and reversal were permitted to vary between 0 and 1; rise / decay times and RF widths were positive. Overall, RF_full was expressed in unitless units in the [0,1] interval, with values closer to one signifying strong activation.

Training of evolutionary models

During the training of evolutionary models aimed at understanding the contribution of BC dynamics to SAC DS (as shown in Figures 2-4), the bipolar population was equally divided into proximally and distally innervating cells according to their proximity to SAC soma.

For each generation, a total of 16 bipolar-SAC models were executed simultaneously on a Dell Precision 5820 workstation, with one model per processor thread. In the first generation, random values were assigned to seed each network model. A separate RF description was instantiated with random values within the specified limits for the two BC populations.

Likewise, random initial values were chosen for passive conductance (ranging from 1e to 1e S/cm), axial resistance (constrained between 50 and 300 Ωcm), N-type calcium conductance and voltage offset (limits = ±30 mV), as well as synaptic conductance (limits = 0.01 – 1 nS). All models were presented with the same stimuli, consisting of bars moving from left to right and from right to left at five different speeds: 0.25, 0.5, 1, 2 and 4 mm/s. The height of the bar was 1 mm, its duration was 2 seconds. The intensity of the stimuli was set to 1 and the background was set to zero (AU).

Calcium signals were recorded from SAC sites located on the distal dendrites, specifically those within proximity of less than 30 µm from the horizontal axis. When a single branch was stimulated, only sites on that dendrites were included in the analysis. For each stimulation speed (comprising two different stimulation directions), we computed the DS index of each site as follows:

Where R_Outward and R_Inward are the peak calcium levels recorded for the outward and inward directions from the perspective of the dendrite. To minimize extreme calcium signals, we computed the directional metric:

Where Ca_opt was the optimal calcium level in the outward direction, set to 500 nM.

Subsequently, the models were ranked based on the average directional metric calculated over the five stimulation speeds. The two best-performing models were retained without any changes, while the parameters of the rest of the models were randomly selected to match either the best or second-best-performing models and then mutated as follows: For each parameter describing the presynaptic RF and postsynaptic SAC properties, a Gaussian distribution with a mean of 1 and a standard deviation of 5% 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. These modified (mutated) models constituted the next generation of candidate solutions. Typically, we evolved the model over 100 generations, as empirical evidence has shown that increasing the number of generations beyond this point does not yield significantly improved results.

RF estimation from recorded glutamate waveforms

To estimate the RF that could mediate the recorded glutamate waveforms presented in Figures 7 and S5, we employed an evolutionary algorithm for training RF models to match the shape of experimentally observed glutamate release clusters. Each cluster was considered separately during the analysis.

The model contained a single RF, and its parameters were described using the same parameter set as mentioned earlier. It is important to note that the fluorescent signal produced by iGluSnFR represents a temporally filtered version of the original glutamatergic drive (Armbruster et al., 2020; Hain & Moser, 2023; Srivastava et al., 2022). In order to mimic the iGluSnFR waveforms in the models, we applied a filtering process to the output of the simulation. On each time step t0, we convolved the value of the simulated response with a wavelet described as the difference between two exponential functions:

Where F_rise = 10 ms, F_decay =50 ms are the rise and decay times of the filter, respectively. The resulting vectors starting at time t0 and lasting till the end of simulation duration, were combined to produce the filtered version of the full response.

Following the initial random instantiation, simulated RF responses were calculated for the five different speeds. The fitness of each model was then evaluated based on the mean square error between the simulated waveforms and the experimentally recorded waveforms. The next generation was formed by introducing potentially mutated offspring of models that exhibited the lowest errors. These mutations were performed according to the methodology described above.

Determination of optimal DS with experimentally recorded glutamatergic waveforms

To investigate the impact of the spatial dependence of the presynaptic innervation waveform on postsynaptic DS, we conducted an initial analysis using a modified version of the bipolar–SAC model. In this derivative model, we systematically replaced the RF descriptions of proximal and distal BCs with fits to experimental data. The focus of the evolution process was on altering postsynaptic parameters to optimize synaptic integration toward the largest directional performance. For the data shown in Figure 9, we employed a different training approach. We divided the input synapses into annuli with a width of 10 µm, centered around the soma. Initially, each annulus was randomly assigned one of the RF fits corresponding to functional clusters.

As the model evolved, its goal was to achieve the best possible directional performance based on the postsynaptic parameters. Additionally, there was a 5% probability for each annulus to replace the assigned functional cluster identity. This allowed for evolution in the spatial configuration of the functional clusters, enabling us to examine what spatial distribution of experimentally recorded waveforms produces the optimal DS.

Virus expression and imaging procedures

All animal procedures were conducted in accordance with US National Institutes of Health guidelines, as approved by the University of Colorado Institutional Animal Care and Use Committee (IACUC). Mice were housed in a 12 light /12 dark cycle at room temperature (∼22°C), 40-60% relative humidity. For intravitreal virus injections, 8–12-week-old ChAT-Cre transgenic mice (Jax strain 031661, https://www.jax.org/strain/031661) were anesthetized with isoflurane; ophthalmic proparacaine and phenylephrine were applied for pupil dilation and analgesia. A small incision at the border between the sclera and the cornea was made with a 30 gauge needle. 1 µL of AAV9.hsyn.FLEX.iGluSnFR.WPRE.SV40 (a gift from Loren Looger, Addgene plasmid # 98931; http://n2t.net/addgene: 98931; RRID:Addgene_98931, 10 vg/mL in water) solution was injected with a blunt tip (30 gauge) modified Hamilton syringe (www.borghuisinstruments.com). Experiments on retinas from all animal groups were performed 2-6 weeks following virus injection on 11–17-week-old animals (4 males and 4 females).

Mice were not dark-adapted to reduce rod-pathway activation. Two hours after enucleation, retina sections were whole mounted on a platinum harp with their photoreceptors facing down, suspended ∼1 mm above the glass bottom of the recording chamber. The retina was kept ∼32°C and continuously perfused with Ames media (Sigma-Aldrich, www.sigmaaldrich.com) equilibrated with 95% O₂ / 5% CO₂.

Visual stimulation

Light stimuli were generated in Igor Pro 8 (Wavemetrics, www.wavemetrics.com) running in Windows 10 and displayed with a DPL projector (Texas Instruments, www.ti.com, model 4710EVM-G2) connected as a second monitor. Only the blue projector LED was used, and its light was further filtered with a 450 nm lowpass filter (Thorlabs, www.thorlabs.com, FEL0450). Light from the visual stimulus was focused by the condenser to illuminate the tissue at the focal plane of the photoreceptors (resolution = 10 µm/pixel, background light intensity = 30,000-60,000 R* rod^-1). Both vertical and horizontal light stimulus positions were checked and centered daily before the start of the experiments. The following light stimulus patterns were used: static flash covering the entire display (1000x1000 µm) presented for 2-4 seconds. A 1 mm-long bar moving either to the left or the right directions (speeds = 0.25, 0.5, 1, 2, 4 mm/s; dwell time over each pixel = 2 s). These stimuli were repeated three times. Typically, stimulus contrast was set to 60% Michelson contrast. To record glutamate signals to OFF-SACs, we reversed the intensity of the moving bar stimulus and the background. In some experiments, we masked different portions of the display to remain at background light levels (Gaynes et al., 2022). To map the spatial RFs, we flashed 10x1000 µm bars for 200 ms every 400 ms. The bars were presented over five evenly spaced (36°) orientations in a pseudorandom sequence over 32 spatial positions in every orientation to densely cover ∼300 µm of visual space centered on the imaged region.

Recording procedures

Glutamate imaging was performed with Throlabs Bergamo galvo-galvo two-photon microscope using the Thorimage 4.1 acquisition software (Throlabs). A pulsed laser light (920 nm, ∼1 µW output at the objective; Chameleon Ultra II, Coherent, www.coherent.com) was used for two-photon excitation projected from an Olympus 20X (1 NA) objective. A descanned (confocal) photomultiplier tube (PMT) was used to acquire fluorescence between 500 and 550 nm. The confocal pinhole (diameter = 1 mm) largely prevented stimulus light (focused on a different focal plane), from reaching the PMT, allowing us to present the visual stimulus during two-photon imaging. A photodiode mounted under the condenser sampled transmitted laser light to generate a reference image of the tissue. Fluorescence signals were collected in a rapid bidirectional frame scan mode (128x64 pixels; ∼50 Hz). The line spacing on the vertical axis was doubled to produce a rectangular imaging window (∼164x164 µm in size; the corresponding pixel size was 1.28 µm). To reduce shot noise, images were subsampled by averaging 2x2 neighboring pixels and filtered by a 20 Hz low pass filter offline. Horizontal and vertical image drifts were corrected online using a reference z-stack acquired before time-series recordings.

Labeled cells in the Chat-Cre/tdTomato line (Jax strain 007909, https://www.jax.org/strain/007909) were targeted for whole-cell recordings using 4-8 MΩ pipettes filled with 90 mM CsCH₃S0₄, 20 mM TEA-Cl, 10 mM HEPES, 10 mM EGTA, 10 mM phosphocreatine disodium salt hydrate, 5 mM QX-314, 4 mM Mg-ATP and 0.4 mM Na-GTP. Recordings were obtained in voltage-clamp configuration with a Double IPA amplifier (Sutter, www.sutter.com), low pass filtered at 2 kHz using a custom acquisition software written in Igor Pro.

Analysis

The fluorescence signals were averaged across multiple presentations of the visual protocol using Igor Pro 8. Specifically, pixels with dF/F values greater than 20% were selected for subsequent clustering analysis, which involved two steps.

In the first step, the extent of contiguous ROIs with similar response kinetics was determined. This was achieved by first averaging the fluorescence signals across repeated presentations of the stimuli. Then, the response shape was measured within a 1-second window centered around the time when the bar motion reached the center of the display during the ten motion trials, which included five speeds and two directions. To perform hierarchical clustering of the responses, the built-in Igor Pro ’FPClustering’ function was utilized. The number of resulting ROIs was adjusted until no ROI spanned more than 10 µm in size. The ROIs were manually curated, and any ROIs with pixel variability exceeding a coefficient of variation threshold of 1 were removed from further analysis.

In the second step, we combined individual ROIs from multiple scan fields and mice into functional clusters. This process was separate for ON- and OFF-SACs. For each ROI, we computed the mean motion response shape by aligning the motion responses to both directions of stimulation by the half-maximum rise time. Secondary hierarchical clustering was performed in R v3.6 (R foundation for statistical computing). Motion responses to 0.5 mm/s stimuli were first converted to data frames in R. The traces were smoothed with a combination of a Butterworth (filter order = 1, critical frequency = 0.1) and rolling average filters, baseline subtracted and normalized to the peak amplitude. For the ON responses, we discarded time points with variance less than 0.5 as this aided with convergence to biologically plausible clusters. To quantify the number of functional clusters, we used within cluster variance (Warren Liao, 2005). Within cluster variance is computed by finding the sum of the squared differences of each observation in a proposed cluster and the mean of that cluster. We used automated assessment of the elbow plots (within-cluster variance as a function of the number of clusters) to identify optimal cluster numbers for a grid search of distance computations and clustering algorithms available in base R and the dynamic time warping library. We then used expert assessment to identify the most plausible clustering from this small set of options (Figure S3). Empirically, we found that the ward.D2 hierarchical clustering method, available in the base R ‘hclust’ function, provides the clearest cluster separation. We used the Euclidean and maximum distances for the ON- and OFF-SAC datasets, respectively.

Tests of statistical significance were performed in Igor Pro using built-in functions.

RF mapping was determined from a two-dimensional Gaussian fit to the responses to the oriented bars, as described previously (Johnston et al., 2014; Poleg-Polsky et al., 2018). The FWHM of the RF was calculated along the axis of motion from the spatial time constant generated by the fit (σ_x) as follows:

The transiency index (TI) was calculated as the ratio between the peak and the mean of the response within the stimulation window. TI=1 indicates a sharp and transient response, TI close to zero is produced by sustained plateaus.

IPL depth was measured from the transmitted light channel extracted from the z-stack taken of the entire width of the retina that accompanied all functional recordings. The curvature of the retina was corrected by measuring the height of the inner limiting membrane at the four corners and the center of the image stack and fitting a curved plane that crossed these 5 points.

Code availability

The code for the visual stimulation and simulations is available on https://github.com/PolegPolskyLab/ This work was supported by NIH grant (R01 EY030841) to Alon Poleg-Polsky.

Author contributions

John Gaynes: Investigation, Writing-Review &Editing

Samuel Budoff: Formal analysis, Software, Writing-Review &Editing

Michael Grybko: Investigation, Writing-Review &Editing

Alon Poleg-Polsky: Conceptualization, Methodology, Software, Formal analysis, Resources, Data Curation, Writing-Original draft, Supervision, Project administration, and Funding acquisition

Competing interests

Authors declare no competing interests