On the potential role of lateral connectivity in retinal anticipation
OOn the potential role of lateral connectivity in retinalanticipation
Selma Souihel and Bruno Cessac Université Côte d’Azur, Inria, Biovision team and Neuromod Institute, France.
Abstract
We analyse the potential effects of lateral connectivity (amacrine cells and gap junctions) onmotion anticipation in the retina. Our main result is that lateral connectivity can - under con-ditions analysed in the paper - trigger a wave of activity enhancing the anticipation mechanismprovided by local gain control [8, 17]. We illustrate these predictions by two examples studiedin the experimental literature: differential motion sensitive cells [1] and direction sensitive cellswhere direction sensitivity is inherited from asymmetry in gap junctions connectivity [73]. Wefinally present reconstructions of retinal responses to 2D visual inputs to assess the ability of ourmodel to anticipate motion in the case of three different 2D stimuli.
Keywords—
Retina, motion anticipation, lateral connectivity, 2D
Our visual system has to constantly handle moving objects. Static images do not exist for it, as the en-vironment, our body, our head, our eyes are constantly moving. A "computational", contemporary view,assimilates the retina to an "encoder", converting the light photons coming from a visual scene into spiketrains sent - via the axons of Ganglion cells (GCells) that constitute the optic nerve - to the thalamus, and thento the visual cortex acting as a "decoder". In this view, comparing the size and the number of neurons in theretina - about million of GCells (humans) - to the size, structure, and number of neurons in the visual cortex(around million per hemisphere in the human visual cortex [19]) the "encoder" has to be quite smart toefficiently compress the visual information coming from a world made of moving objects. Although it haslong been thought that the retina was not more than a simple camera, there are more and more evidencesthat the retina is "smarter than neuroscientists believed" [35]. It is indeed able to perform complex tasks andgeneral motion features extractions such as approaching motion, differential motion, motion anticipation,allowing the visual cortex to process visual stimuli with more efficiency.The process leading from the photons reception in the retina to the cortical response takes about − milliseconds. Most of this delay is due to photo-transduction. Though this might look fast, it is actuallytoo slow. A tennis ball moving at m/s - km/h (the maximum measured speed is about km/h)covers between . and m during this time, so, without a mechanism compensating this delay it wouldn’tbe possible to play tennis (not to speak of survival, a necessary condition for a species to reach the level whereplaying tennis becomes possible). The visual system is indeed able to extrapolate the trajectory of a movingobject to perceive it at its actual location. This corresponds to anticipation mechanisms taking place in thevisual cortex and in the retina, with different modalities [77, 49, 4, 50].In the early visual cortex an object moving across the visual field triggers a wave of activity ahead ofmotion, thanks to the cortical lateral connectivity [7, 70, 39]. Jancke et al. [39] first demonstrated the existence a r X i v : . [ q - b i o . Q M ] S e p f anticipatory mechanisms in the cat primary visual cortex. They recorded cells in the central visual fieldof area 17 (corresponding to the primary visual cortex) of anaesthetized cats, responding to small squares oflight, either flashed or moving in different directions, and with different speeds. When presented with themoving stimulus, cells show a reduction of neural latencies, as compared to the flashed stimulus. Subra-maniyan et al. [70] have reported the existence of similar anticipatory effects in the macaque primary visualcortex, showing that a moving bar is processed faster than a flashed bar. They give two possible explanationsto this phenomenon : either a shift in the cells receptive fields induced by motion, or a faster propagation ofmotion signals as compared to the flash signal.In the retina, anticipation takes a different form. One observes a peak in the firing rate response of GCellsto a moving object, occurring before the peak response to the same object when flashed. This effect can be ex-plained by purely local mechanisms, at individual cells level [8, 17]. To our best knowledge, collective effectssimilar to the cortical ones - that is, a rise in the cell’s activity before the object enters in its receptive field dueto a wave of activity ahead of the moving object - have not been reported yet.In a classical, Hubel-Wiezel-Barlow [37, 5, 53] view of vision, each retinal ganglion cell carries a flow ofinformation with an efficient coding strategy maximizing the available channel capacity by minimizing theredundancy between GCells. From this point of view, the most efficient coding is provided when GCellsare independent encoders (parallel streaming identified by a "I" in Fig. 1). In this setting one can propose asimple and satisfactory mechanism explaining anticipation in the retina, based on gain control at the level ofBipolar cells (BCells) and GCells (label "II" in 1) [8, 17].Yet, some GCells are connected. Either directly, by electric synapses-gap junctions (pathway IV in Fig.1), or indirectly, via specific Amacrine cells (ACells, pathway III in Fig. 1). It is known that these pathwaysare involved in motion processing by the retina. AII ACells play a fundamental role in the interaction be-tween the ON and OFF cone pathway [47]. There are GCells able to detect the differential motion of an objectonto a moving background [1], thanks to ACells lateral connectivity. Some GCells are direction sensitive be-cause they are connected via a specific, asymmetric, gap junctions connectivity [73]. Could lateral connectivityplay a role in motion anticipation, inducing a wave of activity ahead of the motion, similar to the cortical anticipationmechanism ?
While some studies hypothesize that local gain control mechanisms can be explained by theprevalence of inhibition in the retinal connectome [40], the mechanistic aspects of the role of lateral connec-tivity on motion anticipation has not, to the best of our knowledge, been addressed yet on either experimentalor computational grounds.In this paper, we address this question from a modeller, computational neuroscientist, point of view.We propose here a simplified description of the pathways I, II, III, IV of Fig. 1, grounded on biology, butnot sticking at it, to numerically study the potential effects of gain control combined with lateral connectiv-ity - gap junctions or ACells - on motion anticipation. The goal here is not to be biologically realistic but,instead, to propose from biological observations potential mechanisms enhancing the retina’s capacity to an-ticipate motion and compensate the delay introduced by photo-transduction and feed-forward processing inthe cortical response. We want the mechanisms to be as generic as possible, so that the detailed biologicalimplementation is not essential. This has the advantage of making the model more prone to mathematicalanalysis.The first contribution of our work lies in the development of a model of retinal anticipation where GCellshave gain control, orientation selectivity and are laterally connected. It is based on a model introduced byChen et al. in [17] - itself based on [8] - reproducing several motion processing features: anticipation, alertresponse to motion onset and motion reversal. The original model handles one dimensional motions andits cells are not laterally connected (only pathways I and II were considered). The extension proposed herefeatures cells with oriented receptive field, although our numerical simulations do not consider this case (seediscussion). Lateral connectivity is based on biophysical modelling and existing literature [71, 25, 1, 36, 73].In this framework, we study different types of motion. We start with a bar moving with constant speedand study the effect of contrast, bar size, and speed on anticipation, generalizing previous studies by Berry
Synthetic view of the retina model.
A stimulus is perceived by the retina, triggeringdifferent pathways.
Pathway I (blue) corresponds to a feed-forward response where, from topto bottom: The stimulus is first convolved with a spatio-temporal receptive field that mimics theOuter Plexiform Layer (OPL) ("Bipolar receptive field response"). This response is rectified bylow voltage threshold (blue squares). Bipolar cells responses are then pooled (blue circles withblue arrows) and input Ganglion cells. The firing rate response of a Ganglion cell is a sigmoidalfunction of the voltage (blue square). Gain control can be applied at the Bipolar and Ganglion cellslevel (pink circles) triggering anticipation. This corresponds to the label
II (pink) in the figure.Lateral connectivity is featured by pathway III (brown) through ACells, and pathway IV (green) through gap-junctions at the level of GCells. 3 t al [8] and Chen et al [17]. We then extend the analysis to two dimensional motions, investigating e.g.angular motion and curved trajectories. Far from making an exhaustive study of anticipation in complexstimuli, the goal here is to calibrate anticipation, without lateral connectivity, so as to compare the effectwhen connectivity is switched on.The second contribution emphasizes a potential role of lateral connectivity (gap junctions and ACells) onanticipation. For this, we first make a general mathematical analysis concluding that lateral connectivity caninduce a wave triggered by the stimulus which, under specific conditions can improve anticipation. The effectdepends on the connectivity graph and is non linearly tuned by gain control. In the case of gap junctions,the wave propagation depends whether connectivity is symmetric (the standard case) or asymmetric, asproposed by Trenholm et al. in [73] for a specific type of direction sensitive GCells. In the case of ACells, theconnectivity graph is involved in the spectrum of a propagation operator controlling the time evolution of thenetwork response to a moving stimulus. We instantiate this general analysis by studying differential motionsensitive cells [1] with two types of connectivity: nearest neighbours, and a random connectivity, inspiredfrom biology [71], where only numerical results are shown. In general, the anticipation effect depends onthe connectivity graph structure and the intensity of coupling between cells as well as on the respectivecharacteristic times of response of cells, in a way that we analyse mathematically and illustrate numerically.We actually observe two forms of anticipation. The first one, discussed in the beginning of this intro-duction and already observed in [8, 17], is a shift in the peak of a retinal Gcell response, occurring before theobject reaches the center of its receptive field. In our case, lateral connectivity can enhance the shift improvingthe mere effect of gain control. The second anticipation effect we observe is a raise in GCells activity beforethe bar reaches the receptive field of the cell, similarly to what is observed in the cortex [7]. To the best of ourknowledge, this effect has not been studied in the retina and constitutes therefore a prediction of our model.The paper is organized as follows. Section 2 introduces the model of retinal organization and cells typesdynamics, ending up with a system of non linear differential equations driven by a time-dependent stimulus.Section 3 is divided in four parts. The first part analyses mathematically the potential anticipation effects ina general setting, before considering the role of ACells and lateral inhibition on anticipation (section 3.2) andgap junctions (section 3.3). Both sections contain general mathematical results, as well as numerical simula-tions for one dimensional motion. The fourth part investigates examples of two dimensional motions. Thelast section is devoted to discussion and conclusion. In Appendix A, we have added the values of parametersused in simulations, and, in Appendix B the receptive fields mathematical form used in the paper, as wellas the numerical method to compute efficiently the response of oriented two dimensional receptive fields tospatio-temporal stimuli. Appendix C presents a model of random connectivity from Amacrine to Bipolarcells inspired from biological data [71]. Finally, Appendix (D) contains mathematical results which constitutethe skeleton of the work, but whose proof would be too long to integrate in the core of the paper. This workis based on Selma Souihel’s PhD thesis where more extensive results can be found [67]. In particular, there isan analysis of the conjugated effects of retinal and cortical anticipation, subject of a forthcoming paper, andbriefly discussed in the conclusion.In all the following simulations, we use the CImg Library, an open-source C++ tool kit for image process-ing, in order to load the stimuli and reconstruct the retina activity. The source code is available on demand.
In the retinal processing light photons coming from a visual scene are converted into voltage variationsby photoreceptors (cones and rods). The complex hierarchical and layered structure of the retina allows toconvert these variations into spike trains, produced by Ganglion Cells (GCells) and conveyed to the thalamusvia their axons. We considerably simplify this process here. Light response induces a voltage variations ofBipolar cells (BCells), laterally connected via Amacrine cells (ACells), and feeding GCells, as depicted in ig. 1. We describe this structure in details here. Note that neither BCells nor ACells are spiking. They actsynaptically on each other by graded variations of their potential.We assimilate the retina to a flat, two dimensional square of edge length L mm. Therefore, we do notintegrate the dimensional structure of the retina in the model, merely for mathematical convenience. Spatialcoordinates are noted x, y (see Fig. 2 for the whole structure).In the model, each cell population tiles the retina with a regular square lattice. The density of cellsis therefore uniform for convenience but the extension to non uniform density can be afforded. For thepopulation p we note δ p the lattice spacing in mm, and N p the total number of cells. Without loss of gen-erality we assume that L , the retina’s edge size, is a multiple of δ p . We note L p = Lδ p , the number ofcells p per row or column so that N p = L p . Each cell in the population p has thus Cartesian coordinates ( x, y ) = ( i x δ p , i y δ p ) , ( i x , i y ) ∈ { , . . . , L p } . To avoid multiples indices, we associate to each pair ( i x , i y ) aunique index i = i x + ( i y − L p . The cell of population p , located at coordinates ( i x δ p , i y δ p ) is then denotedby p i . We note d (cid:2) p i , p (cid:48) j (cid:3) the Euclidean distance between p i and p (cid:48) j .We use the notation V p i for the membrane potential of cell p i . Cells are coupled. The synaptic weight fromcell p j to cell q i reads W p j q i . Thus, the pre-synaptic neuron is expressed in the upper index; the post-synaptic,in the lower index. Dynamics of cells is voltage-based. This is because our model is constructed from Chen etal model [17] itself derived from Berry et al [8] where a voltage-based description is used. Implicitly, voltageis measured with respect to the rest state of the cell ( V p i = 0 when the cell receives no input). The model consists first of a set of N B BCells, regularly spaced by a distance δ B , with spatial coordinates x i , y i , i = 1 . . . N . Their voltage, a function of the stimulus, is computed as follows. The projection of the visual scene on the retina ("stimulus") is a function S ( x, y, t ) where t is the timecoordinate. As we don’t consider color sensitivity here S characterizes a black and white scene, with a controlon the level of contrast ∈ [0 , . A Receptive Field (RF) is a region of the visual field (the physical space) inwhich stimulation alters the voltage of a cell. Thus, BCell i has a spatio-temporal receptive field K B i , featuringthe biophysical processes occurring at the level of the Outer Plexiform Layer (OPL), that is photo-receptors(rod-cones) response modulated by Horizontal Cells (HCells). As a consequence, in our model, the voltageof BCell i is stimulus-driven by the term: V i drive ( t ) = (cid:104) K B i x,y,t ∗ S (cid:105) ( t ) = (cid:90) + ∞ x = −∞ (cid:90) + ∞ y = −∞ (cid:90) ts = −∞ K ( x − x i , y − y i , t − s ) S ( x, y, s ) dx dy ds, (1)where x,y,t ∗ means space-time convolution. We consider only one family of BCells so that the kernel K is thesame for all BCells. What changes is the center of the RF, located at x i , y i , which also corresponds to thecoordinates of the BCell i . We consider in the paper separable kernel K ( x, y, t ) = K S ( x, y ) K T ( t ) where K S is the spatial part and K T the temporal part. The detailed form of K is given in Appendix B.We have : dV i drive dt = (cid:20) K B i x,y,t ∗ d S dt (cid:21) ( t ) , (2)resulting from the condition K B i ( x, y,
0) = 0 (see Appendix B). Note that the exponential decay of thespatial and temporal part at infinity ensures the existence of the space-time integral. The spatial integral (cid:82) R K S ( x, y ) S ( x, y, u ) dx dy is numerically computed using error function in the case of circular RF, and acomputer vision method from Geusenroek et al. [33] in the case of anisotropic RF, allowing to integrate gen-eralized Gaussians with an efficient computational time. This method is described in the Appendix, sectionB. Example of a retina grid tiling and indexing.
The green and blue ellipses denote re-spectively the positive center and the negative surround of the Bcell receptive field K S . The centerof RF coincides with the position of the cell (blue and green arrows). The red ellipse denotes theganglion cell pooling over bipolar cells (eq. (17)). For explanations purposes, we will often use the approximation of V i drive by a Gaussian pulse, withwidth σ , propagating at constant speed v along the direction (cid:126)e x : V i drive ( t ) = A √ π σ e −
12 ( x − vt )2 σ ≡ V √ π e −
12 ( x − vt )2 σ , (3)where x = k δ B is the horizontal coordinate of BCell i and where σ is in mm , A is in mV.mm (and isproportional to stimulus contrast), V is in mV . In our model, the BCell voltage is the sum of the external drive (1) received by the BCell and of a post-synaptic potential P B i induced by connected ACells: V B i ( t ) = V i drive ( t ) + P B i ( t ) . (4)The form of P B i is given by eq. (11) in the section 2.3.1. P B i ( t ) = 0 when no ACells are considered.BCells have voltage threshold [8]: N B ( V B i ) = (cid:26) , if V B i ≤ θ B ; V B i − θ B , else . (5)Values of parameters are given in appendix A.BCells have gain control, a desensitization when activated by a steady illumination [84]. This desensiti-zation is mediated by a rise in intracellular calcium Ca , at the origin of a feedback inhibition preventing hus prolonged signalling of the ON BCell [66, 17]. Following Chen et al., we introduce the dimensionlessactivity variable A B i obeying the differential equation: dA B i dt = − A B i τ a + h B N ( V B i ( t )) . (6)Assuming an initial condition A B i ( t ) = 0 at initial time t the solution is: A B i ( t ) = h B (cid:90) tt e − t − sτa N ( V B i ( s )) ds. (7)The bipolar output to ACells and GCells is then characterized by a non linear response to its voltagevariation, given by : R B i ( V B i , A B i ) = N B ( V B i ) G B ( A B i ) , (8)where : G B ( A B i ) = (cid:40) , if A B i ≤ A Bi , else . (9)Note that R B i has the physical dimension of a voltage, whereas, from eq. (9), the activity A B i is dimension-less. As a consequence, the parameter h B in eq. (6) must be expressed in ms − mV − . The form (9) and its -th power are based on experimental fits made by Chen et al. Its form is shown in Fig. 3.In the course of the paper we will use the following piecewise linear approximation also represented inFig. 3: G B ( A ) = , if A ∈ ] − ∞ , ∪ [ , + ∞ [ , Silent region ;1 , if A ∈ [0 , ] , Maximal gain ; − A + 2 , if A ∈ [ , ] , Fast decay . (10)Thanks to this approximation we roughly distinguish regions for the gain function G B ( A ) . This shape isuseful to understand the mechanism of anticipation (section 3.1). There is a wide variety of ACells (about 30-40 different types for humans) [57]. Some specific types arewell studied such as Starburst Amacrine Cells, which are involved in direction sensitivity [29, 74, 28], as wellas contrast impression and suppression of GCells response [51], or AII, a central element of the vertebraterod-cone pathway [47].Here, we don’t want to consider specific types of ACells with a detailed biophysical description. Instead,we want to point out the potential role they can play in motion anticipation, thanks to the inhibitory lateralconnectivity they induce. We focus on a specific circuitry involved in differential motion: an object with adifferent motion from its background induces more salient activity. The mechanism, observed in mice andrabbit retinas [54, 35] is featured in Fig. 1, pathway III. When the left pathway receives a different illumi-nation from the right pathway (corresponding e.g. to a moving object), this asymmetry is amplified by theACells’ mutual inhibition, enhancing the response of the left pathway in a "push-pull" effect. We want topropose that such a mutual inhibition circuit, deployed in a lattice through the whole retina, can generate -under specific conditions mathematically analysed - a wave of activity propagation triggered by the movingobject.In the model, ACells tile the retina with a lattice spacing δ A . We index them with j = 1 . . . N A . Gain control (9) as a function of activity A . The function l ( A ) , in dashed line, is apiecewise linear approximation of G B ( A ) from which regions are roughly defined. In the region"Silent" the gain vanishes so the cell does not respond to stimuli; in the region "Max", the gain ismaximal so that cell behaviour does not show any difference with a not gain-controlled cell; theregion "Fast decay" is the one which contributes to anticipation by shifting the peak in the cell’sactivity (see section 3.1). The value A c = corresponds to the value of activity where gain control,in the piecewise linear approximation, becomes effective. We consider here a simple model of ACells. We assimilate them to passive cells (no active ionic channels)acting as a simple relay between BCells. This aspect is further discussed later in the paper. The ACell A j ,connected to the BCell B i , induces on the latter the post-synaptic potential : P A j B i ( t ) = W A j B i ( t ) (cid:90) t −∞ γ B ( t − s ) V A j ( s ) ds ; γ B ( t ) = e − tτB H ( t ) , where the Heaviside function H ensures causality. Thus, the post synaptic potential is the mere convolutionof the pre synaptic ACell voltage, with an exponential α -profile [25]. In addition, we assume the propagationto be instantaneous.Here, the synaptic weight W A j B i < mimics the inhibitory connection from ACell to BCell (glycine orGABA) with the convention that W A j B i = 0 if there is no connection from A j to B i .In general, several ACells input the BCell B i giving a total PSP: P B i ( t ) = N B (cid:88) j =1 W A j B i (cid:90) t −∞ γ B ( t − s ) V A j ( s ) ds. (11)Conversely, the BCell B i connected to A j induces, on this cell, a synaptic response characterized by apost-synaptic potential (PSP) P A j ( t ) . As ACells are passive elements their voltage V A j ( t ) is equal to this PSP.We have thus: V A j ( t ) = N A (cid:88) i =1 W B i A j (cid:90) t −∞ γ A ( t − s ) R B i ( s ) ds, (12) ith γ A ( t ) = e − tτA H ( t ) . Here, W B i A j > corresponding to the excitatory effect of BCells on ACells, througha glutamate release. Note that the voltage of the BCell is rectified and gain-controlled. The coupled dynamics of Bipolar and Amacrine cells can be described by a dynamical system that wederive now.
Bipolar voltage.
By differentiating (11), (4), and introducing: F B i ( t ) = (cid:20) K B i x,y,t ∗ (cid:18) S τ B + d S dt (cid:19) (cid:21) ( t ) = V i drive τ B + dV i drive dt , (13)we end up with the following equation for the bipolar voltage: dV B i dt = − τ B V B i + N A (cid:88) j =1 W A j B i V A j + F B i ( t ) , (14)where we have used (2). This is a differential equation driven by the time dependent term F B i containing thestimulus and its time derivative.To illustrate the role of F B i , let us consider an object moving with a speed (cid:126)v depending on time, thus witha non zero acceleration (cid:126)γ = d(cid:126)vdt . This stimulus has the form S ( t ) = g (cid:16) (cid:126)X − (cid:126)v ( t ) t (cid:17) , with (cid:126)X = (cid:18) xy (cid:19) , so that d S dt = − (cid:126) ∇ g (cid:16) (cid:126)X − (cid:126)v ( t ) t (cid:17) . ( (cid:126)v + (cid:126)γt ) where (cid:126) ∇ denotes the gradient. Therefore, thanks to the eq. (14), BCells are sensitive to changes in directions, thereby justifying a study of dimensional stimuli (Section 3.4). Note thatthis property is inherited from the simple, differential structure of the dynamics, the term dV idrive dt resultingfrom the differentiation of V B i . This term does not appear in the classical formulation (1) of the bipolarresponse, without amacrine connectivity. It appears here because synaptic response involves an implicit timederivative via the convolution (12). Coupled dynamics.
Likewise, differentiating (12) gives: dV A j dt = − τ A V A j + N B (cid:88) i =1 W B i A j R B i . (15)Eq. (6) (activity), (14) and (15) define a set of N B + N A differential equations, ruling the behaviour ofcoupled BCells and ACells, under the drive of the stimulus, appearing in the term F B i ( t ) . We summarize thedifferential system here: dV Bi dt = − τ B V B i + (cid:80) N A j =1 W A j B i V A j + F B i ( t ) , dV Aj dt = − τ A V A j + (cid:80) N B i =1 W B i A j R B i , dA Bi dt = − A Bi τ a + h B N ( V B i ) . (16)We have used the classical dynamical systems convention where time appears explicitly only in the drivingterm F B i ( t ) to emphasize that (16) is non-autonomous. Note that BCells act on ACells via a rectified volt-age (gain control and piecewise linear rectification), in agreement with fig. 1, pathway III. We analyse thisdynamics in section 3.2.1. .3.3 Connectivity graph The way ACells connect to BCells, and reciprocally, have a deep impact on the dynamics (16). In thispaper, we want to point out the role of relative excitation (from BCells to ACells) and inhibition (from ACellsto BCells) as well as the role of the network topology. For mathematical convenience - dealing with squarematrices - we assume from now on that there are as many BCell as ACells and we set N ≡ N A = N B . Atthe core of our mathematical studies is a matrix, L , defined in section 3.2.1, whose spectrum conditions theevolution of the BCells-ACells network under the influence of a stimulus. It is interesting and relevant torelate the spectrum of L to the spectrum of the connectivity matrices ACells to BCells and BCells to ACells.There is not such general relation for arbitrary matrices of connectivity. A simple case holds when the twoconnectivity matrices commute. Here, we choose an even simpler situation, based on the fact that we com-pare the role of the direct feed-forward pathway on anticipation in the presence of ACells lateral connectivity.We feature the direct pathway by assuming that a BCell connects only one ACell with a weight w + uniformfor all BCell, so that W BA = w + I N,N , w + > , where I N,N is the N -dimensional identity matrix. In contrast,we assume that ACells connect to BCells with a connectivity matrix W , not necessarily symmetric, with auniform weight − w − , w − > , so that W AB = − w − W .We consider then two types of network topology for W :1. Nearest neighbours.
An ACell connects its d nearest BCell neighbours where d = 1 , is the latticedimension.2. Random ACell connectivity.
This model is inspired from the paper [71] on the shape and arrangementof starburst ACells in the rabbit retina. Each cell (ACell and BCell) has a random number of branches(dendritic tree), each of which has a random length and a random angle with respect to the horizontalaxis. The length of branches L follow an exponential distribution with spatial scale ξ . The number ofbranches n is also a random variable, Gaussian with mean ¯ n and variance σ n . The angle distributionis taken to be isotropic in the plane, i.e. uniform on [0 , π [ . When a branch of an ACell A intersects abranch of a BCell B there is a chemical synapse from A to B. The probability that two branches intersectfollows a nearly exponential probability distribution that can be analytically computed (see Appendix,section C). There are many different types of GCells in the retina, with different physiologies and functions [3, 63].In the present computational study we focus on specific subtypes associated to the pathways I-II (Fast OFFcells with gain control), III (Differential Motion Sensitive cells), IV (Direction selective cells), in Fig. 1. Allthese have common features: BCells pooling and gain control.
In the retina, GCells of the same type cover the surface of the retina, forming a mosaic. The degree ofoverlap between GCells indicates the extent to which their dendritic arbours are entangled in one another.This overlap remains however very limited between cells of the same type [61]. We note k the index of theGCells, k = 1 . . . N G and δ G the spacing between two consecutive GCells lying on the grid (Fig. 2).In the model, GCell k pools over the output of BCells in its neighbourhood [17]. Its voltage reads: V ( P ) G k = (cid:88) i W B i G k R B i , (17)where the superscript "P" stands for "Pool". We use this notation to differentiate this voltage from the totalGCell voltage, V G k , when they are different. This happens in the case when GCells are directly coupled by ap junctions (sections 2.4.4, 3.3). When there is no ambiguity we will drop the superscript "P". In eq. (17),the weights W B i G k are Gaussian: W B i G k = a p e − d [ Bi, Gk ] σ p . (18)where σ p has the dimension of a distance and a p is dimensionless. The voltage V G k is processed through a gain control loop similar to the BCell layer [17]. As GCells arespiking cells, a non-linearity is fixed so as to impose an upper limit over the firing rate. Here, it is modelledby a sigmoid function, e.g. : N G ( V ) = , if V ≤ θ G ; α G ( V − θ G ) , if θ G ≤ V ≤ N maxG /α G + θ G ; N maxG , else . (19)This function corresponds to a probability of firing in a time interval. Thus, it is expressed in Hz . Conse-quently, α G is expressed in Hz mV − and N maxG in Hz . Parameters values can be found in the appendixA. Gain control is implemented with an activation function A G k , solving the following differential equation: dA G k dt = − A G k τ G + h G N G ( V G k ) , (20)and a gain function : G G ( A ) = (cid:26) , if A < A , else . (21)Note that the origin of this gain control is different from the BCell gain control (9). Indeed, Chen et al.hypothesize that the biophysical mechanisms that could lie behind ganglion gain control are spike-dependentinactivation of Na + and K + channels, while the study by Jacoby et al. [38] hypothesize that GCells gaincontrol is mediated by feed-forward inhibition that they receive from ACells. The specific forms of the non-linearity and the gain control function used in this paper match however the first hypothesis, namely thesuppression of the Na + current [17].Finally, the response function of this GCell type is: R G ( V G k , A G k ) = N G ( V G k ) G G ( A G k ) . (22)In contrast to BCell response R B , (8), which is a voltage, here R G is a firing rate.Gain control has been reported for OFF GCells only [8] [17]. Therefore, we restrict our study to OFF cells,i.e with a negative center of the spatial RF kernel. However, on mathematical grounds, it is easier to carry ourexplanation when the RF center is positive. Thus, for convenience, we have adopted a change in conventionin terms of contrast measurement. We take the reference value 0 of the stimulus to be white rather than black,black corresponding then to 1. The spatial RF kernel is also inverted, with a positive center and a negativesurround. The problem is therefore mathematically equivalent to an ON Cell submitted to positive stimulus. We consider here a class of GCells, connected to ACells according to pathways III in fig. 1, acting as dif-ferential motion detectors. They are able to respond saliently to an object moving over a stationary surround,while being strongly inhibited by global motion. Here, stationary is meant in a general, probabilistic sense.This can be a uniform background, or a noisy background where the probability distribution of the noiseis time-translation invariant. These cells are hence able to filter head and eye movements. Baccus et al. [1] mphasized a pathway accountable for this type of response, involving polyaxonal ACells which selectivelysuppress GCells response to global motion and enhance their response to differential motion as shown in Fig.1, pathway III. The GCell receives an excitatory input from the BCells lying in its receptive field which re-spond to the central object motion, and an indirect inhibitory input from ACells that are connected to BCellswhich respond to the background motion. When motion is global, the excitatory signal is equivalent to the in-hibitory one, resulting in an overall suppression. However, when the object in the center moves distinctivelyfrom the surrounding background, the cell in the center responds strongly.There are here two concomitant effects. When a moving object (say, from left to right) enters the BCellpool connected to a central GCell k D , the BCells in the periphery of the pool respond first, with no significantchange on the GCell response, because of the Gaussian shape (18) of the pooling: weights are small in theperiphery. Those BCells excite however the ACells they are connected to, with the effect of inhibiting theBCells of neighbouring GCells pools. This has the effect of decreasing the voltage of these BCells, which inturn excite less ACells which, in turn, inhibit less the BCells of the pool k D . Thus, the response of the GCell k D is enhanced, while the cells on the background are inhibited. We call this effect "push-pull" effect. Notethat propagation delays ought to play an important role here, although we are not going to consider them inthis paper. These cells correspond to the pathway IV in Fig. 1. They are only coupled via electric synapses (gapjunctions). In several animals, like the mouse, this enables the corresponding GCells to be direction sensitive.Note that other mechanisms, involving lateral inhibition via Starburst Amacrine Cells have also been widelyreported [29, 74, 28, 81, 78, 65, 64]. Here we focus on gap junctions direction sensitive cells (DSGCs). Thereexist four major types of these DSGCs, each responding to edges moving in one of the four cardinal direc-tions. Trenhlom et al. [73] have emphasized the role of these cells coupling in lag normalization: uncoupledcells begin responding when a bar enters their receptive field, i.e, their dendritic field extension, whereascoupled cells start responding before the bar reaches their dendritic field. This anticipated response is due tothe effective propagation of activity from neighbouring cells through gap junctions, and is particularly inter-esting when comparing the responses for different velocities of the bar. Trenhlom et al. have shown that theuncoupled DSGCs detect the bar at a position which is further shifted as the velocity grows, while coupledcells respond at an almost constant position, regardless of the velocity. In our work, we analyse this effectin terms of a propagating wave driven by the stimulus and show that, temporally, this spatial lag normal-ization induces a motion extrapolation that confers to the retina more than just the ability to compensate forprocessing delays, but to anticipate motion.Classical, symmetric bidirectional gap junctions coupling between neighbouring cells would involve acurrent of the form − g ( V G k − V G k − ) − g ( V G k − V G k +1 ) where g is the gap junction conductance. In contrast,here, the current takes the form − g ( V G k − V G k − ) . This is due to the specific asymmetric structure of thedirection selective GCell dendritic tree [73]. The experimental results of these authors suggest that the effectof the possible gap junction input from downstream cells, in the direction of motion, can be neglected due tooffset inhibition and gain control suppression. This, along with the asymmetry of the dendritic arbour, justifythe approximation whereby the cell k+1 doesn’t influence the current in the cell k. This induces a strongdifference in the propagation of a perturbation. Indeed, consider the case V G k − V G k − = V G k − V G k +1 = δ .In the symmetric form the total current vanishes whereas in the asymmetric form the current is − gδ . Still,the current can have both directions depending on the sign of δ . This has a strong consequence on the wayGCells connected by gap junctions respond to a propagating stimulus, as shown in section 3.3.The total GCell voltage is the sum of the pooled BCell voltage V ( P ) G k and of the effect of neighbours GCellsconnected to k by gap junctions: V G k ( t ) = V ( P ) G k − gC (cid:90) t −∞ ( V G k ( s ) − V G k − ( s )) ds here C is the membrane capacitance. Deriving the previous equation with respect to time, we obtain thefollowing differential equation governing the GCell voltage: dV G k dt = dV ( P ) G k dt − w gap (cid:2) V G k ( t ) − V G k − ( t ) (cid:3) , (23)where: w gap = gC . (24)Gain control is then applied on V G k as in (22). An alternative is to consider that gain control occurs beforegap junctions effect. We investigated this effect as well (not shown, see [67]). Mainly, the anticipatory effect isenhanced when the gain control is applied after the gap junction coupling, thus, from now, we focus on theformulation (23) in the paper.Note that our voltage-based model of gap junctions takes a different from as Trenholm et. al (expressedin terms of currents), because we had to adapt it so as to deal with the pooling voltage form (17). Still, ourmodel reproduces the lag normalization as in the original model as we checked (not shown, see [67]). The (smooth) trajectory of a moving object can be extrapolated from its past position and velocity toobtain an estimate of its current location [49, 4, 50]. When human subjects are shown a moving bar travellingat constant velocity, while a second bar is briefly flashed in alignment with the moving bar, the subjects reportseeing the flashed bar trailing behind the moving bar. This led Berry et al [8] to investigate the potential roleof the retina in anticipation mechanisms. Under constraints on the bar’ speed and contrast they were able toexhibit a positive anticipation time, defined as the time lag between the peak in the retinal GCell response toa flashed bar and the corresponding peak when the stimulus is a moving bar.In this paper we adopt a slightly different definition although inspired from it. Indeed, the goal of thismodelling paper is to dissect the various potential stages of retinal anticipation as developed in the nextsubsections.Several layers and mechanisms are involved in the model, each one defining a response time and poten-tially contributing to anticipation, under conditions that we now analyse.
We consider first a single BCell, without lateral connectivity so that V B i = V i drive . The very mechanismof anticipation at this stage is illustrated in Fig. 4. The peak response time of the convolution of the stimuluswith the RF of one BCell occurs at a time t B (dashed line in Fig. 4 a). The increase in V i drive leads to anincrease in activity (Fig. 4, c) and an increase of R B (Fig. 4, e). When activity becomes large enough, gaincontrol switches on (Fig. 4 d) leading to a sharp decrease of the response R B (Fig. 4 e) and a peak in R B occurring at time t B A (dashed line in Fig. 4 e) before t B . The bipolar anticipation time , ∆ B = t B − t B A , istherefore positive.Mathematically, ∆ B > results from the intermediate value theorem using that dV idrive dt ≥ on [ 0 , t B ] and that t B A is defined by: dV i drive dt (cid:12)(cid:12)(cid:12)(cid:12) t = t BA = − V i drive ( t B A ) G (cid:48) B ( A B i ) G B ( A B i ) dA B i dt (cid:12)(cid:12)(cid:12)(cid:12) t = t BA , The mechanism of motion anticipation and the role of gain control.
The figure illus-trates the bipolar anticipation time ∆ B , without lateral connectivity. We see the response of OFFBCells with gain control to a dark moving bar. The curves correspond to three cells spaced by µm . The first line (a) shows the linear filtering of the stimulus, corresponding to V drive ( t ) (eq.(1)). The line (b) corresponds to the threshold non-linearity N B applied to the linear response; (c)represents the adaptation variable (16), and (d) shows the gain control time curse. Finally, the lastline (e) corresponds to the response R B i of the BCell. The two dashed lines correspond respectivelyto t B and t B A , the peak in the response of the (purple) Bcell without pooling.14 here the right hand side is positive provided that the parameters h B , τ a are tuned such that dA Bi dt ≥ on [0 , t B ] . An important consequence is that the amplitude of the response at the peak is smaller in the presenceof gain control (compare the amplitude of the voltage in Fig. 4, a to 4, e).The anticipation time at the BCells level depends on parameters such as h B , τ a . It depends as well oncharacteristics of the stimulus such as contrast, size and speed. An easy way to figure this out is to considerthat the peak in BCell response (Fig. 4 d, e) arises when the gain control function G B ( A B i ) starts to drop off(Fig. 4 e), which, from the piecewise linear approximation (10) of BCell arises when A = . When V i drive has the form (3) this gives, using N ( V i drive ) = V i drive , (7), and letting the initial time t → −∞ (whichcorresponds to assuming that the initial state was taken in a distant past, quite longer than the time scales inthe model): A B i ( t B A ) = A h B v e σ τ a v e τav ( x − v t BA ) (cid:20) − Π (cid:18) x − v t B A σ + στ a v (cid:19) (cid:21) = 23 , (25)where Π( x ) is the cumulative distribution function of the standard Gaussian probability (see definition, eq.(60) in the appendix). This establishes an explicit equation for the time t B A as a function of contrast ( A ), size( σ ), and speed ( v ) as well as the parameters h B and τ a . We do not show the corresponding curves here (see[67] for a detailed study) preferring to illustrate the global anticipation at the level of GCells, illustrated inFig. 5 below. The main effects we want to illustrate in the paper (impact of lateral connectivity on GCells anticipation)are evidenced by the shift of the peak in activity of the BCells pooled voltage, occurring at time t G . Wefocus on this time here, postponing to section 3.1.3 the subsequent effect of GCells gain control. We assumetherefore here that h G = 0 so that A G k = 0 and G G ( A G k ) = 1 in (19). Thus, the firing rate of Gcell k is N G ( V G k ) . For mathematical simplicity we will consider that the firing rate function (5) of G is a smooth,monotonously increasing sigmoid function so that N (cid:48) G ( V G k ) > . We define t G as the time when V G k ismaximum, after the stimulus is switched on. This corresponds to dV Gk dt = 0 and d V Gk dt < . Equivalently,from equations (17), (23): (cid:88) i W B i G k dR B i dt = (cid:88) i W B i G k (cid:20) G B ( A B i ) N (cid:48) B ( V B i ) dV B i dt + N B ( V B i ) G (cid:48) B ( A B i ) dA B i dt (cid:21) = w gap (cid:2) V G k − V G k − (cid:3) , (26)where this equation holds at time t = t G (we have not written explicitly t G to alleviate notation). This is themost general equation for the anticipation time at the level of BCells pooling.In the sum (cid:80) i , there are two types of BCells. The inactive ones where V B i ≤ Θ B , N B ( V B i ) = 0 and dR Bi dt = 0 so they do not contribute to the activity. The active BCells, V B i > Θ B , obey N B ( V B i ) = V B i .For the moment we assume that, at time t G , there is no Bcell switching from one state (active/inactive) to the other ,postponing this case to the end of the section. Then, eq. (26) reduces to: (cid:80) i W B i G k (cid:124) (cid:123)(cid:122) (cid:125) ( V ) G B ( A B i ) (cid:124) (cid:123)(cid:122) (cid:125) ( II ) − τ B V B i + N A (cid:88) j =1 W A j B i V A j (cid:124) (cid:123)(cid:122) (cid:125) ( III ) + F B i ( t ) (cid:124) (cid:123)(cid:122) (cid:125) ( I ) = − (cid:80) i W B i G k (cid:124) (cid:123)(cid:122) (cid:125) ( V ) G (cid:48) B ( A B i ) (cid:124) (cid:123)(cid:122) (cid:125) ( II ) V B i ( t ) dA Bi dt + w gap (cid:2) V G k − V G k − (cid:3)(cid:124) (cid:123)(cid:122) (cid:125) ( IV ) . (27) From (6) dA Bi dt > if A i ( t ) < h B τ a V i drive ( t ) . This essentially requires τ a to be slow enough. his general equation emphasizes the respective role of (I), stimulus (term F B i ( t ) ); (II), gain control (terms G B ( A B i ) , G (cid:48) B ( A B i ) ); (III), ACell lateral connectivity (term W A j B i ); (IV), gap junctions (term w gap (cid:2) V G k ( t (cid:48) G A ) − V G k − ( t (cid:48) G A ) (cid:3) );(V), pooling (terms W B i G k ). Note that we could as well consider a symmetric gap junctions connectivity wherewe would have a term w gap (cid:2) − V G k +1 + 2 V G k − V G k − (cid:3) in IV. The equation terms has been arranged thisway for reasons that become clear in the next lines. It is not possible to solve this equation in full generalitybut it can be used to understand the respective role of each component.In the absence of gain control and lateral connectivity ( W A j B i = 0 , w gap = 0 ) the peak in GCell G k voltage,at time t (cid:48) G is given by: (cid:88) i W B i G k dV i drive dt = 0 , (28)This generalizes the definition of t B , time of peak of a single BCell, to a set of pooled BCells and we willproceed along the same lines as section 3.1.1. We fix as reference time the time when the pooled voltagebecomes positive. It increases then until the time t (cid:48) G when (cid:80) i W B i G k dV idrive dt = 0 . Thus, (cid:80) i W B i G k dV idrive dt ispositive on [0 , t (cid:48) G [ and vanishes at t (cid:48) G .We now show that, in the presence of gain control, the peak occurs at time t G < t (cid:48) G leading to anticipationinduced by gain control and generalizing the effect observed for one Bcell in section 3.1.1. Indeed, equation(27) reads now: (cid:88) i W B i G k G B ( A B i ) dV i drive dt = − (cid:88) i W B i G k G (cid:48) B ( A B i ) V i drive ( t ) dA B i dt . (29)Because ≤ G B ( A B i ) ≤ , (cid:80) i W B i G k G B ( A B i ) dV idrive dt ≤ (cid:80) i W B i G k dV idrive dt so that the left hand side in(29) reaches at a time t G ≤ t (cid:48) G . The right hand side is positive for the same reasons as in section 3.1.1. Thesame mathematical argument holds as well, using the intermediate value theorem, to show that t G < t (cid:48) G .We now investigate eq. (27) with the two terms of lateral connectivity: (III), ACells and, (IV) gap junctions.The effect of gap junctions is straightforward. A positive term w gap (cid:2) V G k − V G k − (cid:3) increases the right handside of eq. (27). As developed in section 3.3 this arises when the stimulus propagates in the preferred directionof the cell inducing a wave of activity propagating ahead of the stimulus. In view of the qualitative argumentdeveloped above using the intermediate value theorem, this can enhance the anticipation time. This deserveshowever a deeper study developed in section 3.3.The effect of ACells cells is less evident, as the term (cid:16) − τ B V B i + (cid:80) N A j =1 W A j B i V A j + F B i ( t ) (cid:17) can have anysign, so that network effect can either anticipate or delay the ganglion response, as illustrated in several ex-amples in the next section. As we show, this term is in general related to a wave of activity, enhancing orweakening the anticipation effect as shown in section 3.2.Let us finally discuss what happens when some BCell switches from one state (active/inactive) to theother (i.e. V B i = Θ B ). In this case, taking into account the definition (5), the derivative N (cid:48) B ( V B i ) = . Thus,when a BCell reaches the lower threshold, there is a big variation in N (cid:48) B ( V B i ) thereby leading to a positivecontribution in (26) and an additional term (cid:80) i W B i G k G B ( A B i ) (cid:16) − τ B V B i + (cid:80) N A j =1 W A j B i V A j + F B i ( t ) (cid:17) in the left hand side of (27), where the sumholds on switching state cells. As we see in section (3.2) this can have an important impact on the anticipationtime. .1.3 Anticipation time at the GCells level We now show that the firing rate of the GCell k , given by (22), reaches its maximum at a time t G A < t G .From (22), at time t G A : dV G k dt = V G k A G k dA G k dt . (30) V G k starts from and increases on the time interval [ 0 , t G ] thus dV Gk dt is positive on [ 0 , t G ] and vanishes at t G . Thus, there is a time t d < t G such that dV Gk dt increases on [0 , t d ] and decreases on [ t d , t G ] . The right handside of (30) starts from at t = 0 and stays strictly positive until, either V G k vanishes which occurs for t > t G ,or until dA Gk dt vanishes. We choose the characteristic time τ G and the intensity h G in (20) so that dA Gk dt > on [ 0 , t G ] . Thus, V Gk A Gk dA Gk dt > on [ 0 , t G ] . Therefore, in the time interval [ t d , t G ] , dV Gk dt decreases to while V Gk A Gk dA Gk dt increases from . From the intermediate value theorem these two curves have to intersect at atime t G A < t G .We finally define the total anticipation time of a Gcell as: ∆ = t B c − t G A , (31)where t B c is the peak of the BCell at the center of the BCells pooling to that GCell. In general, ∆ depends on gain control, lateral connectivity, as well as characteristics of the stimulus suchas speed and contrast. This has been shown mathematically in eq. (25) for a single BCell. Here, we investigatenumerically the dependence of the total anticipation time of a Gcell when the stimulus is a bar of infiniteheight, width σ mm, travelling in one dimension at speed v mm/s with contrast C ∈ [ 0 , . Results areshown in fig. 5. This figure is a calibration later used to compare to the effects induced by lateral connectivity.We first observe that anticipation increases with contrast, as it has experimentally been observed [8].Indeed, increasing the contrast increases V i drive ( t ) thereby accelerating the growth of A i so that gain controltakes place earlier (Fig 5 a). We also notice that anticipation increases with the width of the object until amaximum (Fig 5 b). Finally, the model shows a decrease in anticipation as a function of velocity, as it wasevidenced experimentally [8, 40] (Fig 4 c). Indeed, when the velocity increases, V drive varies faster than thecharacteristic activation time τ a and the adaptation peak value is lower. Consequently, gain control has aweaker effect and the peak activity is less shifted than when the bar is slow.A large part of these effects can be understood from eq. (25). Note however here that simulation of Fig.5 takes into account the convolution of a moving bar with the receptive field, the pooling effect, and gaincontrol at the stage of GCells.In Fig. 5 we also show the evolution of GCells maximum firing rate as a function of the moving barvelocity, contrast and size. We observe that it increases with these parameters, an expected result. In this section we study the potential effect of ACells (pathway III of Fig. 1) on motion anticipation. Werestrict to the case where there are as many BCells as ACells ( N B = N A ≡ N ) so that the matrices W AB and W BA are square matrices. We first derive general mathematical results (for the full derivation, see appendixsection D) before considering the two types of connectivity described in section 2.3.3. Maximum firing rate and anticipation time variability with stimulus parameters in thegain control layer of the model.
Left: contrast (with v = 1 mm/s et size = µ m ); Middle: size(with v = 2 mm/s et contrast = ); Right: speed (with contrast and size = µ m ) . We study mathematically the dynamical system (16) that we write in a more conve-nient form. We use Greek indices α, β, γ = 1 . . . N and define the state vector X as (cid:126) X α = V B i , α = i, i = 1 . . . N ; V A i , α = N + i, i = 1 . . . N ; A i , α = 2 N + i, i = 1 . . . N. Likewise, we define the stimulus vector (cid:126) F α = F B i , if α = 1 . . . N and (cid:126) F α = 0 otherwise. Then, the dynamicalsystem (16) has the general form: d (cid:126) X dt = H ( (cid:126) X ) + (cid:126) F ( t ) . (32)where H ( (cid:126) X ) is a non linear function, via the function R B i ( V B i , A B i ) of eq. (8), featuring gain control andlow voltage threshold. The non linear problem can be simplified using the piecewise linear approximation(10). Indeed, there is a domain of R N : Ω = (cid:26) V B i ≥ θ B , A B i ∈ (cid:20) , (cid:21) , i = 1 . . . N (cid:27) , (33)where R B i ( V B i , A B i ) = V B i so that (16) is linear and can be written in the form: d (cid:126) X dt = L . (cid:126) X + (cid:126) F ( t ) . (34)with: L = − I N,N τ B W AB N,N W BA − I N,N τ A N,N h B I N,N N,N − I N,N τ a (35)where I N,N is the N × N identity matrix and N,N is the N × N zero matrix. This corresponds to intermediateactivity, where neither BCells gain control (9) nor low threshold (5) are active. We first study this case anddescribe then what happens when trajectories of (32) get out of this domain, activating low voltage thresholdor gain control.The idea of using such a phase space decomposition with piecewise linear approximations has been used,in a different context by S. Coombes et al [20] and in [11, 15, 12]. e consider the evolution of the state vector (cid:126) X ( t ) from an initial time t . Typically, t is a reference timewhere the network is at rest, before the stimulus is applied. So, the initial condition (cid:126) X ( t ) will be set to without loss of generality. Linear analysis.
The general solution of (34) is: (cid:126) X ( t ) = (cid:90) tt e L ( t − s ) . (cid:126) F ( s ) ds, (36)The behaviour of the solution (36) depends on the eigenvalues λ β , β = 1 . . . N of L and its eigenvectors, (cid:126) P β , with entries P αβ . The matrix P transforms L in Jordan form ( L is not diagonalizable when h B (cid:54) = 0 , seeappendix D.1). Whatever the form of the connectivity matrices W BA , W AB the N last eigenvalues are always λ β = − τ a , β = 2 N + 1 . . . N .In appendix D.1 we show the following general result (not depending on the specific form of W BA , W AB ,they just need to be square matrices and to be diagonalizable): X α ( t ) = V α drive ( t ) + E BB,α ( t ) + E BA,α ( t ) + E Ba,α ( t ) , α = 1 . . . N, (37)where the drive term (1) is extended here to N -dimensions with V α drive ( t ) = 0 if α > N . The other termshave the following definition and meaning: E BB,α ( t ) = N (cid:88) β =1 (cid:18) τ B + λ β (cid:19) N (cid:88) γ =1 P αβ P − βγ (cid:90) tt e λ β ( t − s ) V γ drive ( s ) ds, α = 1 . . . N, (38)corresponds to the indirect effect, via the ACells connectivity, of the BCells drive on BCells voltages (i.e. thedrive excites BCell i which acts on BCell j via the ACells network); E BA,α ( t ) = N (cid:88) β = N +1 (cid:18) τ B + λ β (cid:19) N (cid:88) γ =1 P αβ P − βγ (cid:90) tt e λ β ( t − s ) V γ drive ( s ) ds, α = N + 1 . . . N, (39)corresponds to the effect of BCell drive on ACell voltages, and: E Ba,α ( t ) = h B N (cid:88) β =1 N (cid:88) γ =1 P α − Nβ P − βγ λ β + τ B λ β + τ a (cid:90) tt e λ β ( t − s ) V γ drive ( s ) ds + − τ B + τ a λ β + τ a A α − N ( t ) , α = 2 N +1 . . . N, (40)corresponds to the effect of the BCells drive on the dynamics of BCell activity variables. The first term of (40)corresponds to the action of BCells and ACells on the activity of BCells, via lateral connectivity. In the secondterm: A α − N ( t ) = (cid:90) tt e − t − sτa V α − N drive ( s ) ds (41)corresponds to the direct effect of the BCell voltage with index α − N on its activity (see eq. (7)).To sum up, equation (37) describes the direct effect of a time dependent stimulus (first term) and theindirect lateral network effects it induces. The term E Ba,α ( t ) is what activates the gain control. In the piecewiselinear approximation (10), the BCell i triggers its gain control when its activity: E Ba,α ( t ) > , α = 2 N + i. (42)This relation extends the computation, made in section 3.1.1 for isolated BCells, to the case of a BCell underthe influence of ACells. On this basis, let us now discuss how the network effect influences the activation of .00.10.20.30.40.50.60.70.80.9-2 -1.5 -1 -0.5 0 0.5 1 1.5 2 f r on t ( ,t ) λ β <00.00.51.01.52.02.53.0-2 -1.5 -1 -0.5 0 0.5 1 1.5 2 f r on t ( ,t ) t (ms) λ β >0 Figure 6:
Front (43) for different values of λ β (Purple) as a function of time, for the cell γ = 0 . Allfigures are drawn with v = 2 mm/s; σ = 0 . mm. Top. λ β = − . ms − ; Bottom. λ β = 0 . ms − . gain control and, thereby, anticipation.The structure of the terms (38), (39) (40) is interpreted as follows. The drive (index γ = 1 . . . N ) excitesthe eigenmodes β = 1 . . . N of L , with a weight proportional to P − βγ . The mode β , in turn excites thevariable α = 1 . . . N with a weight proportional to P αβ . The time dependence and the effect of the drive arecontrolled by the integral (cid:82) tt e λ β ( t − s ) V γ drive ( s ) ds . For example, when the stimulus has the Gaussian form(3) and cells are spaced with a distance δ so that cell γ is located at x = γ δ , we have, taking t → −∞ : (cid:90) t −∞ e λ β ( t − s ) V γ drive ( s ) ds = A v e σ λ βv e − λβv ( γ δ − v t ) Π (cid:20) λ β σv − σ ( γ δ − v t ) (cid:21) , (43)where Π( x ) is the cumulative distribution function of the standard Gaussian probability (see definition, eq.(60) in the appendix). This is actually the same computation as (25) with λ β = − τ a . Eq. (43) corresponds toa front, separating a region where Π [ . . . ] = 0 from a region where
Π [ . . . ] = 1 , propagating at speed v withan interface of width σ multiplied by an exponential factor e − λβv ( γ δ − v t ) . Here, the sign of the real part of λ β , λ β,r is important. If λ β,r < the front has the shape depicted in Fig. 6 top. It decays exponentially fast as t → + ∞ , with a time scale | λ β,r | . On the opposite, it increases exponentially fast, with a time scale λ β,r as t → + ∞ when λ β,r > , thereby enhancing the network effect and accelerating the activation of non lineareffect (low threshold or gain control) leading the trajectory out of Ω . Remark that the peak of the drive occursat γ δ − v t = 0 . The inflexion point of the function Π( x ) is at x = 0 . Thus, when λ β < the front is a bitbehind the drive, whereas it is a bit ahead when λ β > .Having unstable eigenvalues is not the only way to get out of Ω . Indeed, even if all eigenvalues are stablethe drive itself can lead some cells to get out of this set. When the trajectory of the dynamical system (34) getsout of Ω two cases are then possible:(i) Either a BCell i is such that V B i < θ B . In this case, R B i ( V B i , A B i ) = 0 . Then, in the matrix L there is aline of zeros replacing the line i in the matrix W BA , i.e. at the line i + N of L . This corresponds to a stableeigenvalue − τ A for L , controlling the exponential instability observed in Fig. 6 bottom. Thus, too lowBCells voltages trigger a re-stabilisation of the dynamical system. ii) There are BCells such that condition (42) holds, then gain control is activated and the system (32) be-comes non linear. Here we get out of the linear analysis and we have not been able to solve the problemmathematically. There is however a simple qualitative argument. If the cell i enters in the gain controlregion then the corresponding line i + N in the matrix W BA of L is replaced by W B i A G B ( A B i ) whichrapidly decays to (see e.g. Fig. 4 e). From the same argument as in (i) this generates a stable eigenvalue ∼ − τ A controlling as well the exponential instability.Eq. (37) features therefore the direct effect of the stimulus as well as the indirect effect, via the amacrinenetwork, corresponding to a weighted sum of propagating fronts, generated by the stimulus, and influenc-ing a given cell through the connectivity pathways. These fronts interfere, either constructively, inducing awave of activity enhancing the effect of the stimulus and, thereby, anticipation , or destructively somewhatlowering the stimulus effect. The fine tuning between "constructive" and "destructive" interferences dependson the connectivity matrix via the spectrum of L and its projection vectors (cid:126) P β . For example, complex eigen-values introduce time oscillations which are likely to generate destructive interferences, unless some specificresonances conditions exist between the imaginary parts of the eigenvalues λ β . Such resonances are knownto exist e.g. in neural network models exhibiting a Ruelle-Takens transition to chaos [58], and they are closelyrelated to the spectrum of the connectivity matrix [13]. Although we are not in this situation here, our linearanalysis clearly shows the influence of the spectrum of L , itself constrained by W , on the network responseto stimuli and anticipation. Spectrum of L . This argumentation invites us to consider different situations where one can figure outhow connectivity impacts the spectrum of L and thereby anticipation. We therefore provide some generalresults about the spectrum of L and potential linear instabilities before considering specific examples. Theseresults are proved in the appendix D.2. As stated in section 2.3.3, to go further in the analysis we now assumethat a BCell connects only one ACell, with a weight w + uniform for all BCells, so that W BA = w + I N,N , w + > . We also assume that ACells connect to BCells with a connectivity matrix W , not necessarily, symmetric,with a uniform weight − w − , w − > , so that W AB = − w − W .We note κ n , n = 1 . . . N , the eigenvalues of W ordered as | κ | ≤ | κ | ≤ · · · ≤ | κ n | and (cid:126)ψ n is thecorresponding eigenvector. We normalize (cid:126)ψ n so that (cid:126)ψ † n . (cid:126)ψ n = 1 where † is the adjoint. (Note that, as W is not symmetric in general, eigenvectors are complex). From the eigenvalues and eigenvectors of W onecan compute the eigenvalues and eigenvectors of L (see Appendix D.2), and infer stability conditions for thelinear system. The main conclusions are the following:1. The stability of the linear system is controlled by the reduced, a-dimensional parameter: µ = w − w + τ ≥ , (44)where: τ = 1 τ A − τ B , (45)with a degenerate case when τ A = τ B , considered in the appendix.2. If W is symmetric, its eigenvalues κ n are real, but the eigenvalues of L can be real or complex. To each κ n correspond actually to eigenvalues λ ± n of L (see eq. (71)).(a) If κ n < , the two corresponding eigenvalues of L are real and one of the two correspondingeigenmodes of L becomes unstable when: w − w + > − τ A τ B κ n . (46)(b) If κ n > and if τ (cid:54) = 0 the corresponding eigenvalues of L are complex conjugate if: µ > κ n ≡ µ n,c . (47) he corresponding eigenmodes are always stable.3. If W is asymmetric, eigenvalues κ n are complex, κ n = κ n,r + i κ n,i . The eigenvalues of L have theform λ β = λ β,r + i λ β,i , β = 1 . . . N with: λ β,r = − τ AB ± τ √ √ a n + u n ; λ β,i = ± τ √ √ u n − a n , (48)where a n = 1 − µ κ n,r and u n = (cid:113) ( 1 − µ κ n,r ) + 16 µ κ n,i = (cid:113) − µ κ n,r + 16 µ | κ n | . Notethat we recover the real case when κ n,i = 0 by setting u n = a n .Instability occurs if λ β,r > for some β . This gives: a n + u n > τ τ AB , (49)a condition on µ depending on κ n,r and κ n,i . Remarks:
The introduction of the a dimensional parameter µ allows us to simplify the study of the jointinfluence of w − , w + , τ on dynamics because stability is controlled by µ only. In other words, a bifurcationcondition of the form µ = µ c signifies that this bifurcation holds when the parameters w − , w + , τ lay on themanifold defined by w − w + τ = µ c .We now show this in two examples of connectivity and afferent instabilities. We consider the case where the matrix W , connecting ACells toBCells, is a matrix of nearest neighbours symmetric connections. In this case, W can be written in terms ofthe discrete Laplacian ∆ on a d dimensional regular lattice, d = 1 , with lattice spacing δ A = δ B , set hereequal to without loss of generality: W = 2 d I + ∆ . (50)Because of this relation we will often use the terminology Laplacian connectivity for the nearest-neighboursconnectivity. We also assume that dynamics holds on a square lattice with null boundary conditions. That is,ACell and BCells are located on d -dimensional grid with indices i x , i y = 0 . . . L + 1 where, the voltage andactivity of cells with indices i x = 0 , i x = L + 1 , i y = 0 or i y = L + 1 , vanish.The eigenvalues and eigenvectors are explicitly known in this case. They are parametrized by a quantumnumber n = n x ∈ { . . . L = N } in one dimension, and by two quantum numbers ( n x , n y ) ∈ { . . . L = N } in two dimensions. They define a wave vector (cid:126)k n = (cid:16) n x πL +1 , n y πL +1 (cid:17) corresponding to wave lengths (cid:16) L +1 n x , L +1 n y (cid:17) .Hence, the first eigenmode ( n x = 1 , n y = 1) corresponds to the largest space scale (scale of the whole retina)with the smallest eigenvalue (in absolute value) s (1 , = 2 (cid:16) cos (cid:16) πL +1 (cid:17) + cos (cid:16) πL +1 (cid:17) − (cid:17) . To each of theseeigenmodes is related a characteristic time τ n = λ n . The slowest mode is the mode ( 1 , . In contrast, thefastest mode is the mode ( n x = L, n y = L ) , corresponding to the smallest space scale, the scale of the latticespacing δ = 1 .Eigenvalues κ n can be positive or negative. Consider for example the dimensional case, where κ n =2 cos (cid:16) nπL +1 (cid:17) . We choose L even to avoid having a zero eigenvalue κ L . Eigenvalues κ n , n = 1 . . . L arepositive, thus the corresponding eigenvalues λ ± n of L are complex, and stable. The modes with the largestspace scale Ln are therefore stable for the linear dynamical system, with oscillations. Eigenvalues κ n , n = + 1 . . . L are negative, thus the corresponding eigenvalues λ ± n of L are real. From (46) the mode n becomesunstable when : w − w + > − τ A τ B
12 cos (cid:16) nπL +1 (cid:17) . (51)Therefore, the first mode to become unstable is the mode L with the smallest space scale (lattice spacing).For large L , this happens for w − w + ∼
12 1 τ A τ B . This instability induces spatial oscillations at the scale of thelattice spacing. When w − w + further increases the next modes become unstable. This instability results ina wave packet following the drive (as shown in Fig. 6). The width of this wave packet is controlled by theunstable modes and by non linear effects. We now illustrate the relations of these spectral properties with themechanism of anticipation. Numerical results.
In all the following 1D simulations, we consider a bar with a width µm , movingin one dimension at constant speed v = 3 mm/s . We simulate 100 BCells, 100 ACells and 100 GCells placedon a 1D horizontal grid, with a uniform spacing of δ b = δ a = δ g = 30 µm between to consecutive cells. Attime t = 0 , the first cell lies at µm to the right of the leading edge of the moving bar. We set τ B = 300 ms , τ a = 50 ms, τ A = 100 ms , corresponding to τ = 150 ms (eq. (45)). We vary the value of weights w + , w − . For the sake of simplicity, we also choose w + = − w − = w to have only one control parameter.We investigate how the bipolar anticipation time ∆ B and the maximum in the response R B depend on w .This is summarized in Fig. 7 top, where we have shown the effect of gain control alone (blue horizontalline, independent of w ), the effect of ACells lateral connectivity alone (red triangles) and the compoundeffect (white squares). Anticipation time is averaged over all cells. On the same figure (bottom) we see theresponses of two neighbour cells lying at the center of the lattice.As w increases we observe three areas of interest: the first, (A), corresponds to a regime where ACellsconnectivity has a negative effect on anticipation, competing with gain control. As w is small the anticipationis controlled by the direct pathway I, II of Fig. 1, from BCells to GCells, with a small inhibition coming fromACells, thereby decreasing the voltage of BCells and impairing the effect of gain control. This explains whythe anticipation time in the case of lateral connectivity + gain control is smaller than the anticipation timeof gain control alone. The network effect (red triangles) on anticipation time increases with w though. Thiscorresponds to the "push-pull" effect already evoked above in section 2.3. When a BCell B i feels the stimulus,its activity increases favoured by the stimulus, it increases the voltage of the connected ACell, inhibiting thenext Bcell B i +1 thereby inducing a feedback loop, the push-pull effect, enhancing the voltage of B i .In zone (B) the push-pull effect becomes more efficient than gain control alone. In this region, the voltageof the BCell feeling the bar increases fast, while the voltage of its neighbours becomes more and more neg-ative, enhancing the feedback loop. This holds until the voltage rectification (5) takes place. This is the timewhen the dynamical system gets out of Ω . The push pull effect then saturates and V B i reaches a maximum,corresponding to a peak in activity. This peak is reached faster than the peak in the function G B ( A ) . Thus,the peak of R B i ( t ) occurs at the same time as the peak of N B ( V B i ( t )) , and, thus, before the reference peak(time t B for isolated BCells defined in section 3.1.1). In other words, the ACells lateral connectivity allowsthe BCell to outperform the gain control mechanism for anticipation. As w increases in zone B the push-pulleffect (averaged over BCells) reaches a maximum then decreases. This is because the increase in w makes theinhibitory effect of ACells stronger and stronger on silent Bcells which then remain silent longer and longerbecause the ACells voltage increases with w , and it takes longer for it to decrease and de-inhibit the neigh-bours. The silent cells are less and less sensitive to the stimulus, being strongly and durably inhibited.In region C, the anticipation is again dominated by gain control. In this case, the effect on cells depends onthe parity of their index. The response of BCells is either completely suppressed or identical to the responseof the reference case (with gain control alone). This is why the average anticipation time with gain control isabout half of the gain control without network effect. Cells that are inhibited do no participate to anticipation,and the others anticipate in the same way than with gain control alone. Note that this "parity" effect is due tothe nearest neighbours connectivity and the symmetry of interactions. e now interpret and complete these results from the point of view of the spectrum of L and associateddynamics. The fastest mode to destabilize, corresponds to the smallest space scale, the lattice spacing. Thisis a mode with alternate sign, at the scale of the lattice. We call it the "push-pull" mode, as it is precisely whatmakes the push-pull effect. When the push-pull mode becomes unstable, the excited BCell becomes more andmore excited and the next BCell more and more inhibited. However, the time it takes, τ L , has to be comparedto the time where the bar stays in the RF, τ bar (and more generally the time it takes to RF kernel to respondto the bar). In the case of the simulation σ center = 90 µm (see Appendix, table (A)) and v = 3 µm/ms givinga characteristic time τ bar = 270 ms , whereas, as we observed τ L < ms. The push-pull mode is thereforequite faster than τ bar so the push-pull effect takes place fast and lead to a fast exponential increase of the frontdepicted in Fig. 6 right. This explains the rapid increase of network anticipation effect observed in regions A,B of Fig. 7. In this section, we study the behaviour of the model using the more realistic, probabilistic type of con-nectivity presented in section 2.3.3 and more thoroughly studied in Appendix C. Within this framework, agiven ACell A i receives the upstream activity from the BCell lying at the same position, B i , with a constantweight w . The same ACell inhibits BCells with which it is coupled through the random adjacency matrix W ,generated by the probabilistic model of connectivity, and the weight matrix W BA = − w W . We recall that theconnectivity depends on a scale parameter ξ for the branch length), and the mean and variance ¯ n, σ for thedistribution of the number of branches. These parameters can be found in the table (A) in appendix. Eigenmodes of the linear regime.
Similarly to section 3.2.2 we now analyse the spectrum of L when W is a random connectivity matrix. Although a couple of results can be established (using the Perron-Frobeniustheorem) we have not been able to find general mathematical results on the spectrum or eigenvectors of thisfamily of random matrices. We thus performed numerical simulations.The spectrum of L is deduced from the spectrum of W as exposed above. The spectrum of W dependson ¯ n, σ and ξ . In Fig. 8 we have plotted, on the left, an example of such spectrum. This is the spectral density(distributions of eigenvalues in the complex plane) obtained from the diagonalization of matrices × (so the statistics is made over eigenvalues). We note that the largest eigenvalues is always realpositive, a straightforward consequence of Perron-Frobenius theorem [32, 62]. More generally, we observean over-density of real eigenvalues. The same holds for random Gaussian matrices with independent entries N (0 , N ) [27] whose asymptotic density converges to the circular law [34]. The shape of the spectral densityin our model differs from the circular law though and it depends on the parameters ¯ n, σ and ξ .On the same figure we show the corresponding spectral density of L obtained from eq. (48) for w =0 . , . , . We have taken here τ A = 30 , τ B = 10 ms to see better the transitions with w (level lines in Fig.9). There is an evident symmetry with respect to τ AB = − . expected from the mathematical analysis. Wesee that the largest eigenvalue is real (although it is not necessarily related to the largest eigenvalue of W ).We also see that, as w increases, a large number of (complex) eigenvalues become unstable. There is actuallya frontier of instability that we have plotted in the plane w, ξ for different values of ¯ n . This is shown in Fig.9 (dashed line). The level line is the frontier of instability of the linear dynamical system. This frontier hasthe (empirical) form ( ξ − ξ ) . ( w − w
0) = c where c has the dimension of a characteristic speed.What matters here is that there are complex unstable eigenvalues with no specific resonance relationsbetween them. They are therefore prone to generate destructive interferences in (37). Numerical results
In fig. 10 we consider, similarly to Fig. 7 for Laplacian connectivity, the effect ofrandom connectivity on anticipation, compared to pure gain control mechanism. In contrast to the Lapla-cian case, we have here more parameters to handle: ξ , which controls the characteristic length of branchesand ¯ n, σ n which control the number of branches distribution. We present here a few results where ξ varies Anticipation in the Laplacian (nearest-neighbours) case. Top.
Anticipation time andmaximum bipolar response as a function of the connectivity weight w . The blue line correspondsto gain control alone (it does not depend on w ). Red triangles corresponds to the effect of lateralACells connectivity, without gain control. White squares correspond to the compound effect ofACells connectivity and gain control. The regimes A, B, C are commented in the text. Bottom.
Response curves of ACells and BCells corresponding to the three regimes : (A) w = 0 . ms − with a small cross-inhibition, (B) w = 0 . ms − with an opposition in activity between the blue(cell 50) and red cell (51), (C) w = 0 . ms − , where the red cell (51) is completely inhibited by cell50. 25igure 8: Spectral density of eigenvalues for ξ = 2 , ¯ n = 4 , σ n = 1 . Top Left. For the matrix W (density estimated over samples). Density is represented in color plots, in log scale. Thecolor bar refers to powers of . Top Right.
Spectral density of L for w = 0 . . Bottom Left.
Spectral density of L for w = 0 . . Bottom Right.
Spectral density of L for w = 0 . . Unstableeigenvalues are on the right to the vertical, dashed, line x = 0 Heat map for the largest real part eigenvalue in the plane w, ξ , for different values of ¯ n . Left. ¯ n = 1 . Right. ¯ n = 4 . Colour lines are level lines. The level line is the frontier of instabilityof the linear dynamical system. whereas the average number of branches ¯ n = 2 ( σ n = 1 ). A more systematic study is done in [67]. The inter-est of varying ξ is to start from a situation which is close to the Laplacian case (characteristic distance ξ = 1 )and to increase ξ to see how the size of the dendritic tree of ACells may impact anticipation. This is a prelimi-nary step toward considering different physiological ACells type (e.g. narrow-, medium-, or wide-field [26]).Note however that the probability of connection given the distance of cells (fixed by ¯ n, σ n ) implicitly impacts w and the anticipation effects.The main difference with the Laplacian case is the asymmetry of connections. Here, symmetry meansthat if Acell j connects the BCell i , then the Acell i connects the BCell j too. This does not necessarily holdfor random connectivity and this has a strong impact on the push-pull effect and anticipation. So, even if theconnectivity is short-range when ξ is small, mainly connecting nearest neighbours, we observe already a bigdifference with the Laplacian case. This is shown in Figure 10, where ξ = 1 . Similarly to the Laplacian case weobserve 3 main regions depending on w . To have the same representation as Fig. 7 we present V A , N B , R B for two connected cells (here Acell and BCell ). However, in this case, connection is not symmetric:ACell inhibits BCell but ACell does not inhibits BCell .We observe regimes, as in the Laplacian case. In the first region (A) ACells random connectivity has anegative effect on anticipation, as compared to gain control alone. However, since in this case the "push-pull"effect is not evoked, this decay simply comes from the fact that BCell receives an inhibition for the ACell , that reduces the effect of gain control. This inhibition is though not strong enough to significantly shiftthe peak response, as in region (B).Indeed, in region (B), the inhibition of BCell is strong enough to outperform the effect of gain control.In this case, and similarly to the Laplacian case, the peak of R B i ( t ) occurs at the same time as the peak of N B ( V B i ( t )) , and, before the reference peak. However this effect is not consistent over all cells and only occursfor BCells that receive active inhibition. This explains why the performance of the Laplacian connectivity isbetter, on average, in this region.Finally, as w grows higher, the inhibition grows stronger, completely inhibiting BCell . Cells that do notreceive any inhibition, as BCell in this example, keep a response that is identical to the response without Average anticipation in the random connectivity case . Top.
Bipolar anticipation timeand maximum in the response R B as a function of the connectivity weight w , in the case of arandom connectivity graph with ξ = 1 , ¯ n = 2 and σ n = 1 . Bottom.
Response curves of ACellsand BCells − corresponding to the three regimes : (A) w = 0 . ms − , (B) w = 0 . ms − , (C) w = 0 . ms − . 28 Cell connectivity. The fraction of cells receiving inhibition in this case being quite small (about 15), thisexplains why the stationary value of anticipation is fairly close to the value with gain control alone.
The role of the characteristic distance
In figure 11, we analyse the effect of the characteristic length ξ on anticipation. On the top of the figure we represent the joint effect of the random ACell connectivity andgain control on anticipation for three values of ξ . At the bottom we represent the only effect of the randomACell connectivity for the same values of ξ . We observe that performance in anticipation decreases with ξ .More precisely, we observe an anticipatory effect in this case, as shown in figure 11 bottom, but this effect isnot able to compete with gain control alone. Even worse, the compound effect shown in 11 top is disastroussince increasing w renders the anticipation time smaller and smaller.This spurious effect can be interpreted through the analysis made in section 3.2.1, eq. (37). From thespectrum of L , we see that there are unstable complex eigenvalues whose number increases with w . Theseeigenvalues are prone to generate destructive interferences, especially when their number becomes large as w increases, explaining the small peak in region B . The consequence on cells activity and gain control canbe dramatic as seen in the red trace of Fig. 10 B bottom, line R B . This depends on the precise connectivitypattern when long range connections from ACells to BCells induce a desensitization of BCells, which is notcounterbalanced by the push-pull effect as in the Laplacian connectivity case. The two numerical examples considered in this section emphasize the role of symmetry in the synapses,and more, generally the role of complex versus real eigenvalues in the spectrum of L . Recall that, from section3.2.1, if W is symmetric complex eigenvalues are always stable, so, for the type of architecture consideredhere, unstable destructive interferences only occur when W is asymmetric. This leads to several questions,potential subjects for further studies.1. How much does anticipation depend on the degree of asymmetry in the matrix W ? The way wegenerate the random connectivity in the model does not allow us to tune the degree of asymmetry(i.e. the probability that a connection A j → B i exists simultaneously with a connection A i → B j ).Therefore, one has to find a different way to generate the connectivity. From the mathematical analysismade in Appendix C a distribution depending exponentially on the distance, with a tunable probabilityto have a symmetric connection, could be appropriate. We don’t know about any experimental resultscharacterizing this degree of symmetry of the connections in the retina. On mathematical grounds,and from the analogy of the spectrum of W with a circular law, one could expect the spectrum of W to become more and more elongated on the real axis as the degree of symmetry increases, in an ellipticlike law [48].2. Non linear effects.
The destructive interference effect in our model is partly due to the linear nature ofthe ACells dynamics. In non linear dynamics, eigenvalues of the evolution operator can display reso-nances conditions favouring constructive interferences. On biological grounds, it is for example knownthat Starburst Amacrine Cells display periodic bursting activity during development, disappearing afew days after birth [86]. Bursting and its disappearance can be understood in the framework of bifur-cation theory of a non linear dynamical system featuring these cells [43]. In this setting, even if theyare not bursting in the mature stage, SACs remain sensitive to specific stimulation that can temporallysynchronize them, thereby enhancing the network effect, with a potential effect on anticipation.
In this section, we study the network ability to improve anticipation in the presence of gap junctionscoupling, as in eq. (23), and gain control at the level of GCells.We start first with mathematical results and show then simulation results.
Role of the characteristic branch length ξ on anticipation. Top. The joint effect of therandom ACells connectivity and gain control on anticipation for ξ = 1 , , . Left.
Average bipolaranticipation time;
Right.
Maximum value of the bipolar response R B . Bottom.
The single effectof the random ACells connectivity on anticipation with the same representation.
We use a continuous space limit for a one dimensional lattice. The extension to dimension is straight-forward. Here, x corresponds to the preferred direction of the direction sensitive cells. We consider acontinuous spatio-temporal field V ( x, t ) , x ∈ R , such that V G k ≡ V ( kδ G , t ) . We assume likewise that V ( P ) G k ≡ V ( P ) ( kδ G , t ) for some continuous function V ( P ) ( x, t ) corresponding to the GCells bipolar poolinginput (17) and we take the limit δ G → . In this limit eq. (23) becomes: ∂V G ∂t = f ( x, t ) − v gap ∂V G ∂x + O ( δ G ) , (52)where v gap ≡ w gap δ G has the dimension of a speed and ∂V ( P ) ( x,t ) ∂t ≡ f ( x, t ) . Finally, we note C ( x ) the initialprofile so that V ( x, t ) = C ( x ) . olution. Neglecting terms of order δ G the general solution of (52) is : V G ( x, t ) = C ( x − v gap ( t − t )) + (cid:90) tt f ( x − v gap ( t − u ) , u ) du. Eq. (52) is a transport equation of ballistic type [67]. For example, if we consider a stimulation of the form V ( P ) ( x, t ) = h ( x − v t ) , where h is a Gaussian pulse of the form (3), propagating with speed v , and an initialprofile C ( x ) = h ( x − v t ) , the voltage of GCells obeys: V G ( x, t ) = vv − v gap h ( x − vt ) (cid:124) (cid:123)(cid:122) (cid:125) π stim − v gap v − v gap h ( x − v gap t − ( v − v gap ) t ) (cid:124) (cid:123)(cid:122) (cid:125) π gap . (53)When v gap = 0 the GCells voltage follows the stimulation i.e. V G ( x, t ) = h ( x − vt ) . In the presence ofgap junctions there are two pulses: the first one, π stim with amplitude vv − v gap propagating at speed v andfollowing the stimulation; the second one, π gap , with amplitude − v gap v − v gap , propagating at speed v gap .We have the following cases (we take t = 0 for simplicity). An illustration is given in Fig. 12.1. If v and v gap have the same sign:(A) If v gap < v , the front π stim is amplified by a factor vv − v gap , whereas there is a refractory front π gap ,proportional to v gap , behind the excitatory pulse.(B) If v = v gap , V G ( x, t ) = h ( x − v gap t ) + v gap ( t − t ) h (cid:48) ( x − v gap t ) which diverges like t when t → ∞ and x → + ∞ . This divergence is a consequence of the limit δ G → in (52).(C) If v gap > v the amplitude of π stim follows the stimulation with a negative sign (hyper polarization)whereas π gap is ahead of the stimulation, with a positive sign, travelling at speed v gap .2. If v and v gap have the opposite sign, we set v = − α v gap , with α > . Then, the front π stim follows thestimulus but is attenuated by a factor α α . The front π gap propagates in the opposite direction with anattenuated amplitude α .This shows that these gap junctions favour the response to motion in the preferred direction and attenuatethe motion in the opposite direction although the attenuation is weak. The effect is reinforced by gain control[67]. The most interesting case is 1 c where these gap junctions can induce a wave of activation ahead of thestimulation. Effect of gain control.
When the low voltage threshold N G (19) and the gain control G G ( A ) (21) are ap-plied to V G ( x, t ) there are two effects: (i) the hyperpolarized front is cut by N G ; (ii) the positive pulse inducesa raise in activity, which, in turn, triggers the ganglion gain control G G ( A ) inducing an anticipated peak inthe response of the GCell, similar to what happens with BCells, with a different form for the GCell gain con-trol though. Moreover, in contrast to pathway II of Fig. 1 where only gain control generates anticipation, inpathway IV the wave of activity generated by gap junctions increases anticipation by two distinct effects. If v gap < v the cell’s response propagates at the same speed as the stimulus, but its amplitude is larger thanthe case with no gap junction (term π stim ). From eq. (53) this results in an increase of h B to an effectivevalue h B vv − v gap inducing an improvement in the anticipation time (with a saturation of the effect, though,as v gap → v ). If v gap > v the cell’s response propagates at a larger speed than the stimulus (term π gap ),so that the cell responds before the time of response without gaps. This induces as well an increase in theanticipation time. .3.2 Numerical illustrations We consider a bar with a width µm , moving in one dimension at constant speed v = 3 mm/s . Wesimulate here 100 GCells, placed on a 1D horizontal grid, with a spacing of µm between to consecutivecells. At time t = 0 , the first cell lies at µm from the leading edge of the moving bar.We investigate how the GCells anticipation time and GCells firing rate depend on v gap in Fig. 12. The topshows the effect of gain control alone (blue horizontal line, independent of v gap ), the effect of the asymmetricgap junction connectivity alone (red triangles) and the compound effect (white squares). Anticipation time isaveraged over all GCells. On the bottom part of the figure, we show the responses of two GCells of indices30 and 60, spaced by µm .As explained in the section 3.3.1, we observe the 3 regimes A,B,C mathematically anticipated above. Notethat, for these parameter values, the negative trailing front predicted in A is not visible. The asymmetry observed by Trenhlom et al. is due to the specific structure of the direction selectiveGCell dendritic tree [73]. However, in general gap junctions connectivity is expected to be symmetric. So, tobe complete we consider here the effect of symmetric gap junctions on anticipation. It is not difficult to derivethe equivalent of eq. (52) in this case too. This is a diffusion equation of the form: ∂V G ∂t = f ( x, t ) + D gap ∆ V G + O ( δ G ) , where D gap = w gap δG is the diffusion coefficient and ∆ is the Laplacian operator.The response to a Gaussian stimulus of the form (3) reads: V G ( x, y, t ) = (cid:104) H x,y,t ∗ f (cid:105) , (54)where: H ( x, y, t ) = e − x y Dgap t πD gap t (55)which is the heat equation diffusion kernel.Recall that f ( x, t ) ≡ ∂V ( P ) ( x,t ) ∂t . So, if V ( P ) ( x, t ) = h ( x − v t ) , where h is a Gaussian pulse of the form (3)propagating with speed v , f is a bimodal function of the form vσ h ( x − v t ) × ( x − v t ) , the shape of which canbe seen in Fig. 13 bottom, second row. The convolution with the heat kernel leads to a front propagating atthe same rate as the stimulus, with a diffusive spreading whose rate is controlled by D gap . In particular, thereis positive bump ahead of the motion, which can induce anticipation, as shown in Fig. 13 top. The effect isweak, though, essentially because the diffusive spreading makes the amplitude of the response decrease fastas a function of D gap .Although this positive front, for small D gap , increases a bit the anticipation time by accelerating the gaincontrol triggering, rapidly the peak in the response R B is lead by the voltage peak corresponding to the pos-itive bump, with a low voltage. The position of this peak is, roughly, at a distance σ = (cid:112) σ center + σ B fromthe peak of the Gaussian pool, where σ center is the width of the center RF and σ B the width of the bar. Thiscorresponds to a time σv ahead of the peak in the drive, fixing a maximal value to the anticipation time (seethe saturation of the anticipation time curve in Fig. 13 top, left). In our case, σ ∼ given a saturation peakat µm µm/ms = 44 . ms . A consequence of the voltage decay is the corresponding power law ( √ D gap for large D gap ) decay of the firing rate (Fig. 13 top, right).To conclude, the situation with symmetric gap junctions is in high contrast with direction selective gap junc-tions where the response to stimuli was ballistic and was not decreasing with time. On this basis we considerthat, for symmetric gap junctions, the anticipation effect is irrelevant, especially taking into account the small-ness of the voltage response in case C. Anticipation for non symmetric gap junctions. Top.
GCell anticipation time andmaximum firing rate as a function of the gap junction velocity v gap . Bottom : response curves ofGCells corresponding to the three regimes : (A) v gap = 0 . mm/s , (B) v gap = 3 mm/s , (C) v gap =12 mm/s . The curves display 3 main regimes (see text): In (A) v gap < v and the positive frontpropagates at the same speed as the pooling voltage triggered by the stimulus; In (B), v gap = v ,the positive front and the negative fronts both propagate at the speed v gap and the amplitude ofthe positive front ( V G ( t ) ) increases with t ; In (C), v gap > v and the positive front propagates fasterthan the stimulus so that the peak of activity arises earlier. The negative front propagates at thestimulus speed. 33igure 13: Anticipation for symmetric gap junctions. Top.
GCell anticipation time and maximumfiring rate as a function of the gap junction velocity v gap . Bottom. response curves of GCellscorresponding for three values of v gap : (A) v gap = 0 . mm/s , (B) v gap = 3 mm/s , (C) v gap = 12 mm/s . For consistency, we have kept the same values as in the asymmetric case. Here, anticipationtime grows continuously until saturation, while the maximum firing decreases like a power lawas a function of v gap . We investigate in this section how the GCell anticipation time and GCells firing rate depend on the gapjunction conductance in the case of symmetric gap junctions. In figure 13 top, we use the same representationas Fig. 12. For consistency with the direction sensitive case, we choose v gap = D gap δG as control parameter. Wealso take (A) v gap = 0 . mm/s , (B) v gap = 3 mm/s , (C) v gap = 12 mm/s in figure 13 bottom. This correspondsto a diffusion coefficient (A) D gap = 18 × − mm /s , (B) D gap = 90 × − mm /s , (C) D gap = 360 × − mm /s In this section we have shown how gap junctions direction sensitive cells can display anticipation dueto the propagation of a wave of activity ahead of the stimulus. This effect is negligible for symmetric gapjunctions. Note that symmetric gap junctions are known to favour waves propagation, for example in the arly development (stage I, see [42] and reference therein for a recent numerical investigation). Here gapjunctions are considered in a different context, due to the presence of a non stationary stimulus, triggeringthe wave.Let us now comment our computational result. How does it fit to biological reality ? Depending on thegap-junction conductance value the propagation patterns we predict are quite different.What is the typical value of v gap in biology ? It is difficult to make an estimate from the expression v gap = g gap δ G C . The membrane capacity C and gap junctions conductance can be obtained from the literature(for connexins Cx36, g gap ∼ − pS [68]) but the distance δ G is more difficult to evaluate. In the model, thisis the average distance between GCells’ soma which corresponds to ∼ − µm . But, in the computationwith gap junctions what matters is the length of a connexin channel which is quite smaller. Taking δ G as thedistance between GCells assumes a propagation speed between somas at the speed of a connexin, which iswrong because most of the speed is constrained by the propagation of action potential along the dendritictree. So we used a phenomenological argument (we thank O. Marre to point out it to us). The correlationof spiking activity between GCells neighbours is about − ms for cells separated by ∼ − µm [79].This gives a speed v gap in the interval [ 40 ,
150 ] mm/s which is quite fast compared to the bar speed inexperiments.So we are in case 1 b and one should observe an experimental effect ? To the best of our knowledgean effect of DSGC gap junctions on motion anticipation has not been observed. But we don’t know aboutexperiments targeting precisely this effect. It would be interesting to block gap junctions and address Berryet al. [8] or Chen et al. [17] experiments in this case. The difficulty is that blocking gap junctions blocks manyessential retinal pathways. We do not pursue this discussion further concluding that our model proposes acomputational prediction that could be interesting to be experimentally investigated.
In this section, we present some examples of retinal responses and anticipation to trajectories more com-plex than a bar moving in one dimension with a uniform speed. The aim here is not to do an exhaustivestudy but, instead, to assess qualitatively some anticipatory effects not considered in the previous sections.
The flash lag effect is an optical illusion where a bar moving along a smooth trajectory and a flashed barare presented to the subject, and are perceived with a spatial displacement, while they are actually aligned.A variation of this illusion consists of a bar moving in rotation, a bar flashed in angular alignment, givingrise to a perceived angular discrepancy. We have investigated this effect in our model, in the presence of thedifferent anticipatory effects considered in the paper.Fig. 14 shows the response to a bar moving with a smooth motion, while a second bar is flashed inalignment with the first bar at one time frame. The first line shows the stimuli, consisting of 130 frames, of abar moving at 2.7 mm/s, with a refreshment rate of 100 Hz. The second line shows the GCell response withgain control, the third line presents the effect of lateral amacrine connectivity, in the case of a Laplacian graph.Keeping the same values of parameters as in the 1D case, we set w = 0 . ms − corresponding to the case Bin Fig. 7. Finally, the last line shows the effect of asymmetric gap junctions, having a preferred orientation inthe direction of motion, with v gap = 9 mm/s.In the case of the gain control response, the peak of response to the moving bar is shifted by about msin the direction of motion, as compared to the static bar. The flashed bar elicits a lower response, given itsvery short appearance in comparison with the characteristic time of adaptation. We choose this time shortenough to avoid gain control triggering, explaining the difference with the strong response observed by Chenet al. [17] in the presence of a still bar.In the case of amacrine connectivity, the moving bar representation is shrunk as compared to the gaincontrol case, given the prevalence of inhibition, while the level of activity for the flashed bar remains roughly he same. In this case, cells responding to the moving bar reach their peak activity slightly earlier (about ms for these parameters value) than in the gain control case (Fig. 14, (B) top.)Finally, asymmetric gap junction connectivity displays a wave propagating ahead of the bar, increasingthe central blob, which is much larger than the size of the bar in the stimulus, while the flashed bar activityremains similar to the previous cases. Figure 14:
Flash lag effect with different anticipatory mechanisms.
A) Response to a flash lagstimulus : a bar moving in smooth motion, with a second bar flashed in alignment with the first barfor one time frame. The first line shows the stimulus, the second line shows the GCells responsewith gain control, the third line presents the effect of lateral ACells Laplacian connectivity with w = 0 . ms − , and the last line shows the effect of asymmetric gap junctions with v gap = 9 mm/s. B) Time course response of (top) a cell responding to the flashed bar, and (bottom) a cellresponding to the moving bar. Dashed lines indicate the peak of each curve. In this subsection, we assess the effect of the three anticipatory mechanisms on a parabolic trajectory.The interest is to have a trajectory with a change in direction and speed, thus an acceleration. The stimulusconsists of frames, displayed at 10 Hz. The simulations parameters and connectivity weights are the sameas the ones used in the previous section. Fig. 15 shows the response to a dot moving along a parabolictrajectory.In the case of gain control, GCells response is more elongated, which has a distortion effect on the dotrepresentation near the turning point of the trajectory ( − ms). Cells responding near the trajectoryturning point are still anticipating motion, as the peak response of the gain control curve is slightly shifted tothe left, compared to the RF response (Fig. 14, B).In the case of amacrine connectivity, the elicited response is also more localized, as compared to the gaincontrol response, and the flow of activity follows more accurately the stimulus. This is a direct consequenceof the sensitivity of the ACell connectivity model to the stimulus acceleration. In this case, the peak responseis also more shifted as compared to the gain control case.Finally, the gap junction connectivity model performs worse in this case, giving rise to a propagatingwave that doesn’t follow the trajectory, since the latter is not parallel to the direction to which GCells are ensitive. Cells responding near the trajectory turning point have a higher level of activity and an increasedlatency, while the peak response roughly corresponds to the gain control case. Figure 15:
Effect of anticipatory mechanisms on a parabolic trajectory.
A) Response to a dotmoving along a parabolic trajectory. The first line shows the stimulus, the second line shows theGCells response with gain control, the third line presents the effect of lateral ACells connectivitywith w = 0 . ms − , and the last line shows the effect of asymmetric gap junctions with v gap = 9 mm/s . B) Time course response of a cell responding to the dot near the trajectory turning point.Linear response corresponds to the response to the stimulus without gain control. We investigate in this subsection a two dimensional example of motion where angular anticipation takesplace. The stimulus consists here of 72 frames, displayed at 100 Hz, of a bar moving at a constant angularspeed of . rad/ms.Fig. 16 shows the retina response to a rotating bar, with the angular orientation of activity as a functionof time for the different models. We used Matlab to estimate the bar orientation from the displayed activity,fitting the set of activated points by an ellipse whose principal axis determines the response orientation.In the three cases, one can see that the response around the center of the bar is suppressed due to gaincontrol adaptation. While the gain control activity orientation roughly follows the linear response ( Fig.16 (B)), the ACells response shows a slight angular shift (frames : 250-300 ms) which is also visible on theresponse orientation time course. The ACells angular anticipation is however only observed during the firstperiod of the bar. Interestingly, this effect vanishes during the second rotation, due to a persistent effect of theactivation function, generating a sort of a suppressive effect erasing the second occurrence of the bar (frame: 450 ms). We shall point out that the ACells connectivity weight in this simulation has been reduced to w = 0 . ms − , since with a value of w = 0 . ms − used in the previous simulations, the response to thesecond rotation of the bar is completely suppressed.Finally, similarly to the parabolic trajectory, the gap junction connectivity model performs worse due tothe wave propagating from left to right, distorting once more the bar shape. Consequently, the bar activityorientation in this case has been discarded in Fig. 16 B, the orientation estimate giving poor results. Anticipation for a rotating bar.
A) Response to a bar rotating at . rad/ms . The firstline shows the stimulus, the second line shows the GCells response with gain control, the thirdline presents the effect of lateral ACells connectivity with w = 0 . ms − , and the last line showsthe effect of asymmetric gap junctions with v gap = 9 mm/s . B) Time course response of the barorientation in the reconstructed retinal representations. This section shows how lateral connectivity can play a role in motion anticipation of 2D stimuli, both inthe case of the classical flash lag effect, and more complex trajectories. Indeed, for a given network setting,ACells connectivity can noticeably improve anticipation with respect to Gain Control, in all three stimuli, andhas also the advantage of being sensitive to trajectory shifts (sec. 3.4.2).While gap junction connectivity improves anticipation when the trajectory of the bar is parallel to thepreferred GCells direction, it also induces more blur around the bar, and shape distortion in the case ofparabolic motion and rotation, suggesting a trade-off between anticipation and object recognition for thisspecific model.
Using a simplified model, mathematically analysed with numerical simulations examples, we have beenable to give strong evidences that lateral connectivity - inhibition with ACells, gap junctions - could partici-pate to motion anticipation in the retina. The main argument is that a moving stimulus can - under specificconditions mathematically controlled - induce a wave of activity which propagates ahead of the stimulusthanks to lateral connectivity. This suggests that, in addition to local gain control mechanism inducing ananticipated peak of GCells activity, lateral connectivity could induce a mechanism of neural latencies reduc-tion, similar to what is observed in the cortex [7, 70, 39]. This is visible in particular in Fig. 14, where the gapjunction coupling induces a wave which increases the GCell level activity before the bar reaches its RF.Yet, this studies raise several questions and remarks. The first one is, of course, the biological plausibil-ity. At the core of the model, what makes the mathematical analysis tractable is the fact that we can reducedynamics, in some region of the phase space, to a linear dynamical system. This structure is afforded bytwo facts: (i) Synapses are characterized by a simple convolution; (ii) Cells, especially Acells, have a simplepassive dynamics, where non linear effects induced e.g. by ionic channels are neglected, as well as prop- gation delays. As stated in the introduction the goal here is not to be biologically realistic, but instead,to illustrate potential general spatio-temporal response mechanisms taking into account specificities of theretina, as compared to e.g. the cortex. Essentially, most neurons are not spiking (except GCells and sometype of BCells or ACells, not considered here [2]). Yet, synapses follow the same biophysics as their corticalcounterpart. As it is standard to model the whole chain of biophysical machinery triggering a post-synapticpotential upon a sharp increase of the pre-synaptic voltage by a convolution kernel [24], we adopted here thesame approach. Note that it is absolutely not required, in this convolutional approach, for the pre-synapticincrease in voltage to be a spike; it can be a smooth variation of the voltage. Note also that higher orderconvolution kernels can be considered, integrating more details of the biological machinery. These higherorder kernels are represented by higher order linear differential equations [30]. Concerning point (ii) - nonlinear effects are neglected - especially for ACells (BCells have gain control), there are not so many availablemodels of ACells. A linear model for predictive coding using linear ACells has been used by Hosoya et al[36]. We discuss it in more detail below. The non linear models of ACells we know has been developed tostudy the retina in its early stage (retinal waves) and feature either AII ACells [18] or Starburst ACells [43].In section 3.2.4 we have briefly commented how non-linear mechanisms could enhance resonance effects inthe network and, thereby, favour the propagation of a lateral wave of activity induced by a moving stimulus.This would of course deserve more detailed study. Another potentially interesting non linear mechanism isshort term plasticity, discussed below.The second question one may ask about the model is about the robustness of this mechanism with respectto parameters. The model contains many parameters, some of them (BCells, GCells and gain control) comingfrom the previous paper from Berry et al [8] and Chen et al [17]. Although they didn’t perform a structuralstability analysis of their model (i.e. stability of the model with respect to small variations of parameters), webelieve that they are tuned away from bifurcation points so that slight changes in their (isolated cells) modelparameters would not induce big changes. As we have shown the situation changes dramatically when cellsare connected via lateral connectivity. Here, many types of dynamical behaviour can be expected simply bychanging the connectivity patterns in the case of ACells. A more detailed analysis would require a closer in-vestigation of ACells to BCells connectivity and an estimation of synaptic coupling, implying to define morespecifically the type of ACell (AII, Starburst, A17, wide field, medium field, narrow field, ...) and the type offunctional circuit one wants to consider. Note that ACells are difficult to access experimentally due to theirlocation inside the retina. Even more difficult is a measurement of ACells connectivity, especially the degreeof symmetry discussed in our paper. Such studies can be performed at the computational level, though,where the mathematical framework proposed here can be applied and extended. Computational results doesnot tell us what is reality but shed light on what it could be.We would like now to address several possible extensions of this work. The retino-thalamico-cortical pathway.
The retina is only the early stage of the visual system. Visualresponses are then processed via the thalamus and the cortex. As exposed in the introduction, anticipationis also observed in V1 with a different modality than in the retina. In this paper, our main focus was onthe shift of the peak response, while when studying anticipation in the cortex, the main focus lies in theincrease of the response latency, i.e. the delay between the time the bar reaches the receptive field of a corticalcolumn, and the effective time its activity starts rising. How do these two effects combine ? How does retinalanticipation impact cortical anticipation ? To answer these questions at a computational level one would needto propose a model of the retino-thalamico-cortical pathway which, to the best of our knowledge, has neverbeen done. Yet, we have developed a retino-cortical (V1) model - thus, short-cutting the thalamus - based ona mean-field model of the V1 cortex, developed earlier by the groups of F. Chavane and A. Destexhe [85][16],able to reproduce V1 anticipation as observed in VSDI imaging. The aim of this work is to understand,computationally, the effect of retinal anticipation on the cortical one, and more generally the combined effectsof motion extrapolation in the retina and V1. This is the object of a forthcoming paper (S. Souihel, M. di Volo,S. Chemla, A. Destexhe, F. Chavane and B. Cessac., in preparation). See [67] for preliminary results. etinal circuits. BCells, ACells and GCells are organized into multiple, local, functional circuits withspecific connectivity patterns and dynamics in response to stimuli. Each circuit is related to a specific task,such as light intensity or contrast adaptation, motion detection, orientation, motion direction and so on.Here, we have considered a circuit allowing the retina to detect a moving object on a moving background,where motion sensitive retinal cells remain silent under global motion of the visual scene, but fire whenthe image patch in their receptive field moves differently from the background. From our study we haveemitted the hypothesis that this circuit, spread over the retina, could improve motion anticipation thanksto what we have called the "push-pull" effect. Yet, other circuits could be studied in their role to processmotion and anticipation. We especially think of the ON-OFF cone and rod-cone pathways responsible for theseparation of highlights and shadows, allowing to provide information to the GCells concerning brighter thanbackground stimuli (ON-center) or darker than background stimuli (OFF-center) [47]. This circuit involvesboth gap junctions and ACells (AII) connectivity and our model could allow to study its dynamics in thepresence of a moving object.
Adaptation effects.
In a paper from 2005, Hosoya et al [36] have studied dynamic predictive coding inthe retina and shown how spatio-temporal receptive fields of retinal GCells change after a few seconds in anew environment, allowing the retina to adjust its processing dynamically when encountering changes in itsvisual environment. They have shown that an Amacrine network model with plastic synapses can accountfor the large variety of observed adaptations. They feature a linear network model of ACells, similar to ours,with, in addition, anti-Hebbian plasticity. Their mathematical analysis, based on linear algebra, allows todetermine the behaviour of the model in terms of eigenvalues and eigenvectors. However, their analysis doesnot carry out to the gain control introduced by Berry et al, which, as we show, renders quite more complex thespectral analysis. It would therefore be interesting to explore how plasticity in the ACells synaptic network,conjugated with local gain control, contributes to anticipation.
Correlations.
The trajectory of a moving object - which is, in general, quite more complex than a movingbar with constant speed - involves long-range correlations in space and in time. Local information about thismotion is encoded by retinal GCells. Decoders based on the firing rates of these cells can extract some of themotion features [55, 59, 22, 60, 69, 36, 44]. Yet lateral connectivity plays a central role in motion processing (seee.g. [35]). One may expect it to induce spatial and temporal correlations in spiking activity, as an echo, a trace,of the object’s trajectory. These correlations cannot be read in the variations of firing rate; they also cannotbe read in synchronous pairwise correlations as the propagation of information due to lateral connectivitynecessarily involves delays . This example raises the question about what information can be extracted fromspatio-temporal correlations in a network of connected neurons submitted to a transient stimulus. What isthe effect of the stimulus on these correlations? How can one handle this information from data where onehas to measure transient correlations ? This question has been addressed in [14]. The potential impact of thesespatio-temporal correlations on decoding and anticipate a trajectory will be the object of further studies.
Orientation selective cells.
Our model affords the possibility to consider BCells with orientation sen-sitive RF. The potential role of such BCells for predictive coding has been outlined by Johnston et al [41].In their model individual Gcells receive excitatory BCell inputs tuned to different orientations, generating adynamic predictive code, while feed-forward inhibition generates a high-pass filter that only transmits theinitial activation of these inputs, removing redundancy. Should such circuits play a role in motion anticipa-tion ? We didn’t elaborate on this in the present paper, leaving it to a potential forthcoming work. Another,important question, is "how to model a retinal network with cells having different orientation selectivity ?" AV1 cortical model has been proposed by Baspinar et al. [6] for the generation of orientation preference maps,considering both orientation and scale features. Each point (cortical column) is characterized by intrinsicvariables, orientation and scale and the corresponding RF is a rotated Gabor function. The visual stimulusis lifted in a 4-dimensional space, characterized by coordinate variables, position, orientation and scale. The uthors infer, from the V1 connectivity a "natural" geometry from which they can apply methods from differ-ential geometry. This type of mathematical construction could be interesting to investigate in the case of theretina with families of orientation selective cells, although the retinal connectivity between these cells is notthe same as the V1 orientation preference map structure [10, 9, 75, 83]. Biologically inspired vision systems
When the retina receives a visual stimulus, it determines whichcomponent of it is significant and needs to be further transmitted to the brain. This efficient coding heuris-tic has inspired many recent studies in developing biologically inspired systems, both for static image andmotion representations [52, 80, 82]. Two major applications of biologically inspired vision systems are retinalprostheses [46, 56, 72] and navigational robotics [21, 45]. Focusing on the second field of application, theability of a mobile device to navigate in its environment is of utmost interest, especially in order to avoiddangerous situations such as collisions. To be able to move, the robot requires a mapped representation ofits environment, but also the ability to interpret and process this representation. Motion processing mecha-nisms such as anticipation can be thus implemented to assess the efficiency of bio-inspired vision in obstacleavoidance.
Not applicable.
The Authors transfer to Springer the non-exclusive publication rights and they warrant that their contri-bution is original and that they have full power to make this grant.
The model source code is available upon request.
The Authors state that they do not have any conflicts of interest to disclose.
This work was supported by the National Research Agency (ANR), in the project "Trajectory",https://anr.fr/Project-ANR-15-CE37-0011, funding Selma Souihel’s PhD, and by the interdisciplinary In-stitute for Modelling in Neuroscience and Cognition (NeuroMod http://univ-cotedazur.fr/en/idex/projet-structurant/neuromod ) of the Université Côte d’Azur.
Both authors contributed to the final version of the manuscript. Bruno Cessac supervised the project. .7 Acknowledgment We warmly acknowledge Olivier Marre and Frédéric Chavane for their insightful comments, as well asMichael Berry, Matthias Hennig, Benoit Miramond, Stephanie Palmer and Laurent Perrinet for their thoroughfeedback as Jury members of Selma Souihel’s PhD.
References [1] S. Baccus and M. Meister. Fast and slow contrast adaptation in retinal circuitry. Neuron, 36(5):909–919,2002.[2] T. Baden, P. Berens, M. Bethge, and T. Euler. Spikes in mammalian bipolar cells support temporal layer-ing of the inner retina. Current Biology, 23(1):48 – 52, 2013.[3] T. Baden, P. Berens, K. Franke, M. R. Rosón, M. Bethge, and T. Euler. The functional diversity of retinalganglion cells in the mouse. Nature, 2016.[4] M. V. C. Baldo and S. A. Klein. Extrapolation or attention shift? Nature, 378:565–566, 1995.[5] H. Barlow. Possible principles underlying the transformation of sensory messages. Sensorycommunication, pages 217–234, 1961.[6] E. Baspinar, G. Citti, and A. Sarti. A geometric model of multi-scale orientation preference maps viagabor functions. Journal of Mathematical Imaging and Vision, 60:900–912, 2018.[7] G. Benvenuti, S. Chemla, A. Boonman, L. Perrinet, G. S. Masson, and F. Chavane. Anticipatory responsesalong motion trajectories in awake monkey area v1. bioRxiv, 2020.[8] M. Berry, I. Brivanlou, T. Jordan, and M. Meister. Anticipation of moving stimuli by the retina. Nature,398(6725):334—338, 1999.[9] W. Bosking, J. Crowley, D. Fitzpatrick, et al. Spatial coding of position and orientation in primary visualcortex. Nature neuroscience, 5(9):874–882, 2002.[10] W. Bosking, Y. Zhang, B. Schofield, and D. Fitzpatrick. Orientation selectivity and the arrangement ofhorizontal connections in tree shrew striate cortex. The Journal of Neuroscience, 17(6):2112–2127, 1997.[11] B. Cessac. A discrete time neural network model with spiking neurons. rigorous results on the sponta-neous dynamics. J. Math. Biol., 56(3):311–345, 2008.[12] B. Cessac. A discrete time neural network model with spiking neurons ii. dynamics with noise. Journalof Mathematical Biology, 62(6):863–900, 2011.[13] B. Cessac. Linear response in neuronal networks: From neurons dynamics to collective response. Chaos:An Interdisciplinary Journal of Nonlinear Science, 29(10):103105, 2019.[14] B. Cessac, I. Ampuero, and R. Cofre. Linear response for spiking neuronal networks with unboundedmemory, 2020. Submitted to J. Math. Neuro.[15] B. Cessac and T. Viéville. On dynamics of integrate-and-fire neural networks with adaptive conduc-tances. Frontiers in neuroscience, 2(2), July 2008.[16] S. Chemla, A. Reynaud, M. di Volo, Y. Zerlaut, L. Perrinet, A. Destexhe, and F. Chavane. Suppressivetraveling waves shape representations of illusory motion in primary visual cortex of awake primate.Journal of Neuroscience, 39(22):4282–4298, 2019.[17] E. Y. Chen, O. Marre, C. Fisher, G. Schwartz, J. Levy, R. A. da Silviera, and M. J. I. Berry. Alert responseto motion onset in the retina. Journal of Neuroscience, 33(1):120–132, 2013.[18] H. Choi, L. Zhang, M. S. Cembrowski, C. F. Sabottke, A. L. Markowitz, D. A. Butts, W. L. Kath, J. H.Singer, and H. Riecke. Intrinsic bursting of aii amacrine cells underlies oscillations in the rd1 mouseretina. Journal of Neurophysiology, 2014.
19] M. Colonnier and J. O’Kusky. Number of neurons and synapses in the visual cortex of different species.Rev Can Biol, 1981.[20] S. Coombes, Y. M. Lai, M. ¸Sayli, and R. Thul. Networks of piecewise linear neural mass models.European Journal of Applied Mathematics, 29(5):869–890, 2018.[21] T. Delbruck, K. Martin, S.-C. Liu, B. Linares-Barranco, and D. P. Moeys. Analog And DigitalImplementations Of Retinal Processing For Robot Navigation Systems. PhD thesis, ETH Zurich, 2016.[22] S. Deny, U. Ferrari, E. Macé, P. Yger, R. Caplette, S. Picaud, G. Tkaˇcik, and O. Marre. Multiplexedcomputations in retinal ganglion cells of a single type. Nature Communications, 2017.[23] R. Deriche. Using Canny’s criteria to derive a recursively implemented optimal edge detector.International Journal of Conputer Vision, 1(2):167–187, May 1987.[24] A. Destexhe, Z. Mainen, and T. Sejnowski. An efficient method for computing synaptic conductancesbased on a kinetic model of receptor binding. Neural Computation, 6(1):14—18, 1994.[25] A. Destexhe, Z. Mainen, and T. Sejnowski. Synthesis of models for excitable membranes, synap-tic transmission and neuromodulation using a common kinetic formalism. Journal of ComputationalNeuroscience, 1(3):195—230, 1994.[26] J. Dowling. Retina: An overview. In Reference Module in Biomedical Sciences. Elsevier, 2015.[27] A. Edelman. The probability that a random real gaussian matrix has k real eigenvalues, related distri-butions, and the circular law. Journal of Multivariate Analysis, 60(2):203 – 232, 1997.[28] G. A. Enciso, M. Rempe, A. V. Dmitriev, K. E. Gavrikov, D. Terman, and S. C. Mangel. A model ofdirection selectivity in the starburst amacrine cell network. Journal of Computational Neuroscience,18(3):567–578, June 2010.[29] T. Euler, P. Detwiler, and W. Denk. Directionally selective calcium signals in dendrites of starburstamacrine cells. Nature, 418:845–852, 2002.[30] O. Faugeras, J. Touboul, and B. Cessac. A constructive mean field analysis of multi populationneural networks with random synaptic weights and stochastic inputs. Frontiers in ComputationalNeuroscience, 3(1), 2009.[31] W. Freeman and E. Adelson. The design and use of steerable filters. IEEE Transactions on PatternAnalysis and Machine Intelligence, 13(9):891–906, 1991.[32] F. R. Gantmacher. the theory of matrices. AMS Chelsea Publishing, 1998.[33] J.-M. Geusebroek, A. Smeulders, and J. Weijer. Fast anisotropic gauss filtering. IEEE Transactions onImage Processing, 2350, 2003.[34] V. Girko. Circular law. Theor. Prob. Appl, 29:694–706, 1984.[35] T. Gollisch and M. Meister. Eye smarter than scientists believed: neural computations in circuits of theretina. Neuron, 65(2):150–164, Jan. 2010.[36] T. Hosoya, S. A. Baccus, and M. Meister. Dynamic predictive coding by the retina. Nature, 436:71–77,2005.[37] D. H. Hubel and T. N. Wiesel. Receptive fields of optic nerve fibres in the spider monkey. J. Physiol.,154:572–80, Dec. 1960.[38] J. Jacoby, Y. Zhu, S. H. DeVries, and G. W. Schwartz. An amacrine cell circuit for signaling steadyillumination in the retina. Cell reports, 13(12):2663–2670, 2015.[39] D. Jancke, W. Erlaghen, G. Schöner, and H. Dinse. Shorter latencies for motion trajectories than forflashes in population responses of primary visual cortex. Journal of Physiology, 556:971–982, 2004.[40] J. Johnston and L. Lagnado. General features of the retinal connectome determine the computation ofmotion anticipation. Elife, 2015.
41] J. Johnston, S.-H. Seibel, L. S. A. Darnet, S. Renninger, M. Orger, and L. Lagnado. A retinal circuitgenerating a dynamic predictive code for oriented features. Neuron, 102(6):1211 – 1222.e3, 2019.[42] M. Kähne, S. Rüdiger, A. H. Kihara, and B. Lindner. Gap junctions set the speed and nucleation rate ofstage i retinal waves. PLOS Computational Biology, 15(4):1–15, 2019.[43] D. Karvouniari, L. Gil, O. Marre, S. Picaud, and B. Cessac. A biophysical model explains the oscillatorybehaviour of immature starburst amacrine cells. Scientific Reports, 9:1859, 2019.[44] D. B. Kastner and S. A. Baccus. Spatial segregation of adaptation and predictive sensitization in retinalganglion cells. Neuron, 79(3):541 – 554, 2013.[45] H. Lehnert, M. Escobar, and M. Araya. Retina-inspired visual module for robot navigation in complexenvironments. In 2019 International Joint Conference on Neural Networks (IJCNN), pages 1–8, 2019.[46] C. Morillas, S. Romero, A. Martínez, F. Pelayo, L. Reyneri, M. Bongard, and E. Fernández. A neuro-engineering suite of computational tools for visual prostheses. Neurocomputing, 70(16):2817 – 2827,2007. Neural Network Applications in Electrical Engineering Selected papers from the 3rd InternationalWork-Conference on Artificial Neural Networks (IWANN 2005).[47] R. Nelson and H. Kolb. On and off pathways in the vertebrate retina and visual system. The VisualNeurosciences, 1:260–278, 2004.[48] H. H. Nguyen and S. O’Rourke. The Elliptic Law. International Mathematics Research Notices,2015(17):7620–7689, 2014.[49] R. Nijhawan. Motion extrapolation in catching. Nature, 370:256–257, 1994.[50] R. Nijhawan. Visual decomposition of colour through motion extrapolation. Nature, 386:66–69, 1997.[51] W.-Q. Niu and J.-Q. Yuan. Recurrent network simulations of two types of non-concentric retinal ganglioncells. Neurocomputing, 70(13):2576 – 2580, 2007. Selected papers from the 3rd International Conferenceon Development and Learning (ICDL 2004) Time series prediction competition: the CATS benchmark.[52] R. F. Oliveira and A. C. Roque. A biologically plausible neural network model of the primate primaryvisual system. Neurocomputing, 44-46:957 – 963, 2002. Computational Neuroscience Trends in Research2002.[53] B. Olshausen and D. Field. Sparse coding with an overcomplete basis set: A strategy employed by V1?Vision Research, 37:3311–3325, 1998.[54] B. Ölveczky, S. Baccus, and M. Meister. Segregation of object and background motion in the retina.Nature, 423:401–408, 2003.[55] S. E. Palmer, O. Marre, M. J. Berry, and W. Bialek. Predictive information in a sensory population.Proceedings of the National Academy of Sciences, 112(22):6908–6913, 2015.[56] N. Parikh, L. Itti, and J. Weiland. Saliency-based image processing for retinal prostheses. Journal ofNeural Engineering, 7(1):016006–10, Jan. 2010.[57] L. A. Remington. Chapter 4 - retina. In L. A. Remington, editor, Clinical Anatomy and Physiology of theVisual System (Third Edition), pages 61 – 92. Butterworth-Heinemann, Saint Louis, third edition edition,2012.[58] D. Ruelle and F. Takens. On the nature of turbulence. Comm. Math. Phys., 20:167–192, 1971.[59] J. Salisbury and S. Palmer. Optimal prediction in the retina and natural motion statistics. Journal ofStatistical Physics, 162, 2016.[60] A. J. Sederberg, J. N. MacLean, and S. E. Palmer. Learning to make external sensory stimulus predictionsusing internal correlations in populations of neurons. Proceedings of the National Academy of Sciences,115(5):1105–1110, 2018.
61] R. Segev, J. Puchalla, and M. J. B. II. Functional organization of ganglion cells in the salamander retina.J Neurophysiol, 95:2277–2292, 2006.[62] E. Seneta. Non-negative Matrices and Markov Chains. Springer, 2006.[63] E. Sernagor and M. Hennig. Chapter 49 - retinal waves: Underlying cellular mechanisms and theoreticalconsiderations. In J. L. Rubenstein and P. Rakic, editors, Cellular Migration and Formation of NeuronalConnections, pages 909 – 920. Academic Press, Oxford, 2013.[64] S. Sethuramanujam, G. B. Awatramani, and M. M. Slaughter. Cholinergic excitation complements glu-tamate in coding visual information in retinal ganglion cells. The journal of physiology, pages 464–475,May 2018.[65] S. Sethuramanujam, A. J. McLaughlin, G. deRosenroll, A. Hoggarth, D. J. Schwab, and G. B. Awatramani.A central role for mixed acetylcholine/gaba transmission in direction coding in the retina. Neuron, 2016.[66] J. Snellman, T. Kaur, Y. Shen, and S. Nawy. Regulation of on bipolar cell activity. Progress in retinal andeye research, 27(4):450–63, 2008.[67] S. Souihel. Generic and specific computational principles for visual anticipation of motion trajectories.Phd thesis, Université Nice Côte d’Azur ; EDSTIC, Dec. 2019.[68] M. Srinivas, R. Rozental, T. Kojima, R. Dermietzel, M. Mehler, D. F. Condorelli, J. A. Kessler, and D. C.Spray. Functional properties of channels formed by the neuronal gap junction protein connexin36. J.Neurosci., 19(22):9848–9855, 1999.[69] M. Srinivasan, S. Laughlin, and A. Dubs. Predictive coding: A fresh view of inhibition in the retina.Proceedings of the Royal Society of London. Series B, Biological Sciences, 216(1205):427–459, 1982.[70] M. Subramaniyan, A. S. Ecker, S. S. Patel, R. J. Cotton, M. Bethge, X. Pitkow, P. Berens, and A. S. To-lias. Faster processing of moving compared with flashed bars in awake macaque v1 provides a neuralcorrelate of the flash lag illusion. Journal of Neurophysiology, 2018.[71] M. Tauchi and R. Masland. The shape and arrangement of the cholinergic neurons in the rabbit retina.Proceedings of the Royal Society of London. Series B, Containing papers of a Biological character. RoyalSociety (Great Britain), 223:101–19, 1984.[72] T. K. Tran. Large scale retinal modeling for the design of new generation retinal prostheses, 2015.[73] S. Trenholm, D. Schwab, V. Balasubramanian, and G. Awatramani. Lag normalization in an electricallycoupled neural network. Nature medicine, 16, 2013.[74] J. J. Tukker, W. R. Taylor, and R. G. Smith. Direction selectivity in a model of the starburst amacrine cell.Visual Neuroscience, 2004.[75] T. Tversky and R. Miikkulainen. Modeling directional selectivity using self-organizing delay-adaptationmaps. Neurocomputing, 44-46:679 – 684, 2002. Computational Neuroscience Trends in Research 2002.[76] M. Unser. Fast gabor-like windowed fourier and continuous wavelet transforms. IEEE Signal ProcessingLetters, 1(5):76–79, 1994.[77] R. L. D. Valois and K. K. D. Valois. Vernier acuity with stationary moving gabors. Vision Research,31(9):1619 – 1626, 1991.[78] D. I. Vaney, B. Sivyer, and W. R. Taylor. Direction selectivity in the retina: symmetry and asymmetry instructure and function. Nature Reviews Neuroscience, 2012.[79] B. Völgyi, F. Pan, D. L. Paul, J. T. Wang, A. D. Huberman, and S. a. Bloomfield. Gap Junctions AreEssential for Generating the Correlated Spike Activity of Neighboring Retinal Ganglion Cells. PLoSOne, 8(7), 2013.[80] H. Wei and Q. Zuo. A biologically inspired neurocomputing circuit for image representation.Neurocomputing, 164:96 – 111, 2015.
81] W. Wei, A. Hamby, K. Zhou, and M. Feller. Development of asymmetric inhibition underlying directionselectivity in the retina. Nature, 469(7330):402–406, 2010.[82] J. Xu, S. H. Park, and X. Zhang. A bio-inspired motion sensitive model and its application to estimatinghuman gaze positions under classified driving conditions. Neurocomputing, 345:23 – 35, 2019. DeepLearning for Intelligent Sensing, Decision-Making and Control.[83] X. Xu, W. Bosking, G. Sáry, J. Stefansic, D. Shima, and V. Casagrande. Functional organization of visualcortex in the owl monkey. The Journal of neuroscience, 24(28):6237, 2004.[84] Y. Yu and T. Sing Lee. Adaptive contrast gain control and information maximization. Neurocomputing,65-66:111 – 116, 2005. Computational Neuroscience: Trends in Research 2005.[85] Y. Zerlaut, S. Chemla, F. Chavane, and A. Destexhe. Modeling mesoscopic cortical dynamics using amean-field model of conductance-based networks of adaptive exponential integrate-and-fire neurons.Journal of Computational Neuroscience, 2018.[86] J. Zheng, S. Lee, and Z. J. Zhou. A transient network of intrinsically bursting starburst cells underliesthe generation of retinal waves. Nat Neurosci, 9(3):363–371, 2006. ppendices A Parameters of the model
Function Parameter Value Unit K B,S (eq. (56)) σ (center) 90 µmσ (surround) 290 µmA (center) 1.2 mVA (surround) 0.2 mV K T ( t ) (eq. (57)) µ msµ msσ msσ msK K BCells dynamics τ a msh B . e − mV − .ms − θ B mVτ B ms ACells dynamics τ A msw [0, 1] ms − ACells connectivity ξ { , , , } mm ¯ n { , , , } unitless σ n GCells dynamics a p . unitless σ p µmτ G msh G . e − unitless α G Hz/mVθ G mVN maxG Hz B Spatio-temporal filtering
B.1 Receptive Fields
The spatial kernel of the Bcell i is modelled with a difference of Gaussians (DOG): K B i ,S ( x, y ) = A π √ det C e − ˜ X i .C − .X i − A π √ det C e − ˜ X i .C − .X i , (56) here X i = (cid:18) x − x i y − y i (cid:19) , (cid:101) denotes the transpose, x i and y i are the coordinates of the receptive field centerwhich coincide with the coordinates of the cell, C , C are positive definite matrix whose main principal axisrepresent the preferred orientation. For circular DOGs(no preferred orientation) C ≡ σ I , C ≡ σ I where I is the identity matrix in 2 dimensions. The two Gaussians of the DOG are thus concentric. They have thesame principal axes. X i has the physical dimension of a length ( mm ) thus the entries of C a , a = 1 , areexpressed in mm . A a , a = 1 . . . have the dimension of mV so that the convolution (1) has the dimensionof a voltage.We model the temporal part of the RF with a difference of non concentric Gaussians whose integral onthe time domain is zero. This kernel well fits the shape of the temporal projection of the bipolar RF observedin experiments [67]. K T ( t ) = (cid:32) K √ πσ e − ( t − µ σ − K √ πσ e − ( t − µ σ (cid:33) H ( t ) (57)where H ( t ) is the Heaviside function. The parameters µ b , σ b , b = 1 , have the dimension of a time ( s )whereas K b are dimensionless. The following condition must hold to ensure the continuity of K T ( t ) at zero: K σ e − µ σ = K σ e − µ σ . (58)Thus, K B i ( x, y,
0) = 0 . In addition, we require that the integral of a constant stimulus converges to zero, sothat the cell is only reactive to changes. This reads: K Π (cid:18) µ σ (cid:19) = K Π (cid:18) µ σ (cid:19) , (59)where: Π( x ) = 1 √ π (cid:90) x −∞ e − y dy, (60)is the cumulative distribution function of the standard Gaussian probability. B.2 Numerical convolution
Here we describe the method use to numerically integrate the convolution (1) of the receptive field K B i ,S (56) in the spatial domain with a stimulus S . For the sake of clarity we restrict the computation to oneGaussian in the DOG. The extension to a difference of Gaussian is straightforward. In the following, weconsider a spatially discretized stimulus. When dealing with a 2D stimulus, we have to integrate over twoaxis. In the case where the eigenvectors of the 2D of Gaussians are the axis of integration, the spatial filter isseparable in the stimulus coordinate system. Considering the stimulus as a grid of pixels, we can integrateusing the following discretization : let L x be the size of the stimulus along the x axis in pixels, L y its sizealong the y axis, and δ the pixel length. We set S ij ( t ) ≡ S ( iδ, jδ, u ) , with i = 0 , . . . , L x δ and j = 0 , . . . , L y δ .The spatial convolution becomes then : (cid:104) K B i x,y ∗ S (cid:105) ( t ) = 12 πσ x σ y (cid:90) (cid:90) R S ( x, y, t ) e − ( x − x σ x − ( y − y σ y dxdy = (cid:88) i,j S ij ( t ) (cid:20) erf ( i + δ − x √ σ x ) − erf ( i − x √ σ x ) (cid:21) (cid:20) erf ( j + δ − y √ σ y ) − erf ( j − y √ σ y ) (cid:21) In the case where the eigenvectors of the 2D of Gaussians are not the axes of integration, the spatial filteris not separable in the stimulus coordinates system. There exists methods that perform the computation y making a linear combination of basis filters [31], others that use Fourier based deconvolution techniques[76], and others using recursive filtering techniques [23]. However, these methods are of high computationalcomplexity. We choose instead to use a computer vision method from Geusenroek et al. [33].It is based on a projection in a non-orthogonal basis, where the first axis is x and the second is parametrizedby a angle φ (see Fig. 17). The new standard deviations read : σ x (cid:48) = σ x σ y (cid:112) σ x cos θ + σ y sin θ σ φ = (cid:112) σ y cos θ + σ x sin θ sin φ with tan( φ ) = σ y cos θ + σ x sin θ ( σ x − σ y ) cos θ sin θ and σ x (cid:54) = σ y (in the orientation sensitive case). Figure 17:
Filter transformation description [33]. The original system of axes is represented byx and y, and the ellipse system of axes by u and v. φ represents the angle of the second axis ofthe non-orthogonal basis. The integration domain of a pixel is limited by four lines of equations : x = iδ , x = ( i + 1) δ , y = jδ and y = ( j + 1) δ . Rewriting these fours equations in the new system ofaxes through a coordinate change enables us to write the equation (61) . We adapt the implementation to the spatially discretized stimulus, using an integration scheme similarto the one introduced in the separable case. The spatial convolution -reads now: (cid:104) K B i x,y ∗ S (cid:105) ( x , y , t ) = σ x (cid:48) (cid:112) π (cid:80) ( i ; j ) ∈ [0 ,s x ] × [0 ,s y ] (cid:82) ( y +1) δ sin( φ ) y δ sin( φ ) C ij e ( y (cid:48)− y σ φ [ erf ( ( − cos( φ ) y (cid:48) + x +1) δ − x √ σ x (cid:48) ) − erf ( ( − cos( φ ) y (cid:48) + x ) δ − x √ σ x (cid:48) )] dy (cid:48) (61)The integral is then computed numerically. The advantage of this formulation is to replace a two dimensionintegration by a one dimensional. Random connectivity
Here we define the random connectivity matrix from ACell to BCells considered in section 2.3.3. Eachcell (ACell and BCell) has a random number of branches (dendritic tree), each of which has a random lengthand a random angle with respect to the horizontal axis. The length of branches L follow an exponentialdistribution: f L ( l ) = 1 ξ e − lξ , l ≥ . (62)with spatial scale ξ . The number of branches n is also a random variable, Gaussian with mean ¯ n and variance σ n . The angle distribution is taken to be isotropic in the plane, i.e. uniform on [0 , π [ . When a branch of anACell A intersects a branch of a BCell B there is a chemical synapse from A to B.Here, we assume that both cells types have the same probability distributions for branches, thus neglect-ing the actual shape of ACell and BCell dendritic trees. On biological grounds, this assumption is relevant ifwe consider the shape of BCell dendritic tree in the Inner Plexiform Layer (IPL) (see e.g https://webvision.med.utah.edu/book/part-iii-retinal-circuits/roles-of-ACell-cells/ Fig. 5, 16, 17). Whileout of the IPL BCells have the form of a dipole, in the IPL their dendrites have a form well approximated byour -dimensional model. A potential refinement would consist of considering different set of parameters inthe probability laws respectively defining BCell and ACell dendritic tree.We show, in Fig. 18 an example of connectivity matrix produced this way, as well as the probability thattwo branches intersect as a function of the distance of the two cells.We now compute this probability. We use the standard notation in probability theory where the randomvariable is written in capitals and its realization in small letter. Thus, F X ( x ) = P [ X < x ] is the cumulativedistribution function of the random variable X and f X ( x ) = dF X dx its density.We consider the oriented connection between the cell A (ACell), of coordinates ( x A , y A ) to a cell B (BCell),of coordinates ( x B , y B ) , so that the distance between the two cells is d AB = (cid:113) ( x B − x A ) + ( y B − y A ) .The vector (cid:126)AB makes an oriented angle η = (cid:92) (cid:16) (cid:126)AB, (cid:126)Ax (cid:17) with the positive horizontal axis, where: η = (cid:40) arctan y B − y A x B − x A , if x B > x A π + arctan y B − y A x B − x A , if x B < x A . (63)Here, we neglect the effects of boundaries - taking, e.g., an infinite lattice, or periodic boundary conditions- so that the probability to connect A to B is invariant by rotation. Thus, we compute this probability in thefirst quadrant x B > x A , y B > y A . In this case η = arctan y B − y A x B − x A .Each cell has a random number of branches (dendritic tree), each of which has a random length and arandom angle with respect to the horizontal axis (Fig. 19). The length of branches L follow the exponentialdistribution (62): f L ( l ) = ξ e − lξ , l ≥ , with repartition function: F L ( l ) = 1 − e − lξ . (64)The spatial scale ξ favours short range connections. The number of branches N distribution follows an nor-mal distribution with mean ¯ n and variance σ n . The angle distribution is taken to be isotropic in the plane, i.e.uniform on [0 , π [ .We compute the probability that a branch of ACell A , of length L A , intersects, at point C , a branch ofBCell B , of length L B . We note α , the oriented angle (cid:92) (cid:16) (cid:126)Ax, (cid:126)AC (cid:17) ; β , the oriented angle (cid:92) (cid:16) (cid:126)Bx, (cid:126)BC (cid:17) ; θ , theoriented angle (cid:92) (cid:16) (cid:126)AB, (cid:126)AC (cid:17) . In the first quadrant, α = θ + η . Note that the condition to be in the first quadrantconstraints η but not α . B i A j P ( d ) d ExpTh Figure 18:
Random connectivity. Left.
Example of a random connectivity matrix from ACells A j to BCells cells B i . White points correspond to connection from A j to B i . Right.
Probability P ( d ) that two branches intersect as a function of the distance between two cells. ’Exp’ corresponds tonumerical estimation and ’Th’ corresponds to the theoretical prediction. Here, ξ = 2 .Figure 19: Geometry of connection between neurons. α ( β ) is the angle of the neuron A ’s branchwith length L A (neuron B ’s branch with length L B ) with respect to the horizontal axis. θ is theangle between the segment connecting AB and the branch A . C represents the virtual point thatlies at the intersection of the branches of length L A and L B . d AB (resp. d AC , d BC ) denotes thedistance between A and B (resp. A, C and B,C). Note that d AC ≤ L A , d BC ≤ L B .51 rom the sin rule we have: sin( β − α + θ ) d AC = sin θd BC = sin( β − α ) d AB . This holds however if A, B, C is a triangle,that is, if the two branches are long enough to intersect at C which reads: ≤ d AC = sin( β − α + θ )sin( β − α ) d AB ≤ L A and ≤ d BC = sin θ sin( β − α ) d AB ≤ L B . These are necessary and sufficient conditions for the branches tointersect.Note that the positivity of these quantities imposes conditions linking the angles α, β, η with θ = α − η .1. If sin ( β − α ) > ⇔ < β − α < π ⇔ α < β < π + α we must have sin θ > ⇔ < θ = α − η < π so that η < α < π + η , and, sin ( β − α + θ ) > ⇔ < β − α + θ = β − η < π so that η < β < π + η (because η ≥ ). All these constraints are satisfied if η < α < β < min ( π + α, π + η ) = π + η .2. If sin ( β − α ) < ⇔ − π < β − α < ⇔ − π + α < β < α we must have sin θ < ⇔ − π < θ = α − η < so that − π + η < α < η , and, sin ( β − α + θ ) < ⇔ − π < β − α + θ = β − η < so that − π + η < β < η .All these constraints are satisfied if max ( − π + α, − π + η ) = − π + η < β < α < η .Modulo these conditions, the conditional probability ρ c | ( α,β ) to have intersection given the angles α, β is: ρ c | ( α,β ) = P [ Connection | α, β ]= P (cid:20) ≤ sin ( β − α + θ )sin ( β − α ) d AB ≤ L A , ≤ sin θ sin ( β − α ) d AB ≤ L B | α, β (cid:21) . Using the cumulative distribution function (64) of the exponential distribution, and the independence of L A , L B this gives: ρ c | ( α,β ) = (cid:18) − F L (cid:18) sin ( β − α + θ )sin ( β − α ) d AB (cid:19) (cid:19) (cid:18) − F L (cid:18) sin θ sin ( β − α ) d AB (cid:19) (cid:19) = e − dABξ sin ( α + β − η ) sin ( β − α ) The probability to connect the two branches in the first quadrant is then: ρ c ( d AB , η ) = 14 π (cid:90) π + ηα = η (cid:90) π + ηβ = α e − dABξ sin ( α + β − η ) sin ( β − α ) dα dβ + (cid:90) ηα = − π + η (cid:90) αβ = − π + η e − dABξ sin ( α + β − η ) sin ( β − α ) dα dβ, (65)which depends on the distance between the two cells and their angle η , depending parametrically on thecharacteristic length ξ . Note that the condition of positivity of the sine ratio ensures an exponential decay ofthe probability as d AB increases. Remark.
The positivity of arguments in the exponential implies that: π ρ c ( d AB , η ) ≤ (cid:90) π + ηα = η (cid:90) π + ηβ = α dα dβ + (cid:90) ηα = − π + η (cid:90) αβ = − π + η dα dβ = π , so that ρ c ( d AB , η ) ≤ . D Linear analysis
D.1 General solution of the linear dynamical system
Here we consider the dynamical system (34), d (cid:126) X dt = L . (cid:126) X + (cid:126) F ( t ) , whose general solution is: (cid:126) X ( t ) = (cid:90) tt e L ( t − s ) . (cid:126) F ( s ) ds, he behaviour of this integral depends on the spectrum of L . The difficulty is that L is not diagonalisable(because of the activity term h B I N,N ). We write it in the form: L = M (cid:122) (cid:125)(cid:124) (cid:123)(cid:32) − I N,N τ B W AB W BA − I N,N τ A (cid:33) N,N N,N N,N N,N − I N,N τ a (cid:124) (cid:123)(cid:122) (cid:125) D + N,N N,N N,N N,N N,N N,N h B I N,N N,N N,N (cid:124) (cid:123)(cid:122) (cid:125) J We assume that the matrix M is diagonalisable. Even in this case, L is not diagonalisable because of theJordan matrix J . We note (cid:126)φ β , the normalized eigenvectors of M and λ β the corresponding eigenvalues, with β = 1 . . . N . The eigenvalues of D are then the N eigenvalues of M plus N eigenvalues − τ a . We notethem λ β too, with λ β = − τ a , β = 2 N + 1 . . . N . The eigenvectors of D have the form: (cid:126) P β = (cid:18) (cid:126)φ β (cid:126) N (cid:19) , β = 1 . . . N ; (cid:126)e β , β = 2 N + 1 . . . N, where (cid:126) N is the N dimensional vector with entries and (cid:126)e β is the canonical basis vector in direction β . Thematrix P made by the columns (cid:126) P β is the matrix which diagonalizes D . We note Λ = P − DP the diagonalform where Λ =
Diag { λ β , β = 1 . . . N } . P and P − have the block form: P = (cid:18) Φ 0 N,N N, N I N,N (cid:19) ; P − = (cid:18) Φ − N,N N, N I N,N (cid:19) , (66)where Φ is the matrix whose columns are the eigenvectors (cid:126)φ β of M . This form implies that P αβ = P − αβ = δ αβ ,for β = 2 N + 1 . . . N .We now compute e L t using the series expansion e L t = (cid:80) + ∞ n =0 t n n ! ( D + J ) n . Using the relations: J = 0 N, N ; D n . J = (cid:18) − τ a (cid:19) n J , one proves that: ( D + J ) n = D n + J . n − (cid:88) k =0 (cid:18) − τ a (cid:19) k D n − − k . Therefore: e L t = e D t + J . + ∞ (cid:88) n =1 t n n ! n − (cid:88) k =0 (cid:18) − τ a (cid:19) k D n − − k . We use the matrices P , P − to write it in the form: e L t = P .e Λ t . P − + J . + ∞ (cid:88) n =1 t n n ! n − (cid:88) k =0 (cid:18) − τ a (cid:19) k P . Λ n − − k . P − . From relation (36) we obtain, for the entries of (cid:126) X ( t ) : X α ( t ) = N (cid:88) β,γ =1 P αβ P − βγ (cid:90) tt e λ β ( t − s ) F γ ( s ) ds + N (cid:88) δ,β,γ =1 J α,δ P δβ P − βγ (cid:90) tt + ∞ (cid:88) n =1 ( t − s ) n n ! n − (cid:88) k =0 (cid:18) − τ a (cid:19) k λ n − − kβ F γ ( s ) ds. (67) e consider the first term of this equation. We use F γ = F B i , γ = i = 1 . . . N (BCells). We recall that,from (13), F B i ( t ) = V idrive τ B + dV idrive dt , so that: (cid:90) tt e λ β ( t − s ) F γ ( s ) ds = V γ drive ( t ) + (cid:18) τ B + λ β (cid:19) (cid:90) tt e λ β ( t − s ) V γ drive ( s ) ds. (68)Moreover, N (cid:88) β =1 3 N (cid:88) γ =1 P αβ P − βγ V γ drive ( t ) = N (cid:88) γ =1 V γ drive ( t ) N (cid:88) β =1 P αβ P − βγ = N (cid:88) γ =1 V γ drive ( t ) δ αγ = V α drive ( t ) , We extend the definition of the drive term (1) to N -dimensions such that V α drive ( t ) = 0 if α > N . Thus: N (cid:88) β,γ =1 P αβ P − βγ (cid:90) tt e λ β ( t − s ) F γ ( s ) ds = V α drive ( t ) + N (cid:88) β =1 (cid:18) τ B + λ β (cid:19) N (cid:88) γ =1 P αβ P − βγ (cid:90) t −∞ e λ β ( t − s ) V γ drive ( s ) ds. We decompose the sum over β in sums: β = 1 . . . N corresponding to BCells; β = N + 1 . . . N corresponding to ACells; β = 2 N + 1 . . . N corresponding to activities of BCells. We define (eq. (38) in thetext): E BB,α ( t ) = N (cid:88) β =1 (cid:18) τ B + λ β (cid:19) N (cid:88) γ =1 P αβ P − βγ (cid:90) tt e λ β ( t − s ) V γ drive ( s ) ds, α = 1 . . . N, corresponding to the indirect effect, via the ACells connectivity, of the drive on BCells voltages. The term E BA,α ( t ) = N (cid:88) β = N +1 (cid:18) τ B + λ β (cid:19) N (cid:88) γ =1 P αβ P − βγ (cid:90) tt e λ β ( t − s ) V γ drive ( s ) ds, α = N + 1 . . . N, (eq. (39) in the text) corresponds to the effect of BCell drive on ACell voltages. The third term : N (cid:88) β =2 N +1 (cid:18) τ B + λ β (cid:19) N (cid:88) γ =1 P αβ P − βγ (cid:90) tt e λ β ( t − s ) V γ drive ( s ) ds = 0 , α = 2 N + 1 . . . N. because P − βγ = δ βγ .To compute the second term in eq. (67), we first first remark that J α,δ = 0 , if α = 1 . . . N , and J α,δ = h B δ α − N,δ , if α = 2 N + 1 . . . N , so that this term is non zero only if α = 2 N + 1 . . . N (BCells activities).Also F γ (cid:54) = 0 for γ = 1 . . . N , while P − βγ = δ βγ for β = 2 N + 1 . . . N . Therefore, for α = 2 N + 1 . . . N thesecond term in X α ( t ) is: h B N (cid:88) β =1 N (cid:88) γ =1 P α − Nβ P − βγ (cid:90) tt + ∞ (cid:88) n =1 ( t − s ) n n ! n − (cid:88) k =0 (cid:18) − τ a (cid:19) k λ n − − kβ F γ ( s ) ds We now simplify the series: + ∞ (cid:88) n =1 ( t − s ) n n ! n − (cid:88) k =0 (cid:18) − τ a (cid:19) k λ n − − kβ = + ∞ (cid:88) n =1 ( t − s ) n n ! λ n − β n − (cid:88) k =0 (cid:18) − τ a λ β (cid:19) k = + ∞ (cid:88) n =1 ( t − s ) n n ! λ n − β − (cid:16) − τ a λ β (cid:17) n τ a λ β λ β
11 + τ a λ β (cid:34) + ∞ (cid:88) n =1 ( λ β ( t − s ) ) n n ! − + ∞ (cid:88) n =1 ( − t − sτ a ) n n ! (cid:35) = 1 λ β + τ a (cid:104) e λ β ( t − s ) − e − t − sτa (cid:105) The time integral is computed the same way as eq. (68): (cid:90) tt + ∞ (cid:88) n =1 ( t − s ) n n ! n − (cid:88) k =0 (cid:18) − τ a (cid:19) k λ n − − kβ F γ ( s ) ds = 1 λ β + τ a (cid:20) (cid:90) tt e λ β ( t − s ) F γ ( s ) ds − (cid:90) tt e − t − sτa F γ ( s ) ds (cid:21) = 1 λ β + τ a (cid:20) (cid:18) τ B + λ β (cid:19) (cid:90) tt e λ β ( t − s ) V γ drive ( s ) ds − (cid:18) τ B − τ a (cid:19) (cid:90) tt e − t − sτa V γ drive ( s ) ds (cid:21) . Similarly to eq. (38), (39) in the text we introduce: E Ba,α ( t ) = h B N (cid:88) β =1 N (cid:88) γ =1 P α − Nβ P − βγ λ β + τ a (cid:16) τ B + λ β (cid:17) (cid:82) tt e λ β ( t − s ) V γ drive ( s ) ds − (cid:16) τ B − τ a (cid:17) (cid:82) tt e − t − sτa V γ drive ( s ) ds. , α = 2 N +1 . . . N, corresponding to the action of BCells and ACells on the activity of BCells, via the network effect. Let usconsider in more detail the second term. From (66), (cid:80) Nβ =1 P α − Nβ P − βγ = δ α − Nγ thus: N (cid:88) β =1 N (cid:88) γ =1 P α − Nβ P − βγ (cid:90) tt e − t − sτa V γ drive ( s ) ds = N (cid:88) γ =1 (cid:90) tt e − t − sτa V γ drive ( s ) ds N (cid:88) β =1 P α − Nβ P − βγ (cid:124) (cid:123)(cid:122) (cid:125) δ α − Nγ = (cid:90) tt e − t − sτa V α − N drive ( s ) ds ≡ A α − N ( t ) , (eq. (41) in the text).This finally leads to ((40) in the text): E Ba,α ( t ) = h B N (cid:88) β =1 N (cid:88) γ =1 P α − Nβ P − βγ λ β + τ B λ β + τ a (cid:90) tt e λ β ( t − s ) V γ drive ( s ) ds + − τ B + τ a λ β + τ a A α − N ( t ) , α = 2 N +1 . . . N, D.2 Spectrum of L and stability of the dynamical system (34) Here, we assume that a BCell connects only one ACell, with a weight w + uniform for all BCells, so that W BA = w + I N,N , w + > . We also assume that ACell connect to BCell with a connectivity matrix W , notnecessarily, symmetric, with a uniform weight − w − , w − > , so that W AB = − w − W . We have shown inthe previous section that the N first eigenvalues and eigenvectors of L are given by the N eigenvalues andeigenvectors of M which reads now: M = (cid:32) − I N,N τ B − w − W w + I N,N − I N,N τ A (cid:33) . (69)We now show that this specific structure allows compute the spectrum of M in terms of the spectrum of W . .2.1 Eigenvalues and eigenvectors of M We note κ n , n = 1 . . . N , the eigenvalues of W ordered as | κ | ≤ | κ | ≤ · · · ≤ | κ n | and (cid:126)ψ n is thecorresponding eigenvector. We normalize (cid:126)ψ n so that (cid:126)ψ † n . (cid:126)ψ n = 1 where † is the adjoint. (Note that, as W isnot symmetric in general, eigenvectors are complex).We shall neglect the case where, simultaneously, τ = 0 and κ n = 0 for some n . Proposition.
For each n , there is a pair of eigenvalues λ ± n and eigenvectors (cid:126)φ ± n = c ± n (cid:32) (cid:126)ψ n ρ ± n (cid:126)ψ n (cid:33) of M with c ± n = (cid:113) ( ρ ± n ) (normalisation factor), and: ρ ± n = τ w − κ n (cid:0) ± √ − µ κ n (cid:1) , κ n (cid:54) = 0 , τ (cid:54) = 0; w + τ, κ n = 0 , τ (cid:54) = 0; ± (cid:113) − w + w − κ n , τ = 0 . (70)where: τ = (cid:18) τ A − τ B (cid:19) . and: τ AB = (cid:18) τ A + 1 τ B (cid:19) . Eigenvalues are given by: λ ± n = − τ AB ∓ τ √ − µ κ n , τ (cid:54) = 0; − τ A ∓ √− w − w + κ n , τ = 0 . (71)with: µ = w − w + τ ≥ , As a consequence, in addition to the N last eigenvalues − τ A , L admits N eigenvalues given by (71),while the N first columns of the matrix P (eigenvectors of L ) are: (cid:126) P β = 1 (cid:113) (cid:0) ρ − n (cid:1) (cid:126)ψ n ρ − n (cid:126)ψ n (cid:126) N ; (cid:126) P β + N = 1 (cid:113) (cid:0) ρ + n (cid:1) (cid:126)ψ n ρ + n (cid:126)ψ n (cid:126) N , β = n = 1 . . . N. (72)For the N last eigenvectors (cid:126) P β = (cid:126)e β , β = 2 N + 1 . . . N . Remark.
The structure of these eigenvectors is quite instructive. Indeed, the factors ρ ± n control the projec-tion of the eigenvectors (cid:126) P β on the space of ACells, thereby tuning the influence of ACells via lateral connectivity . Proof.
We use here the generic notation λ β , (cid:126)φ β , β = 1 . . . N for the eigenvalues and associated eigenvec-tors of M . If we assume that (cid:126)φ β is of the form (cid:126)φ β = (cid:32) (cid:126)ψ n ρ (cid:126)ψ n (cid:33) for some n , then, we have: M .(cid:126)φ β = (cid:32) − Iτ B − w − W w + I − Iτ A (cid:33) . (cid:32) (cid:126)ψ n ρ (cid:126)ψ n (cid:33) = (cid:16) − τ B − w − ρ κ n (cid:17) . (cid:126)ψ n (cid:16) − ρτ A + w + (cid:17) . (cid:126)ψ n = λ β (cid:32) (cid:126)ψ n ρ (cid:126)ψ n (cid:33) , hich gives: (cid:16) − τ B − w − ρ κ n (cid:17) = λ β (cid:16) − ρτ A + w + (cid:17) = λ β ρ, leading to: w − κ n ρ − τ ρ + w + = 0 , where: τ = 1 τ A − τ B . This gives, if κ n (cid:54) = 0 and τ (cid:54) = 0 : ρ ± n = 12 τ w − κ n (cid:16) ± (cid:112) − µ κ n (cid:17) , where: µ = w − w + τ ≥ . Thus, for each n , there are two eigenvalues: λ ± n = − τ AB ∓ τ (cid:112) − µ κ n , with: τ AB = 1 τ A + 1 τ B . Note that τ AB ≥ τ .If κ n = 0 , τ (cid:54) = 0 , ρ ± n = w + τ . Then: λ β = − τ B − w − ρ κ n Finally, if κ n (cid:54) = 0 , τ = 0 , ( τ A = τ B ), ρ ± n = − w + w − κ n and λ ± n = − τ B ± √− w − w + κ n . If κ n = 0 , τ = 0 there isno solution for ρ . End of proof.
Remark.
When µ = 0 , M is diagonal: the N first eigenvalues are − τ B , the N next eigenvalues are − τ A .We have, in this case: λ + n = − τ B and λ − n = − τ A . Therefore, in order to be coherent with this diagonalform of L when µ = 0 we order eigenvalues and eigenvectors of M such that the N first eigenvalues are λ β = λ + n , β = 1 . . . N , and the N next are λ β = λ − n , β = N + 1 . . . N. D.2.2 Stability of eigenmodesStability of eigenmodes when W is symmetric. If W is symmetric, its eigenvalues κ n are real, but λ β , β = 1 . . . N can be real or complex, depending on κ n , as µ is positive.We have four cases: • κ n < . Then, from (71), λ β s are real and there are two cases. If τ > the eigenvalues λ β , β = 1 . . . N can have a positive real part (unstable) while λ β , β = N + 1 . . . N has always a negative real part(stable); for τ < the situation is inverted. In both case, the eigenvalue λ β has a positive real part if: µ > − κ n τ A τ B ( τ B − τ A ) ≡ µ n,u , hich reads as well, using the definition of µ : w − w + > − τ A τ B κ n (73)Thus, τ A , τ B play a symmetric role. If τ = 0 ( τ A = τ B ), all eigenvalues are real. Eigenvalues λ − n are allstable. The eigenvalue λ + n becomes unstable if: w − w + > − τ A κ n , corresponding to (73). Proof
There are two cases. – τ > ⇔ τ A < τ B . λ β = − τ AB ± τ (cid:112) − µ κ n > ⇔ ± (cid:112) − µ κ n > ττ AB Only + is possible (because τ and τ AB are positive). This gives: − µ κ n > τ τ AB − τ τ AB = − τ A τ B ( τ B − τ A ) > µ κ n which is possible because κ n < . Thus, λ β is unstable if: µ > − κ n τ A τ B ( τ B − τ A ) ≡ µ n,u . – τ < ⇔ τ A > τ B . − τ AB ± τ (cid:112) − µ κ n > ⇔ ± (cid:112) − µ κ n < ττ AB Only − is possible (because τ < ). This gives: − µ κ n > τ τ AB , the same condition as in the previous item. – If τ = 0 , λ β = − τ A ∓ √− w − w + κ n so that eigenvalues are real. The eigenvalue with the minussign ( λ − n ) are all stable. The eigenvalue λ + n becomes unstable if: w − w + > − τ A κ n . End of proof. • κ n > . Then λ β , β = 1 . . . N are real or complex. If τ (cid:54) = 0 they are complex if: µ > κ n ≡ µ n,c . In this case the real part is − τ AB , the imaginary part is ± τ √ − µ κ n , and all eigenvalues arestable. If µ ≤ µ n,c eigenvalues λ β are real and all modes are stable as well. Indeed: If τ = 0 , all eigenvalues are equal to − τ AB , hence are stable. – If τ > ⇔ τ A < τ B . − τ AB ± τ (cid:112) − µ κ n > ⇔ ± (cid:112) − µ κ n > ττ AB which is not possible because ττ AB > whereas √ − µ κ n < . – If τ < ⇔ τ A > τ B . − τ AB ± τ (cid:112) − µ κ n > ⇔ ± (cid:112) − µ κ n < ττ AB Only − is possible because τ < . − µ κ n > τ τ AB which is not possible because ττ AB > whereas √ − µ κ n < . Stability of eigenmodes when W is asymmetric. If W is asymmetric, eigenvalues κ n are complex, κ n = κ n,r + i κ n,i . We write λ β = λ β,r + i λ β,i , β = 1 . . . N with: (cid:40) λ β,r = − τ AB ± τ √ √ a n + u n ; λ β,i = ± τ √ √ u n − a n , where a n = 1 − µ κ n,r and u n = (cid:113) ( 1 − µ κ n,r ) + 16 µ κ n,i = (cid:113) − µ κ n,r + 16 µ | κ n | . Note thatwe recover the real case when κ n,i = 0 by setting u n = a n .Instability occurs if a n + u n > τ τ AB , a condition depending on κ n,r and κ n,i ..