## 1

Evgenia Kartsaki  Université Côte d’Azur, Inria, Biovision team and Neuromod
Institute, Sophia Antipolis, France Biosciences Institute, Newcastle
University, Newcastle upon Tyne NE2 4HH, UK  Gerrit Hilgen  Biosciences
Institute, Newcastle University, Newcastle upon Tyne NE2 4HH, UK  Health and
Life Sciences, Applied Sciences, Northumbria University, Newcastle upon Tyne
NE1 8ST, UK  Evelyne Sernagor  Biosciences Institute, Newcastle University,
Newcastle upon Tyne NE2 4HH, UK  Bruno Cessac  Université Côte d’Azur, Inria,
Biovision team and Neuromod Institute, Sophia Antipolis, France

## How does the inner retinal network shape the ganglion cells receptive field
: a computational study

Evgenia Kartsaki  Université Côte d’Azur, Inria, Biovision team and Neuromod
Institute, Sophia Antipolis, France Biosciences Institute, Newcastle
University, Newcastle upon Tyne NE2 4HH, UK  Gerrit Hilgen  Biosciences
Institute, Newcastle University, Newcastle upon Tyne NE2 4HH, UK  Health and
Life Sciences, Applied Sciences, Northumbria University, Newcastle upon Tyne
NE1 8ST, UK  Evelyne Sernagor  Biosciences Institute, Newcastle University,
Newcastle upon Tyne NE2 4HH, UK  Bruno Cessac  Université Côte d’Azur, Inria,
Biovision team and Neuromod Institute, Sophia Antipolis, France

###### Abstract

We consider a model of basic inner retinal connectivity where bipolar and
amacrine cells interconnect, and both cell types project onto ganglion cells,
modulating their response output to the brain visual areas. We derive an
analytical formula for the spatio-temporal response of retinal ganglion cells
to stimuli taking into account the effects of amacrine cells inhibition. This
analysis reveals two important functional parameters of the network: (i) the
intensity of the interactions between bipolar and amacrine cells, and, (ii)
the characteristic time scale of these responses. Both parameters have a
profound combined impact on the spatiotemporal features of RGC responses to
light. The validity of the model is confirmed by faithfully reproducing
pharmacogenetic experimental results obtained by stimulating excitatory
DREADDs (Designer Receptors Exclusively Activated by Designer Drugs) expressed
on ganglion cells and amacrine cells subclasses, thereby modifying the inner
retinal network activity to visual stimuli in a complex, entangled manner. Our
mathematical model allows us to explore and decipher these complex effects in
a manner that would not be feasible experimentally and provides novel insights
in retinal dynamics.

##  2 Introduction

Visual processing involves some of the most complex neural networks in the
vertebrate central nervous system (Marr, , 1982; Chalupa and Werner, , 2004;
Besharse and Bok, , 2011; Daw, , 2012). The retina is the entry point to our
visual system. Located at the back of the eye, this thin neural tissue
receives the light that the cornea and lens have captured from different parts
of the visual scene, converts it into electrical signals and finally,
transmits these signals to the brain visual areas. In particular, light
follows a vertical excitatory pathway in the retina, from photoreceptors (PRs)
to bipolar cells (BCs) and onwards to retinal ganglion cells (RGCs), modulated
laterally by inhibitory interneurons; horizontal (HCs) and amacrine cells
(ACs) (Figure 1, left). RGCs serve as a bridge between the retina and the
brain, conveying highly processed and integrated signals from the upstream
retinal neurons to downstream visual processing cortical areas. Amazingly, the
human brain can recreate images from interpreting parallel streams of
information of about one million RGCs, the sole retinal output neurons. This
ability is partially due to the astonishing functional, anatomical and
molecular diversity across the retinal layers. However, how the different cell
classes interact to ultimately encode, via RGCs, a visual scene into spike
trains (action potentials) deciphered by the brain remains largely a mystery.

RGCs are indeed embedded in a complex network and their response is roughly
driven by two "controllers": 1) The output of BCs that includes both the
intrinsic response properties of these cells and the actions of HCs and ACs
upon them (lateral connectivity); 2) The direct input from ACs via chemical
synapses or gap junctions, helping spike synchrony between neighbor RGCs
(Masland, 2012b, ; Demb and Singer, , 2015). As a consequence, the response of
RGCs to visual stimuli does not only depend on the physiological
characteristics of the cells (Sanes and Masland, , 2015), but also on the
network they are embedded in, and, on the stimulus itself (Cessac, , 2022).
Previous studies have thus attempted to investigate how the inner retinal
neurons are organised into parallel circuits across different cell types and
converge onto RGCs (Wässle, , 2004; Gollisch and Meister, , 2010). This has
been studied extensively at the level of bipolar cells, leading to a fairly
good understanding of their function (Euler et al., , 2014). Other studies
have investigated the functional role of amacrine cell (ACs) types in retinal
processing (Asari and Meister, , 2012; Franke and Baden, , 2017; Diamond, ,
2017; Schröder et al., , 2020), suggesting either specific functions such as
direction selectivity (starbust ACs) or more general computations, like motion
anticipation (Berry et al., , 1999; Souihel and Cessac, , 2021). A number of
studies have also focused on how retinal interneurons contribute towards the
sensitivity to activity patterns of RGCs. For example, de Vries et al.,
(2011) studied the spatio-temporal effects of ACs on RGCs outputs and proposed
mechanisms related to predictive coding of the retina. Manu et al.,  (2022),
Ichinose et al.,  (2014), manipulated HCs and ACs, the first by injecting
patterns of current and the latter by using pharmacology, to directly measure
the effect on retinal receptive fields surround. Likewise, Protti et al.,
(2014) studied the spatial organisation of excitatory and inhibitory synaptic
inputs onto ON and OFF RGCs in the primate retina using pharmacology.
Nevertheless, the potential role of the BCs-ACs network on the RGCs response,
both from a theoretical and experimental perspective, hasn’t been sufficiently
explored yet. The scope of the present study is to make one step further in
this direction.

Numerous models of the retina have been proposed, with different levels of
biological detail and across multiple spatial and temporal scales.
Phenomenological models have become quite popular due to their simplicity and
efficiency to fit to experimental data and explain observations, such as the
linear-nonlinear model (Chichilnisky, , 2001; Paninski, , 2003; Simoncelli et
al., , 2004; Schwartz et al., , 2006), the generalized linear models (Pillow
et al., , 2008) or black-box convolutional neural networks of retinal function
(McIntosh et al., , 2017). However, these models are agnostic to the
biophysical details that underlie the input-output transform of the visual
system. On the other hand, mechanistic models exhibit a direct relationship
between the underlying biophysics of neurons and the parameters of the model.
Well known examples are the Hudgkin-Huxley (Hodgkin and Huxley, , 1952), the
Fitzhugh-Nagumo (FitzHugh, , 1969; Nagumo et al., , 1962), the Morris-Lecar
(Morris and Lecar, , 1981) models but they were mainly designed to describe
the spike generation while most of the neurons in the retina (except RGCs) are
not spiking. In addition, this type of models turns out to be cumbersome to
analyze and simulate, especially when dealing with collective network
dynamics.

The model developed in the present work falls into both categories and it is
definitely inspired by some previous work (Berry et al., , 1999; Chen et al.,
2013b, ; Souihel and Cessac, , 2021). It is a multi-stage, phenomenological
model, with simplifications regarding the complex retinal structure.
Nevertheless, it aims to be relatively precise from a biological perspective,
by realistically reproducing RGCs’ responses to light from the experimental
recordings. Our primary goal is to gain insights about the underlying
biophysical processes giving rise to certain experimentally observed
phenomena, rather than proposing an exact model of the retina.

On a theoretical ground, we consider a network of BCs connected via ACs, both
cell types being connected with chemical synapses to RGCs. Based on a
mathematical study, we derive an analytical formula of the RGCs receptive
field (RF) that takes into account the lateral ACs connectivity and shows how
the response of RGCs to spatio-temporal stimuli is shaped. Especially, we
emphasize the role of: (i), the average intensities of the interactions BCs -
ACs and ACs-BCs, and, (ii), the characteristic time scales of these cells
response. Varying these parameters acts on the shape of the response to light
with potential prominent effects such as a switch from monophasic to biphasic
in the temporal RF.

We illustrate our predictions by analysing experimental data obtained from the
pharmacological action of excitatory DREADDs (Designer Receptors Exclusively
Activated by Designer Drugs) in genetically modified mice. DREADDs are
activated by "designer drugs" such as clozapine-n-oxide (CNO), resulting in an
increase in free cytoplasmic calcium and profound increase in excitability in
the cells that express these DREADDs (Roth, , 2016). We found DREADD
expression both in subsets of RGCs and in ACs in the inner nuclear layer (INL)
(Hilgen et al., , 2022). In these conditions, CNO would act both on ACs and
RGCs providing a way to unravel the entangled effects of: (i) direct
excitation of DREADD-expressing RGCs and increase inhibitory input onto RGCs
originating from DREADD-expressing ACs; (ii) change in cells response time
scale (via a change in the membrane conductance), thereby providing an
experimental set up to validate our theoretical predictions.

In the following, we propose a simplified model for the BCs - ACs - RGCs
network. The concerted activity of BCs - ACs -RGCs in response to a visual
stimulus is described by a large dimensional dynamical system whose
mathematical study allows an explicit computation of the RF of BCs, ACs, RGCs
and, more generally, to anticipate the effects resulting from light
stimulation conjugated with the network activity due to ACs lateral
inhibition, in control conditions and in the presence of CNO. Computing the
receptive fields of all cell types (especially RGCs) allows us to to
disentangle the concerted effect of ACs lateral connectivity and CNO on the RF
of RGCs and provides an excellent agreement with experimental data. We argue,
on the basis of the model and analytical computations, that the BCs-ACs
network shapes the RF of RGCs via two main reduced dimensionless parameters,
one characterizing the ratio in the time scales of ACs and BCs response, and
the other, characterizing the ratio between the synaptic weight BCs $\to$ ACs
and ACs $\to$ BCs. This defines a two dimensional map, called the "RFs map".
We show how the experimental data settle in this map and we link the observed
changes in the shape of the RF, when applying CNO, to a displacement in this
map. Finally, we extrapolate our model predictions to situations where the
RGCs response could not be explained by the physiological properties of the
cell but would involve the ACs network.

##  3 Methods

###  3.1 The retina model

####  3.1.1 Structure

We assimilate the retina to a superimposition of $3$ layers, each one
corresponding to a cell type (BCs, ACs, RGCs), and being a flat, two-
dimensional square of edge length $L$ mm where spatial coordinates are noted
$x,y$ (Figure 1).

We consider a simplified form of connectivity, inspired from the real
connectivity in the retina, illustrated in Figure 1. First, there are as many
BCs as ACs and RGCs ($N$ cells per type so that the total number of cells is
$3N$). BCs are labelled with an index $i=1\dots N$, ACs with an index
$j=1\dots N$, RGCs with an index $k=1\dots N$. BCs are connected to ACs. We
note $W^{{A}_{j}}_{{B}_{i}}\leq 0$ the weight of the inhibitory synapse from
AC $j$ to BC $i$. We use the convention that $W^{{A}_{j}}_{{B}_{i}}=0$ when AC
$j$ is not connected to BC $i$. Likewise, we note $W^{{B}_{i}}_{{A}_{j}}\geq
0$ the weight of the excitatory synapse from BC $i$ to AC $j$,
$W^{{B}_{i}}_{{G}_{k}}\geq 0$ the weight of the excitatory synapse from BC $i$
to RGC $k$ and $W^{{A}_{j}}_{{G}_{k}}\leq 0$ the weight of the inhibitory
synapse from AC $j$ to RGC $k$. We note $W^{{A}}_{{B}}$ the (square) matrix of
connections from ACs to BCs and so on. Note that the model affords as well
connections between the same cell types (e.g. ACs to ACs or RGCs to RGCs via
gap junctions) but it would consequently complicate the theoretical analysis.

Figure 1: Translation of the retinal circuit to a computational network model.
Left. Schematic of the retina. Light activates the photoreceptor cells (PRs),
that transduce the input into a cascade of biochemical and electrical events
that can stimulate BCs and onwards RGCs. This vertical excitatory pathway is
modulated by inhibitory interneurons comprising two groups; horizontal (HCs)
and amacrine cells (ACs). All these neural signals are integrated by RGCs and
finally converted into action potentials going to the brain. Right. Schematic
view of the model. The joint integration of PRs and HCs in a limited region of
space is modelled by a spatio-temporal kernel mathematically defining the OPL
input to a BC, corresponding to eq. (1). The convolution of this kernel with
the stimulus drives the BC evolution, modulated by inhibitory connections with
ACs. BCs indeed make excitatory synaptic connections with ACs and ACs inhibit
BCs. Finally, RGCs pool over many BCs and ACs in their neighbourhood. Note
that the cells connections on the right panel are just shown for illustration
and do not necessarily correspond to the connectivity used in the model. Cells
are organised along a square lattice viewed here in a 3D perspective.

####  3.1.2 Visual input

Each neuron has a receptive field (RF), a specific region of the visual field
where light stimulus will modify the activity of this neuron. The term RF is
not limited only to the spatial region but it is often extended to include the
temporal structure within this region. The RF usually exhibits a centre-
surround organisation shared by many cell types in the retina (Masland, 2012a,
). Each BC receives synaptic inputs from its upstream circuitry, a combination
of dendritic excitatory inputs from photoreceptors (rods and cones) and
inhibitory inputs from horizontal cells (Franke and Baden, , 2017) occurring
at the level of the Outer Plexiform Layer (OPL) and this interaction emerges
on their RF.

A popular approach to simplify the complex process involved is to model the RF
as a single spatio-temporal linear filter that essentially represents the
opposition between the centre of the receptive field, driven by
photoreceptors, and the surround signal transmitted by horizontal cells. The
membrane potential of the BC’s soma can then be linearly approximated by the
convolution of the spatio-temporal kernel ${\mathcal{K}}_{{B}_{i}}$, featuring
the biophysical processes at the OPL, with the visual stimulus
${\mathcal{S}}(x,y,t)$. As we do not consider color sensitivity here,
${\mathcal{S}}$ characterizes a black and white scene, with a control on the
level of contrast $\in[0,1]$. The voltage of BC $i$ is stimulus-driven by the
term ("OPL input"):

$\begin{split}{\mathcal{V}}^{(drive)}_{{i}}(t)&=\left[\,{\mathcal{K}}_{{B}_{i}}\stackrel{{\scriptstyle
x,y,t}}{{\ast}}{\mathcal{S}}\,\right](t)\\\
&=\int_{x=-\infty}^{+\infty}\,\int_{y=-\infty}^{+\infty}\,\int_{s=-\infty}^{t}\,{\mathcal{K}}(x-x_{i},y-y_{i},t-s)\,{\mathcal{S}}(x,y,s)\,dx\,dy\,ds,\end{split}$
(1)

where $\stackrel{{\scriptstyle x,y,t}}{{\ast}}$ means space-time convolution.
For simplicity we consider only one family of BCs so that the kernel
${\mathcal{K}}$ is the same for all BCs. What changes is the center of the RF,
located at $x_{i},y_{i}$, which also corresponds to the coordinates of the BC
$i$. We consider in the paper a separable kernel
${\mathcal{K}}(x,y,t)={\mathcal{K}}_{S}(x,y)\,{\mathcal{K}}_{T}(t)$ where
${\mathcal{K}}_{S}$ is the spatial part and ${\mathcal{K}}_{T}$ the temporal
part. We restrict to ON or OFF BCs with a monophasic spatio-temporal kernel of
the form:

${\mathcal{K}}_{{}_{T}}=\left(\,A_{0}\,\frac{t^{2}}{2\,\tau_{RF}^{3}}e^{-\frac{t}{\tau_{RF}}}\,+\,b_{0}\,\right)\,H(t).$
(2)

where $H(t)$ is the Heaviside function. The parameter $A_{0}$ controls the
amplitude of the OPL input and can be positive (ON BC) or negative (OFF BC),
while $b_{0}$ controls the level of residual polarisation observed in
experiments, when light stimulation has stopped. This parameter is further
explained below. Experimentally, the BC temporal response can be biphasic in
time, because of receptor desensitization and because HCs feedback on to PRs
terminals. Note however, that the BCs response can also be monophasic
(Schreyer, , 2018). As we precisely want to point out the effect of lateral
inhibition on the BCs response, we consider a monophasic temporal OPL input
and show that, tuning the lateral connectivity parameters in the model, one
can indeed obtain biphasic shapes, as expected, but also more complex
dynamical responses.

The spatial part, ${\mathcal{K}}_{{}_{S}}$, is a classical difference of
Gaussians. This spatio temporal kernel is illustrated in Figure 2.

Figure 2: OPL kernel. Left. Spatial part in one spatial dimension. Right.
Temporal part.

####  3.1.3 Voltage dynamics

In the model, neurons are characterized by their voltage. Although RGCs are
spiking, we will be indeed interested only in their voltage variations, as a
function of the stimulus and network effects. Spiking could be obtained from
voltage e.g. using a Linear Nonlinear Poisson process (LNP) for spiking
neurons (Berry et al., , 1999; Chen et al., 2013a, ; Souihel and Cessac, ,
2021). We note $V_{{B}_{i}}$ the voltage of BC $i$, $V_{{A}_{j}}$ the voltage
of AC $j$, $V_{{G}_{k}}$ the voltage of RGC $k$. The joint dynamics of
voltages is given by the dynamical system:

$\left\\{\begin{array}[]{lll}\frac{dV_{{B}_{i}}}{dt}&=&-\frac{1}{\tau_{B_{i}}}V_{{B}_{i}}+\sum_{j=1}^{N_{A}}W^{{A}_{j}}_{{B}_{i}}\,{\mathcal{N}}^{(A)}\left(\,V_{{A}_{j}}\,\right)+F_{{B}_{i}}(t),\quad
i=1\dots N\\\ &&\\\
\frac{dV_{{A}_{j}}}{dt}&=&-\frac{1}{\tau_{A_{j}}}V_{{A}_{j}}+\sum_{i=1}^{N_{B}}W^{{B}_{i}}_{{A}_{j}}\,{\mathcal{N}}^{(B)}\left(\,V_{{B}_{i}}\,\right)+\zeta_{{A}},\quad
j=1\dots N,\\\ &&\\\
\frac{dV_{{G}_{k}}}{dt}&=&-\frac{1}{\tau_{G_{k}}}V_{{G}_{k}}+\sum_{i=1}^{N_{B}}W^{{B}_{i}}_{{G}_{k}}{\mathcal{N}}_{B}(V_{{B}_{i}})+\sum_{j=1}^{N_{A}}W^{{A}_{j}}_{{G}_{k}}{\mathcal{N}}_{A}(V_{{A}_{j}})+\zeta_{{G}},\quad
k=1\dots N.\end{array}\right.$ (3)

This form is derived, from first principles and using simplifying hypotheses,
in the section 6.1 of the supplementary material. The rectification term:

${\mathcal{N}}_{A}(V_{{A}_{j}})=\left\\{\begin{array}[]{lll}V_{{A}_{j}}-\theta_{A},&\mbox{if}\,\,V_{{A}_{j}}>\theta_{A};\\\
0,&\mbox{otherwise}\end{array}\right.,$ (4)

ensures that the synapse $j\to i$ is inactive when the pre-synaptic voltage
$V_{{A}_{j}}$ is smaller than the threshold $\theta_{A}$. The same holds for
the rectification term ${\mathcal{N}}_{B}(V_{{B}_{i}})$ with a threshold
$\theta_{B}$.

The term:

$F_{{B}_{i}}(t)=\frac{{\mathcal{V}}^{(drive)}_{{i}}}{\tau_{B_{i}}}\,+\,\frac{d}{dt}{\mathcal{V}}^{(drive)}_{{i}},$
(5)

is the stimulus driven input to BCs. It takes this form to ensure that, in the
absence of ACs coupling, $V_{B_{i}}(t)={\mathcal{V}}^{(drive)}_{{i}}(t)$.

CNO binds to specific designed receptors leading either to an excitatory or
inhibitory response, depending on the receptors ($G_{i}$ or $G_{q}$) (Roth, ,
2016; Urban and Roth, , 2015). Here, we propose a simplification of the actual
biological process for mathematical convenience. We model the CNO effect by a
current of the form
$-g_{CNO_{T}}\,\left(\,V-{\mathcal{E}}_{CNO_{T}}\,\right)$, where $V$ is the
voltage of a cell sensitive to CNO, $T$ is the cell type (i.e. only ACs or
RGCs according to the experimental set-up); $g_{CNO_{T}}$ is the conductance
of channels sensitive to CNO that is zero in the absence of CNO while it
increases with CNO concentration and ${\mathcal{E}}_{CNO_{T}}$ is the
corresponding reversal potential, positive for excitatory CNO and negative for
inhibitory.

As CNO modifies the membrane conductance it also induces a change of
polarization, which is characterized by the parameter $\zeta_{A}$ for ACs and
$\zeta_{G}$ for RGCs, with a general form (see section 6.1):

$\zeta_{T}=\frac{{\mathcal{E}}_{CNO_{T}}}{C}\,g_{CNO_{T}}$ (6)

where $T=A,G$ and where $C$ is the membrane capacitance assumed to be the same
for all cells.

In eq. (3), $\tau_{B_{i}},\tau_{A_{j}},\tau_{G_{k}}$ are the characteristic
integration times of the BC $i$, AC $j$, RGC $k$. As discussed in the next
section and in section 6.2 (supplementary), they depend on CNO in a complex,
non linear way.

####  3.1.4 Parameters reduction

The model depends on many parameters that constrain the dynamical evolution of
the system (3). Although some of our results, like the explicit form for RFs
(15), are quite general, it is easier, for analytic derivations and for
simulations to further simplify these parameters.

First, we consider that the synaptic weights from ACs to BCs, are controlled
by a unique parameter, $w^{-}>0$, $W^{{A}}_{{B}}=-w^{-}\,\Gamma^{{A}}_{{B}}$
(inhibitory synapse) where $\Gamma^{{A}}_{{B}}$ is an adjacency matrix
($\Gamma^{{A}_{j}}_{{B}_{i}}=1$ if the $j$-th AC connects to the $i$-th BC,
and $\Gamma^{{A}_{j}}_{{B}_{i}}=0$ otherwise). Likewise,
$W^{{B}}_{{A}}=w^{+}\,\Gamma^{{B}}_{{A}}$ with $w^{+}>0$ (excitatory synapse).
We use a simple form of connectivity, where
$\Gamma^{{B}}_{{A}}=\Gamma^{{A}}_{{B}}$ are nearest neighbors adjacency
matrices with null boundary conditions. The retina is regularly tiled by
different cell types so this approximation is reasonable, although in reality
one cell connects to more than 4 neighbors. Our "cells" must actually be
considered as "effective" cells with "effective" interactions. In particular,
our parameters $w^{-},w^{+}$ correspond to many synaptic contacts. This
simplification allows to considerably reduce the number of parameters shaping
the synaptic interactions and affords analytic computations (see
supplementary).

Following e.g. (Berry et al., , 1999; Chen et al., 2013a, ; Souihel and
Cessac, , 2021), the synaptic weight matrices:

$W^{{B}_{i}}_{{G}_{k}}=w^{B}_{G}\,\,\frac{e^{-\frac{d^{2}\left[\,{B}_{i},\,{G}_{k}\,\right]}{2\,\sigma_{p}^{2}}}}{\sqrt{2\pi}\,\sigma_{p}},\quad
W^{{A}_{j}}_{{G}_{k}}=w^{B}_{G}\,\,\frac{e^{-\frac{d^{2}\left[\,{A}_{j},\,{G}_{k}\,\right]}{2\,\sigma_{p}^{2}}}}{\sqrt{2\pi}\,\sigma_{p}},$
(7)

are Gaussian "pooling" matrices with $w^{B}_{G}>0$, $w^{A}_{G}<0$, and with
$d\left[\,{B}_{i},\,{G}_{k}\,\right]$, the $2$-d Euclidean distance between BC
$i$ and RGC $k$. The pooling standard deviation, $\sigma_{p}$ was fixed to
adapt to the spatial RF of observed cells (see section 4.1.2).

In this context, we can further understand the effect of CNO on the model. We
assume, for simplicity, that the lattice is so large that we can neglect
boundary effects, so that the characteristic times of BCs, ACs, RGCs and their
rest state noted $V_{B}^{\ast}$ and $V_{A}^{\ast}$, $V_{G}^{\ast}$, are
uniform in space (the general case is discussed in section 6.2). We set
$M^{{B}}_{{G}}=\sum_{i=1}^{N}W^{{B}_{i}}_{{G}_{k}}$ (resp.
$M^{{A}}_{{G}}=\sum_{j=1}^{N}W^{{A}_{j}}_{{G}_{k}}$). Strictly speaking, these
numbers depend on $k$, due to boundary conditions, but we neglect here this
dependence.

The characteristic times are then given by:

$\left\\{\begin{array}[]{lll}\tau_{B}&=&\frac{\tau_{L}}{1-\frac{2dw^{-}\tau_{L}}{{\mathcal{E}}_{A}}(V_{A}^{\ast}-\theta_{A})}\\\
\tau_{A}&=&\frac{\tau_{L}}{1+\tau_{L}\left(\,\frac{\zeta_{A}}{{\mathcal{E}}_{CNO_{A}}}+\frac{2dw^{+}}{{\mathcal{E}}_{B}}(V_{B}^{\ast}-\theta_{B})\,\right)}\\\
\tau_{G}&=&\frac{\tau_{L}}{1+\tau_{L}\left(\,\frac{\zeta_{G}}{{\mathcal{E}}_{CNO_{G}}}\,+\,\frac{M^{{B}}_{{G}}}{{\mathcal{E}}_{B}}\,\left(\,V^{\ast}_{B}\,-\,\theta_{B}\,\right)\,+\,\frac{M^{{A}}_{{G}}}{{\mathcal{E}}_{A}}\,\left(\,V^{\ast}_{A}\,-\,\theta_{A}\,\right)\,\right)},\end{array}\right.$
(8)

where $d$ is the lattice dimension, $\tau_{L}$ the leak characteristic time,
${\mathcal{E}}_{B}$ and ${\mathcal{E}}_{A}$ are respectively the reversal
potential for the synaptic connection from BCs to ACs, and for ACs to BCs (see
sections 6.1 and 6.2 for detail). The rest states are given by:

$\left\\{\begin{array}[]{lll}V^{\ast}_{B}&=&\frac{4d^{2}\tau_{A}\tau_{B}w^{-}w^{+}\theta_{B}-2d\tau_{B}w^{-}\tau_{A}\,\zeta_{{A}}\,+\,2d\tau_{B}w^{-}\theta_{A}}{1+4d^{2}\tau_{A}\tau_{B}w^{-}w^{+}}\\\
V^{\ast}_{A}&=&\frac{4d^{2}\tau_{A}\tau_{B}w^{-}w^{+}\theta_{A}\,-\,2d\tau_{A}w^{+}\theta_{B}\,+\,\tau_{A}\,\zeta_{{A}}}{1+4d^{2}\tau_{A}\tau_{B}w^{-}w^{+}}\\\
V^{\ast}_{G}&=&\tau_{G}\left[\,M^{{B}}_{{G}}\,\left(\,V^{\ast}_{B}\,-\,\theta_{B}\,\right)\,+\,M^{{A}}_{{G}}\,\left(\,V^{\ast}_{A}\,-\,\theta_{A}\,\right)\,+\,\zeta_{{G}}\,\right].\end{array}\right.$
(9)

Equations (8) and (9) hold if
$V^{\ast}_{B}>\theta_{B},V^{\ast}_{A}>\theta_{A}$, i.e. the rest state is not
rectified. This entails conditions on $\zeta_{A}$ that we computed explicitly.
However, these computations where too lengthy to be included in the paper. The
interested reader will find them on the web page
https://team.inria.fr/biovision/mathematicalderivationscno/. In particular,
considering a positive $\zeta_{A}$ (excitatory CNO), having
$V^{\ast}_{B}>\theta_{B},V^{\ast}_{A}>\theta_{A}$ require that $\theta_{B}<0$
and that it is compatible with $\underline{\theta_{A}=0}$, a condition that we
will assume from now on.

On this web page, we also address the question of the variation of
$\tau_{B},\tau_{A},V^{\ast}_{B},V^{\ast}_{A}$ as $\zeta_{A}$ varies. We found
that $\frac{\partial
V_{B}^{\ast}}{\partial\zeta_{A}}=-\frac{2d\tau_{A}\tau_{B}w^{-}}{1+4d^{2}\tau_{A}\tau_{B}w^{-}w^{+}}\leq
0$, meaning that BCs get more hyperpolarized when the excitatory effect of CNO
on ACs increases, while $\frac{\partial
V_{A}^{\ast}}{\partial\zeta_{A}}=\frac{\tau_{A}}{1+4d^{2}\tau_{A}\tau_{B}w^{-}w^{+}}\geq
0$ (ACs get more depolarized when $\zeta_{A}$ increases). This result is the
one qualitatively expected, but note that our modelling is also quantitative.
In strong contrast, the computation of
$\frac{\partial\tau_{B}}{\partial\zeta_{A}}$ and
$\frac{\partial\tau_{A}}{\partial\zeta_{A}}$ reveals that the sign of these
derivative depends on specific domains in the space of parameters
$w^{-},w^{+},\theta_{B}$, with nonlinear frontiers. That is, even for this
reduced model, the behaviour of the characteristic times as CNO increase
depends on network parameters in a non trivial way. There are actually domains
where these characteristic times increase with $\zeta_{A}$.

Thus, while a rapid argumentation would predict that, the characteristic times
of ACs should increase because the membrane conductance increases as
$\zeta_{A}$ increases, the opposite can happen as well. This is because, as
revealed by eq. (8), the characteristic time $\tau_{A}$ depends on the rest
states (eq. (9) ) which depend themselves on $\zeta_{A}$ via the
characteristic times $\tau_{A},\tau_{B}$. This results in a two dimensional
implicit non linear system studied on the web page
https://team.inria.fr/biovision/mathematicalderivationscno/.

As explained in the introduction, our main goal is to confront this modelling
to experiments and extract general statements from it. One major difficulty,
however, is that we cannot fine tune the parameter $\zeta_{A}$ (the CNO
conductance) in experiments because it is not possible to establish a CNO dose
response curve of the changes in retinal activity. Therefore, even under
modelling simplifications it is not possible to fit the model from experiments
with the non linear equations (8) and (9). When fitting experimental results
we computed the characteristic integration times of cells in control (CTL) and
CNO conditions independently (see section 4.1.2). Likewise, we cannot access
the rest state of BCs and ACs. We fixed the thresholds $\theta_{A},\theta_{B}$
to zero. (Mathematically, $\theta_{B}$ should be strictly negative, but we can
consider a very small absolute value.) Thus, the rest states (9) are zero in
the absence of CNO. Nevertheless, in experiments, one can observe a small
residual polarisation, which changes in the presence of CNO. That is why we
added the polarization parameter, $b_{0}$, appearing in eq. (2), a huge, but
necessary simplification of the entangled equations (9).

Most of the analysis below will be done considering that no rectification
takes place, even under light simulation, so that we essentially consider a
linear model. This is supported by the paper of Baccus et al.,  (2008) showing
that the voltage of bipolar cells is essentially a linear function of the
stimulus for white noise. For a more general analysis dealing with
rectification please check the supplementary, as well as the discussion
section.

####  3.1.5 Experimental set-up

In our experiments, excitatory DREADDs (hM3Dq) were activated using CNO on
RGCs and ACs co-expressing a certain gene (Scnn1a or Grik4), triggering a
calcium release from organelles and thus, leading to increase of intracellular
concentration of free calcium. This resulted in membrane depolarisation and
higher neuronal excitability. Our experiments suggested that subclasses of ACs
and RGCs could be simultaneously sensitive to CNO but we did not observe any
evidence of an effect on BCs.

Detailed experimental details can be found in our recent publication (Hilgen
et al., , 2022). All experimental procedures were approved by the ethics
committee at Newcastle University and carried out in accordance with the
guidelines of the UK Home Office, under the control of the Animals (Scientific
Procedures) Act 1986. Recordings were performed on the BioCamX platform with
high-density-multielectrode array (HD-MEA) Arena chips (3Brain GmbH, Lanquart,
Switzerland), integrating 4096 square microelectrodes in a 2.67 × 2.67
$\text{\,}{\mathrm{mm}}^{2}$ area and aligned in a square grid with
$42\text{\,}\mathrm{\SIUnitSymbolMicro m}$ spacing. Light stimuli were
projected onto the retina using a LED projector. Briefly, the projector
irradiance was attenuated using neutral density filters to mesopic light
levels.

##  4 Results

In this section we present the theoretical and numerical results based on our
retina model. We provide only the main conclusions of the mathematical
derivations, which are presented in detail in the supplementary section.

###  4.1 Model fitting of ganglion cell receptive fields characterized from
experimental data

RGC responses emanate from a dynamic balance of synaptic excitation and
inhibition, originating from the interactions of BCs and ACs. We believe that
such network connectivity gives rise to various response patterns and we show
that our model can capture these joint effects, by providing an analytic form
of the RF of the cells. As we demonstrate, this computation provides us an
algorithmic way to fit the model parameters to the light responses recorded
from mouse RGCs. One can then infer the possible behaviour of ACs and BCs
leading to this RGC response, even if we do not measure them experimentally.

####  4.1.1 Mathematical form of the RF of retinal cells

The results presented below hold for all cell types. Thus, we label cells with
a generic index $\alpha=1\dots 3N$. BCs have an index $\alpha=1\dots N$, ACs
have an index $\alpha=N+1\dots 2N$, RGCs have an index $\alpha=2N+1\dots 3N$
and we write ${\mathcal{X}}_{\alpha}$ the voltage of cell $\alpha$.

The time evolution of the dynamical system in eq. (3) is controlled by a
matrix, ${\mathcal{L}}$, called "transport operator" and explicitly written in
the supplementary section 6.3. ${\mathcal{L}}$ depends on the connectivity
matrices $W^{{A}}_{{B}},W^{{B}}_{{A}},W^{{A}}_{{G}},W^{{B}}_{{G}}$ and on all
the parameters controlling the dynamics. The form of ${\mathcal{L}}$ also
depends on the set of rectified cells. In the following, we assume that cells
are not rectified i.e. hyperpolarised BCs do not reach the rectification
threshold (the rectified case is discussed in the conclusion section and in
the supplementary material). Consequently, the dynamical system (3) is linear.

In this case, the eigenvalues $\lambda_{\beta}$, $\beta=1\dots 3N$ and the
eigenvectors ${\mathcal{P}}_{\beta}$ of ${\mathcal{L}}$ characterize the
evolution of cells’ voltages. We note ${\mathcal{P}}$ the matrix transforming
${\mathcal{L}}$ in diagonal form (the columns of ${\mathcal{P}}$ are the
eigenvectors ${\mathcal{P}}_{\beta}$) and ${\mathcal{P}}^{-1}$ its inverse.

In this context, we show in the supplementary section 6.3.3 that the voltage
of a cell with index $\alpha$ is the sum of $4$ terms:

${\mathcal{X}}_{\alpha}(t)={\mathcal{V}}^{(drive)}_{{\alpha}}(t)\,+{\mathcal{E}}^{(drive)}_{{\alpha}}(t)\,+\,{\mathcal{E}}^{(CNO_{A})}_{{\alpha}}\,+\,{\mathcal{E}}^{(CNO_{G})}_{{\alpha}},\quad\alpha=1\dots
3N.$ (10)

##### Stimulus drive.

The first term, ${\mathcal{V}}^{(drive)}_{{\alpha}}(t)$ corresponds to (1),
and is non zero for BCs only. It corresponds to the BCs response in the
absence of the ACs network.

##### CNO effects.

The terms:

${\mathcal{E}}^{(CNO_{A})}_{{\alpha}}=\zeta_{A}\,\sum_{\beta=1}^{3N}\sum_{\gamma=N+1}^{2N}\frac{{\mathcal{P}}_{\alpha\beta}\,{\mathcal{P}}^{-1}_{\beta\gamma}}{\lambda_{\beta}},\quad\alpha=N+1\dots
2N;$ (11)

and

${\mathcal{E}}^{(CNO_{G})}_{{\alpha}}=\zeta_{G}\,\tau_{G_{\alpha}},\quad\alpha=2N+1\dots
3N.$ (12)

correspond, respectively, to the impact of CNO on the voltages of ACs and
RGCs. There are important nonlinear effects hidden in the terms
$\frac{{\mathcal{P}}_{\alpha\beta}\,{\mathcal{P}}^{-1}_{\beta\gamma}}{\lambda_{\beta}}$
(eq. (11)) which depend themselves on the characteristic times (8). Thus, the
polarization level of ACs and RGCs is not only fixed by the direct effect of
CNO on the cell, but is also tuned by entangled network effects.

##### Stimulus-network interaction term.

In eq. (10), the term :

${\mathcal{E}}^{(drive)}_{{\alpha}}(t)=\sum_{\beta=1}^{3N}\sum_{\gamma=1}^{N}{\mathcal{P}}_{\alpha\beta}\,{\mathcal{P}}^{-1}_{\beta\gamma}\varpi_{\beta\gamma}\,\int_{-\infty}^{t}e^{\lambda_{\beta}(t-s)}\,{\mathcal{V}}^{(drive)}_{{\gamma}}(s)\,ds,\quad\alpha=1\dots
3N,$ (13)

where $\varpi_{\beta\gamma}=\frac{1}{\tau_{B_{\gamma}}}+\lambda_{\beta}$,
corresponds to the indirect effect, via the network connectivity, of the
stimulus drive on (i) BCs, for $\alpha=1\dots N$; (ii) ACs for
$\alpha=N+1\dots 2N$; (iii) RGCs $\alpha=2N+1\dots 3N$. Thus, this equation
describes the response of all cells to the stimulus. Especially, it tells us
how the direct input (1) to BCs is modulated by the concerted activity of BCs
and ACs.

Mathematically, the term (13) can be interpreted as follows. The drive (index
$\gamma=1\dots N$) triggers the eigenmodes $\beta=1\dots 3N$ of
${\mathcal{L}}$, with a weight proportional to
${\mathcal{P}}^{-1}_{\beta\gamma}$. The mode $\beta$, in turn, acts on the
voltage of cell $\alpha=1\dots 3N$ with a weight proportional to
${\mathcal{P}}_{\alpha\beta}$. The time dependence and the effect of the drive
are controlled by the integral
$\int_{-\infty}^{t}e^{\lambda_{\beta}(t-s)}\,{\mathcal{V}}^{(drive)}_{{\gamma}}(s)\,ds$.
The eigenvalues $\lambda_{\beta}$ introduce characteristic time scales in the
response: exponential decay rate for the real part, frequency of oscillations
for the imaginary part, if any. The right eigenvectors (columns of
${\mathcal{P}}$) and the left eigenvectors (rows of ${\mathcal{P}}^{-1}$)
introduce characteristic space scales in the response. The easiest case to
figure this out is the nearest neighbours connectivity considered in section
5.3.5. Here, eigenvectors are labelled by an index $n$ which actually
corresponds to a wave vector, (eq. (49)), the inverse of a space scale. In
this context, the eigenmode with $n=1$ has the largest space scale (the size
of the lattice) whereas the eigenmode with $n=N$ corresponds to the minimum
space scale (the distance between two nodes). Thus, eq. (13) illustrates that
the response to a spatio-temporal stimulus is a combination of multiple time
scales and space scales, leading to interesting phenomena such as resonances
as commented below, or waves of anticipation as studied in Souihel and Cessac,
(2021).

##### The Receptive Field of all cell types.

Introducing the function $e_{\beta}(t)\equiv e^{\lambda_{\beta}\,t}\,H(t)$ so
that
$\int_{-\infty}^{t}e^{\lambda_{\beta}(t-s)}\,{\mathcal{V}}^{(drive)}_{{\gamma}}(s)\,ds\equiv\left[\,e_{\beta}\stackrel{{\scriptstyle
t}}{{\ast}}{\mathcal{K}}_{{B}_{T}}\stackrel{{\scriptstyle
t}}{{\ast}}\left(\,{\mathcal{K}}_{{B}_{S_{\gamma}}}\stackrel{{\scriptstyle
x,y}}{{\ast}}{\mathcal{S}}\,\right)\,\right](t)$, and using the separated
kernel form (2), the response (13) reads:

${\mathcal{E}}^{(drive)}_{{\alpha}}(t)=\left[\,{\mathcal{K}}_{{}_{\alpha}}\stackrel{{\scriptstyle
x,y,t}}{{\ast}}S\,\right](t),$ (14)

with:

${\mathcal{K}}_{{}_{\alpha}}(x,y,t)=\sum_{\beta=1}^{3N}\left(\,{\mathcal{P}}_{\alpha\beta}\,U_{\beta}(t)\,\times\,\sum_{\gamma=1}^{N}{\mathcal{P}}^{-1}_{\beta\gamma}\,\varpi_{\beta\gamma}\,{\mathcal{K}}_{{B}_{S_{\gamma}}}(x,y)\,\right),$
(15)

where we have set $U_{\beta}(t)\equiv\left[\,e_{\beta}\stackrel{{\scriptstyle
t}}{{\ast}}{\mathcal{K}}_{{B}_{T}}\,\right](t)$. The response of cell $\alpha$
is thus expressed as a convolution of the stimulus with a spatio-temporal
kernel ${\mathcal{K}}_{{}_{\alpha}}(x,y,t)$, an expected result from the
linear response. Nevertheless, it’s important to point out that the expression
(15) holds for all cell types, not only RGCs and that it contains the network
effects induced by the BCs-ACs network. Thus, for $\alpha=1\dots N$, equation
(15) characterizes the indirect (network induced) response of BCs to the
stimulus drive, in addition to the direct response (1). For $\alpha=N+1\dots
2N$, equation (15) represents the RF of ACs. Finally, for $\alpha=2N+1\dots
3N$ we obtain the RF of RGCs. We focus on this last case from now on,
essentially because this predicted RF can be confronted to experiments,
whereas we have no experimental access to BCs or ACs RF.

##### The Receptive Field of RGCs.

Henceforth, we will refer to ${\mathcal{K}}_{{}_{\alpha}}(t)$ as
${\mathcal{K}}_{{G}_{\alpha}}(t)$, to make explicit that we are dealing with
RGCs. The RF of RGCs can be often written as a product of a space dependent
term and a time dependent term (separability). In our case, this would
correspond to write ${\mathcal{K}}_{{G}_{\alpha}}(x,y,t)$ in the form of a
product
${\mathcal{K}}_{{G}_{\alpha}}(x,y,t)={\mathcal{K}}_{{G}_{T_{\alpha}}}(t)\times{\mathcal{K}}_{{G}_{S_{\alpha}}}(x,y)$
where:

${\mathcal{K}}_{{G}_{T_{\alpha}}}(t)=\sum_{\beta=1}^{3N}{\mathcal{P}}_{\alpha\beta}\,U_{\beta}(t),$
(16)

is a temporal kernel and:

${\mathcal{K}}_{{G}_{S_{\alpha}(x,y)}}=\sum_{\gamma=1}^{N}{\mathcal{P}}^{-1}_{\beta\gamma}\,\varpi_{\beta\gamma}\,{\mathcal{K}}_{{B}_{S_{\gamma}}}(x,y),$
(17)

is a spatial kernel.

This separation is not strictly possible in eq. (15), because there is a
dependency on $\beta$ on the term
$\sum_{\gamma=1}^{N}{\mathcal{P}}^{-1}_{\beta\gamma}\,\varpi_{\beta\gamma}\,{\mathcal{K}}_{{B}_{S_{\gamma}}}(x,y)$.
Thus, from this computation the RF of RGCs is not expected to be separable in
general. However, it is important to remark that the stimulus used in
experiments for determining RFs, the white noise, is uniform in probability in
the space domain and in the time domain. This induces an effective separation
of the response which might not hold for more complex stimuli (e.g. moving
objects). As the experimental RFs considered here have been obtained from
white noise we will assume separability from now on.

##### The spatial part of the RGCs RF.

Equation (17) appears as an overlap of spatial RFs of BCs. In such naive
overlaps approximations, spatial RFs of BCs are just summed up with a uniform
weight. However, here the RF of each BC is constrained by ACs lateral
connectivity, via the term
${\mathcal{P}}^{-1}_{\beta\gamma}\,\varpi_{\beta\gamma}$. In particular,
equation (17) is not necessarily circular even if BCs RFs are, and the center
of the RGC cell RF is not necessarily at the barycentre of connected BCs RFs.
This holds, for example, if AC connectivity is not invariant by rotation.

##### The temporal part of the RGCs RF,

(16). As we consider monophasic temporal kernels ${\mathcal{K}}_{{B}_{T}}$ of
BCs with the form (2) we have:

$U_{\beta}(t)=A_{0}\,\left(\,\frac{\left[\,2\tau_{RF}^{2}\,e^{\lambda_{\beta}\,t}-\left(\,t^{2}\,\left(\,\lambda_{\beta}\tau_{RF}+1\,\right)^{2}+2t\tau_{RF}\left(\,\lambda_{\beta}\tau_{RF}+1\,\right)+2\tau_{RF}^{2}\,\right)\,e^{-\frac{t}{\tau_{RF}}}\,\right]}{2\,\left(\,\lambda_{\beta}\tau_{RF}+1\,\right)^{3}\,\tau_{RF}^{2}}\,\right)\,H(t)$

(18)

$U_{\beta}$, and, thereby, ${\mathcal{K}}_{{G}_{T_{\alpha}}}(t)$ is the
temporal part of the RGC receptive field, that changes their shape due to
variations in the eigenvalues of ${\mathcal{L}}$, who are themselves
controlled by model parameters. A striking effect arises when some eigenvalues
become complex, leading to temporal oscillations of $U_{\beta}$. This remark
is at the core of the analysis exposed in section 4.2.2.

####  4.1.2 Fitting the RFs of ganglion cells

Experimentally, RGCs RFs were reconstructed from Spike Triggered Average (STA)
in response to Shifted White Noise (SWN). The SWN, introduced by Hilgen et
al.,  (2017); Pamplona et al.,  (2021), is a spatially uniform noise where the
images of random black and white squares in the classical "White Noise"
stimulus have, in addition, random spatial shifts, improving the resolution of
the STA. This allowed us to fit the spatial and temporal part of the RF. Note,
however, that it is difficult to obtain the surround in the spatial part. From
the center part of the spatial RF, we fixed the parameter $\sigma_{p}$ in the
Gaussian pooling, eq. (7). More extended results can be found in E. Kartsaki
thesis (Kartsaki, , 2022), where we display experimental spatial RFs, see e.g.
section 3.4.1 of the thesis.

We mainly focused on the temporal part of the RF. The time RF estimation
resulted in temporal traces with duration $600$ ms sampled with a rate
$33/4=8.25$ ms. In order to assess the validity of the model, we have fitted
these time traces in CTL and CNO conditions. We have a data base of $117$
cells sensitive to CNO, i.e. exhibiting increase or decrease in firing rate
beyond a certain threshold.

These reconstructed temporal RFs provide the linear response of a RGC to a
spatially uniform flashed stimulus, mathematically corresponding to a Dirac
distribution. As we wanted to compare our model’s output, the RGC voltage,
$V_{G}$, to this experimental RF, computed from firing rates, we neglected the
effect of non linearities and assumed that the experimental response is
proportional to the RGC voltage. We considered a one dimensional model (chain)
with $N=60$ cells of each type where the cells at the boundaries have a fixed,
zero, voltage (zero boundary conditions), corresponding to the reference rest
state. To reduce the boundaries effect, we made the fit for the RGC in the
center of the network.

We performed simulations of the model (3) using a spatially uniform Dirac
pulse as the stimulus and compute the cell responses, using two modalities:
simulation of the differential equations (3) (green traces labelled "Sim" in
Figure 3) and analytic computation (15) (black traces labelled "Th" in Figure
3). We observe that these two traces are always identical confirming the
goodness of the simulation scheme.

We recall that the parameters shaping this response are: $A_{0}$, the
intensity of the OPL input; $b_{0}$, controlling a residual
depolarization/hyper polarization observed in experimental responses; $w^{+}$,
controlling the synaptic intensity from BCs to ACs; $w^{-}$, controlling the
synaptic intensity from ACs to BCs; $w^{B}_{G}$, controlling the synaptic
intensity from BCs to RGCs; $w^{A}_{G}$, controlling the synaptic intensity
from ACs to RGCs; $\tau_{B}$, the characteristic membrane time scale of BCs;
$\tau_{A}$, the characteristic membrane time scale of ACs; $\tau_{RF}$, the
characteristic membrane time scale of the OPL drive; $\tau_{G}$, the
characteristic membrane time scale of RGCs.

We note $\vec{\eta}$ the set of all the parameters shaping response.
$\vec{\eta}$ is therefore a point in a $10$-dimensional space.

The fit was then done by a gradient descent to minimize the $L^{2}$-distance
$D_{2}(\vec{\eta})$ between the experimental trace of the time STA, $STA(s)$
and the theoretical temporal RF (16) which depends on $\vec{\eta}$. The
minimization is done by iterating the differential equation:

$\frac{d\vec{\eta}}{du}=-\vec{\nabla}_{\vec{\eta}}D_{2}.$

The gradient of $\vec{\nabla}_{\vec{\eta}}D_{2}$ involves
$\vec{\nabla}_{\vec{\eta}}{\mathcal{K}}_{{G}_{T_{\alpha}}}(s)$ which can be
explicitly computed when we have the analytic form of $RF$, or numerically.
Note that having the analytic form gives better results especially because it
allows second order corrections (Hessian). Our minimisation is done only on
the temporal trace. Nevertheless, the simulation allows us to draw the
corresponding spatial RF.

Although the experimental temporal RFs were quite diverse among cells, we were
able to fit all of them with a very good accuracy (final error smaller than
$1\%$). We rejected fits where some parameters became unrealistic (e.g.
$\tau_{A}$ larger than $1$ s or $\left|\,w^{-}\,\right|>1$ kHz). We rejected
about $4\%$ of the fits. An example of fit is shown in Figure 3.

In this synthetic representation the simulated responses of the OPL, BCs and
ACs connected to the RGC located at the centre of the network appears in the
top left figure. We also show the simulated response of the RGC vs the
experimental temporal STA of this cell (bottom left). On the top right we see
the numerical spatio-temporal RF using a color map. Finally, at the bottom
right, we show the power spectrum (modulus of the Fourier transform) of the
temporal response. As developed below, this spectrum provides important
information on the cell response.

RGCs RFs in the presence of CNO could be fit equally well with our model. All
fitted cells, in CTL and CNO conditions, can be found on the web page
https://gitlab.inria.fr/biovision/dreadds. The C-code used for simulations can
also be found there.

Figure 3: Summary panel. The top left panel (a) illustrates the simulated
responses of the OPL (term ${\mathcal{V}}^{(drive)}_{{\alpha}}(t)$ in (1),
black trace), BCs (red trace), ACs (blue trace) connected to the RGC located
at the centre of the network. The bottom left panel (b) shows the simulated
response of the RGC (green, "Sim", and black, "Th", trace) vs the experimental
temporal STA of this cell (orange dots). The green trace ("Sim") is the result
of a numerical simulation of the dynamical system (3) under a spatially
uniform flashed stimulus, whereas the black trace ("Th") is the result given
by the analytic expression (15). The top right panel (c) shows the spatio-
temporal RF of the RGC, time in abscissa, space in ordinate. The bottom right
panel (d) displays the power spectrum of the time response, experimental
(orange dots) and theoretical (black lines).

###  4.2 Network connectivity shapes the receptive fields of ganglion cells

Throughout our analysis of the experimental data, we observed great
variability in the effect of DREADD activation with CNO on the RGCs responses.
In this section, we develop the consequences of our mathematical analysis in
an attempt to explain this observed diversity of CNO effects on RF features.
We propose here an explanation purely based on network effects. There are
certainly other possible interpretations based on single cell characteristics
such as non linear effects due to changes in conductance etc, discussed in the
discussion section. The main advantages of our analysis is that it determines
network effects on the RGCs RF, controlled by two main parameters, and that it
predicts the response to more complex stimuli than full field flashes.

####  4.2.1 Two main parameters constrain the RF shape of ganglion cells

The entangled, feedback effects of ACs-BCs can be characterized by two
dimensionless parameters. The first one, $r=\frac{\tau_{A}}{\tau_{B}}$,
characterizes the ratio between the ACs and BCs membrane integration times.
The second, $s=\frac{w^{-}}{w^{+}}$, characterizes the ratio between the ACs
$\to$ BCs interaction ($w^{-}$) and the BCs $\to$ ACs interaction ($w^{+}$).
Of course, the other parameters play an important role when fitting a specific
RF. But, what we argue here is that the shape of RF and its space-time scaling
essentially depend on the value of $r,s$.

The theoretical explanation is that the RF of a RGC is given by the formula
(14), which is a cascade of convolutions involving the BC response to the
stimulus (OPL input) and the network effects expressed in terms of eigenvalues
$\lambda_{\beta}$ and eigenvectors components ${\mathcal{P}}_{\alpha\beta}$
appearing in equation (15). As explained in the supplementary section 6.3.4,
these eigenvalues and eigenvectors are essentially tuned by the two parameters
$r,s$. There is also a dependence on other parameters discussed in section
4.2.3. Depending on the location in the space $r,s$, some eigenvalues are
real, some others are complex. All eigenvalues have a negative real part,
ensuring the stability of the linear system. Imaginary parts in eigenvalues
introduce oscillations in the response, whereas the real part fixes a
characteristic decay time. The RF formula (15), involving a sum of
exponentials $e^{\lambda_{\beta}\,t}$ mixes these effects. As we considered
monophasic OPL response here, the time RF of RGCs is monophasic when all
eigenvalues are real. In contrast, oscillations in this RF can appear when
some eigenvalues are complex. However, the shape of this RF depends in more
detail on the period of oscillations, brought by the imaginary part of complex
eigenvalues, and on the characteristic decay times, brought by the real parts.

When moving in the $(r,s)$ plane, the eigenvalue $n$ switches from real to
complex conjugate pair when crossing a critical line, depending on $n$, whose
equation (48) is given in the supplementary section 6.3.4. There are $2N$
eigenvalues associated with the BCs-ACs network each one determining a
critical line in the plane $(r,s)$. The set of all these lines is what we call
the "skeleton". An example of this skeleton is shown in Figure 4, where we
only show some of the critical lines. These lines delimit color regions
corresponding to the number of complex eigenvalues (see colorbar legend on the
right of the figure).

####  4.2.2 The RFs map

The existence of this skeleton determines regions in the $(r,s)$ plane with
specific shapes for the temporal RF, given by eq. (16), a linear combination
of functions $U_{\beta}(t)$ given by (18). The Fourier transform
$\hat{U}_{\beta}(\omega)$ of $U_{\beta}(t)$ is:

$\hat{U}_{\beta}(\omega)=\frac{1}{\left(\,1+i\omega\tau_{RF}\,\right)^{3}}\,\frac{1}{i\omega-\lambda_{\beta}}.$
(19)

Thus, the Fourier transform of $\eqref{eq:KGT}$ is:

${\mathcal{K}}_{{G}_{T_{\alpha}}}(\omega)=\sum_{\beta=1}^{3N}{\mathcal{P}}_{\alpha\beta}\,\hat{U}_{\beta}(\omega),$
(20)

a linear combination of rational fractions. Extending to complex $\omega$s,
$\hat{U}_{\beta}(\omega)$ has two poles: $\omega=\frac{i}{\tau_{RF}}$ and
$\omega=-i\,\lambda_{\beta}$, corresponding to complex resonances. The
contributions of all these poles (for $\beta=1\dots 3N$) are combined in eq.
(20) with weights ${\mathcal{P}}_{\alpha\beta}$.

As we move in the $(r,s)$ plane, we notice the following. When $(r,s)$ are
small, eigenvalues are real and the terms ${\mathcal{P}}_{\alpha\beta}$ are
close to diagonal. In this case, the dominant pole contribution in (20) is the
pole $\omega=i\frac{1}{\tau_{RF}}$ corresponding to the OPL contribution.
Equation (20) has a single peak centered at $\omega=0$, corresponding to a
monophasic response. For larger values of $r,s$ some eigenvalues become
complex, giving potential additional peaks in the power spectrum. Actually, we
observed two cases mutually compatible. First, the central peak at $\omega=0$
switches to a non zero value. This corresponds to the appearance of an
exponentially damped oscillation in the RF, giving a biphasic response.
However, secondary peaks may appear leading to residual oscillations, in
addition to the main trend (monophasic or biphasic). This gives what we call a
polyphasic response. Such residual oscillations were observed in our
experiments and were relatively numerous (about $40\%$). There are, of course,
other hypotheses explaining these residual oscillations, but here, we will
support the hypothesis that they are generated by a network effect. An example
is given in Figure 3 where we observe, at the bottom left, residual
oscillations after the main biphasic response, and, at the bottom right, the
power spectrum with a main peak not centered at zero and a secondary peak
corresponding to the residual oscillations. Note that this secondary peak is
observed on experimental data but we failed to reproduce it in the fits. This
is further explained in the discussion section.

This analysis leads us to broadly decompose the $(r,s)$ plane into $3$ regions
corresponding to cells response phases: monophasic, biphasic, polyphasic. One
switches from one phase to the other when some peaks in the power spectrum
appear or disappear, driven by the spectrum of ${\mathcal{L}}$.

The corresponding "phase diagram", obtained numerically, constitutes what we
call the "RFs map" shown in Figure 4 (top right). There are four points,
labelled A,B,C,D, on this map, each representing a different cell’s response
phase. For each point, we have plotted the RGC temporal RF, as computed with
the model (bottom panels). A more general representation of what is going on
when moving along a specific pathway in this map can be found at the web page
https://team.inria.fr/biovision/cno_paper_supplementary/, where one can see
movies showing how the network effects shape the RGCs RFs when $(r,s)$ vary.
It is important to note that this map is a projection from a $10$ dimensional
space in $2$ dimensions. It has been obtained, from the skeleton, by fixing
the other parameters values, based on the average values obtained from fits
over the cell populations (see e.g. Fig. 5).

Figure 4: Top Left. The skeleton. Map in the $(r,s)$ plane showing the
"skeleton" of eigenvalues structure (white lines). When crossing line $n$ the
eigenvalue $n$ switches from real to complex. The color corresponds to the
number of complex eigenvalues (see the color bar on the right side). Top
right. The RFs Map. This figures summarizes how the entangled effects of the
BCs-ACs network act on the RF of RGCs. In this map in the $(r,s)$ plane we
distinguish three main regions (see text for their determination): monophasic,
biphasic, polyphasic. The points labelled "A", "B, "C", "D" correspond to the
temporal RF plotted in the bottom panels respectively called "Panel A", "Panel
B", "Panel C", "Panel D". We use the same representation as in Figure 3.

Figure 4 illustrates how the BCs-ACs network shapes the RF of a RGC by a
subtle balance between BCs-ACs, BCs-ACs interactions (parameter s) and the
time scale of their response (parameter r). For simplicity, we consider here
ON BCs, but the explanation holds also for OFF BCs.

The OPL drive (black line) induces a depolarisation of BCs (red line) within a
time scale of order $\tau_{B}$. This excites the connected ACs (blue line),
with an intensity $w^{+}$. The excitation of ACs hyperpolarises BCs with an
intensity $w^{-}$ within a time scale of order $\tau_{A}$. RGCs receive a
combination of excitatory and inhibitory inputs from their afferent circuit
with respective weights $w^{B}_{G}$ and $w^{A}_{G}$.

In the monophasic region, ACs respond in the same time scale as BCs. One
observes, for RGCs a monophasic response (Panel "A", middle left) whose
intensity depends on the ratio between excitation, provided by BCs, with a
weight $w^{B}_{G}$ and inhibition, with a weight $w^{A}_{G}$. When moving to
the biphasic region, ACs respond slower than BCS, leading to the biphasic
response illustrated with panel "B". When moving upward in the RFs map
(increasing $s$) this biphasic response becomes polyphasic with oscillations.
This is illustrated in panel "C". BCs start to raise due to the OPL drive,
leading to a rising of ACs, slower than BCs, leading to a hyper-polarisation
of BCs. This leads to a decrease of ACs voltage, thereby, to a rising of BCs
which still respond to the OPL drive. This cycle can be repeated several
times, depending again on the parameters $r,s$. Note that the amplitude is
always decreasing exponentially fast. The period of the observed oscillations
and the damping characteristic time depend on the location in the Map.
Finally, panel "D" is the point at the intersection of the $3$ phases regions.
We show it for completeness. Note that increasing $s$ leads to a decrease of
the RGCs response. When $s$ is too high, the response becomes too weak to be
observed experimentally.

####  4.2.3 Experimental cells spread in the RFs map

To confront our theoretical insight with experiments, we have placed the
recorded cells in the RFs map as shown in Figure 5. That is, for each
experimentally recorded RGC, we fit the model parameters $\vec{\eta}$ as
explained in section 4.1.2 thereby providing an estimation of $r,s$,
independently in CTL and in CNO conditions. This defines a virtual network,
made of identical cells in each layer, where all RGCs are responding like the
experimental RGC. As mentioned above, the map is only a projection of
$\vec{\eta}$, which exists in a $10$ dimensional space, in the two dimensional
plane $r,s$. Some parameters are linked together though. The mathematical
analysis in the supplementary section 6.3.4 shows us that the skeletons
obtained for a fixed value of $\tau_{B},w^{+}$, can be extrapolated to other
values ${\tau^{\prime}}_{B},{w^{\prime}}^{+}$ by the simple rescaling
$r^{\prime}=r,s^{\prime}=\left(\,\frac{{w^{\prime}}^{+}\,{\tau^{\prime}}_{B}}{w^{+}\,\tau_{B}}\,\right)^{2}\,s$.
Using this rescaling, the map of Figure 5 has been drawn for a specific value
of $w^{+}=8.5$ Hz and $\tau_{B}=30$ ms. This corresponds to mean values of
these parameters, averaged over the set of experimental cells (in CTL
conditions). In this two dimensional representation, extra information coming
from the other parameters $\tau_{G},w^{B}_{G},w^{A}_{G},A_{0},b_{0}$ is lost.

Figure 5, left, shows us the repartition of cells in the map, in CTL
conditions. A few of them are monophasic, but many of them are biphasic, with
a significant proportion close to the polyphasic region and showing residual
oscillations. The figure 5, right, shows the same cells in CNO conditions. The
main observation is that CNO (right panel), does not dramatically change the
repartition of cells in the RFs map. This is made more explicit in the bottom
panels of figure 5. We show the mean and standard deviation of the main model
parameters: $\tau_{A},\tau_{B},w^{-},w^{+},b_{0}$ in CTL and CNO conditions,
separating the two subclasses of investigated genes: Grik4 and Scnn1a, and
separating ON or OFF cells. These parameters are essentially constant showing
that there is no statistical trend induced by CNO. This is in agreement with
our previous paper (Hilgen et al., , 2022) where we showed that Scnn1a or
Grik4 groups actually include multiple anatomical cell types, which suggest
that our model or fitting does not introduce any obvious bias.

Figure 5: Top. Repartition of fitted experimental cells in the RFs map. Each
point corresponds to an experimental cell. Top left. CTL conditions. Top
right. CNO conditions Middle. Mean and standard deviation of $\tau_{A}$
(left), $\tau_{B}$ (center), $r$ (right), fitted from experiments, for genes
Grik4 and Scnn1a, in CTL and CNO conditions. We have separated the estimation
for OFF cells (blue), ON Cells (green) and all cells (red). Bottom. Mean and
standard deviation of $w^{-}$ (left), $w^{+}$ (center), $s$ (right), fitted
from experiments. The representation is the same as the previous row.

The situation is radically different when investigating the effect on
individual cells. Indeed, the application of CNO makes some cells to move
their representative point from one region in the RF map to the other, thereby
drastically changing the cell’s response. Two examples are shown in Figure 6
and 7. In Figure 6, the application of CNO induces a motion of the
representative point in the RFs map. However, as the cell is close to the area
separating the monophasic from the biphasic phase, this motion impacts
dramatically the shape of the time response. Note that this motion looks tiny
because it is plotted in log scale, but in reality reflects a relatively large
change in the RGC properties. This corresponds therefore to a relatively big
change in the properties of the ganglion cell. Actually, the RFs map can be
refined by plotting the value of the main period $T_{1}$ (corresponding to the
main peak in the power spectrum) as shown in the top right figure. The dashed
black lines correspond to the separation between the monophasic, biphasic, and
polyphasic regions. Note that the switch from bi- to poly-phasic corresponds
to additional peaks in the spectrum that are not visible when considering the
period of the first peak, $T_{1}$. (This transition is seen in Fig. 7 where
the period of the second peak is shown). In CTL conditions, the cell is
located in the green region with a high $T_{1}$ of order $600$ ms, in the
limit of experimental resolution. With CNO, the cell switches to the grey
region where $T_{1}$ is of order $300$ ms. At the bottom of the figures, the
synthesis panel (same representation as in Figure 3) for CTL (left) and CNO
conditions (right) is presented.

Figure 6: CNO may change the cell response. The upper left panels shows the
representative point of a cell, labelled $83_{2}$, in the RF map, in CTL
(purple point) and CNO (red point) conditions. This shows the displacement of
the representative point of the cell when CNO is applied. The upper right
panel shows the same motion, in the same space $(r,s)$, but considering here
the main peak in the power spectrum, corresponding to a characteristic period
$T_{1}$. The dashed black lines correspond to the separation between the
monophasic, biphasic, and polyphasic regions.

One switches from a region with a high $T_{1}$ of order $600$ ms (green area)
in CTL, to a region where $T_{1}$ is of order $300$ ms (grey area) in CNO
conditions. As shown in the bottom panel this displacement corresponds to a
switch from monophasic to biphasic region.

In Figure 7 we show a motion inducing a switch from polyphasic to biphasic.
Here, CNO drives the cell from the boundary of the polyphasic region to the
biphasic region. The representation is the same as for Figure 6 except that
the top right figure displays the period of the second period $T_{2}$
(secondary peak in the power spectrum).

Figure 7: CNO may change the cell response. Left top. Application of CNO moves
the representative point of the cell (here label $52_{1}$) in from polyphasic
to biphasic region resulting in a change in the cell response (Bottom panel).
The figure Right top shows how the period of polyphasic oscillations depends
on the place in the RFs map.

Remark. The range of parameters r, s is very broad, and seems to be stretched
in terms of biological plausibility. In Fig. 5, $s$ can raise up to $s=200$
which means that ACs inhibitory feedback is $200$ times stronger than BCs
excitation onto ACs which looks quite implausible. In the same vein, we
observe values of $r=30$ corresponding to ACs having a $30$-fold longer time
constant than BCs, which might be questionable. These values must actually not
be taken literally as they correspond to several artifacts. First, as
explained above, $s$ has been rescaled to match a RF map drawn from a set of
reference parameters with e.g. $w^{+}=8.5$ Hz, and $\tau_{B}=30$ ms. In
addition, as mentioned earlier, these synaptic weights are actually effective
interactions that mimic the effect of cells populations. These extreme values
of $(r,s)$ are therefore somewhat reflecting the fact that we are dealing with
interactions among multiple cell types in the inner retina.

##  5 Discussion

In this paper, we investigated the role of AC-mediated lateral connectivity in
the response of RGCs to visual stimuli. Our conjecture was that these
responses are strongly constrained by such lateral connectivity. Based on the
mathematical analysis of a network featuring the interaction of BCs-ACs-RGCs
we were able to produce an analytic form for the spatio-temporal response
(receptive fields) of all cell types in the model. This finding has
significant implications for the usefulness, identifiability (i.e. its
parameters could be obtained from experimental data) and interpretability of
the model. First, it provides an algorithmic way to fit the model parameters
to the light responses recorded from mouse RGCs, using the analytical formula
of the RF. This means that we are able not only to find the parameters that
best fit the variables concerning the RGCs responses, but also to infer the
possible behaviour of ACs and BCs leading to the RGCs responses, even if we
don’t measure them experimentally. Second, it provides an intuitive
understanding of the role of various model variables and highlights the impact
of two phenomenological parameters (with a physical meaning) on the spatio-
temporal response (i), the intensity of the interactions BCs-ACs, and, (ii),
the characteristic time scale of these cells responses. This can be summarized
in the two dimensional RF maps, where one observes phases corresponding to
different modalities in the response. We were able to validate experimentally
these modelling results, based on the ability to pharmacologically modify the
level of ACs and RGCs neural activity using pharmacogenetics (DREADD-CNO). We
would like now to comment some caveats and potential extensions of this work.

##### Polyphasic phase.

Although we observe about $40\%$ of polyphasic cells in the experimental plots
(characterized by secondary peaks in the power spectrum) the model has
difficulty to properly fit them. This is visible in figure 5 where no cell is
within the polyphasic region whereas the secondary peaks are clearly visible
in the experimental power spectra (see Figure 7). This can be explained by
several factors. First, the fitting method, trying to estimate a model with
$10$ parameters from a trace with a few hundreds of points, has clearly
limits. In particular, the secondary peaks, having a few points in the power
spectrum, are hard to capture and require a patient fine tuning. The main
limitations may also come from the model itself, as developed in the
following.

##### Inhomogeneities.

Our model is quite homogeneous. It assumes that there is only one subclass for
each cell type, represented by a unique set of parameters, the connectivity is
fairly regular … This homogeneity was the key to answer a specific question:
understanding better the potential role of the BCs-ACs network on the RGCs
response, both from a theoretical and experimental perspective. This study was
actually done in the spirit of theoretical physics: find a set of canonical
equations, grounded on reality (here biophysics), and mathematically analyse
their behaviour with a minimal set of reduced (dimensionless) parameters.
Then, confront the model predictions to experiments and propose a simple
representation of the observed effects. This type of approach is useful,
because it allows one to understand the action of cells within a large retinal
network at a level of detail and comprehensiveness generally not possible by
numerical simulation alone. Yet, one may wonder what would happen if more
parameters were taken into consideration, to get closer to real retinal
networks. For example, what would be the consequences of introducing more
sustained and transient bipolar and amacrine cells ? Or considering more
general forms of connectivity, e.g. connections between amacrine cells ?

As mentioned at the end of section 3.1.1, equation (3) extends to more cell
types and general connectivity. The same holds for equations (13) and (15).
For example, considering $M_{B}$ types of BCs, each with a specific RF
${\mathcal{K}}_{m}(x,y,t)={\mathcal{K}}_{{B}_{S_{m}}}(x,y)\,{\mathcal{K}}_{{B}_{T_{m}}}(t)$,
equation (15) would extend to:

${\mathcal{K}}_{{}_{\alpha}}(x,y,t)=\sum_{\beta=1}^{3N}\left(\,{\mathcal{P}}_{\alpha\beta}\,\sum_{m=1}^{M_{B}}U_{\beta
m}(t)\,\times\,\sum_{\gamma\in
m}{\mathcal{P}}^{-1}_{\beta\gamma}\,\varpi_{\beta\gamma}\,{\mathcal{K}}_{{B}_{S_{\gamma}}}(x,y)\,\right),$
(21)

where $U_{\beta m}(t)\equiv\left[\,e_{\beta}\stackrel{{\scriptstyle
t}}{{\ast}}{\mathcal{K}}_{{B}_{T_{m}}}\,\right](t)$ and where $\gamma\in m$
stands for "the BC $\gamma$ belongs to population m". Beyond that, we are, at
the moment not able to tell what would be the impact of such heterogeneity on
the model. The reason is fairly simple. Adding heterogeneity will increase the
number of parameters, and, in general, will prevent us to go as far as we did
in this paper, constructed in the spirit of having a mathematical control of
the results. On technical grounds, the structure of eigenvalues and
eigenvectors is expected to dramatically change in a way that we cannot, in
general, determine analytically, in contrast to the present work, and it would
require numerical explorations in a huge parameters space.

Adding heterogeneity clearly necessitates to precisely define which question
we want to answer, and, from this point of view, specific forms of
heterogeneity could be studied in an extension of this work. A simple idea
would be to introduce some randomness in the parameters characterizing BCs and
ACs. For example, considering random leak times in order to have a
distribution of characteristic times $\tau_{A},\tau_{B}$ to study e.g. how we
depart from the result obtained here when increasing the variance of these
times. Or, considering random versions of the incidence matrices
$\Gamma^{{A}}_{{B}},\Gamma^{{B}}_{{A}}$. An example of this has been
considered in the paper Souihel and Cessac,  (2021) to analyse anticipation
effects. In any case, such questions would require long investigations beyond
the scope of this paper.

##### Rectification.

The model includes weak non linearities (rectification) that were neglected in
the mathematical computations. The effects of such rectification can be
mathematically investigated (Cessac, , 2022). Mainly, rectification projects
dynamics on the subspace of non rectified cells. This means that the
dimensionality of the dynamical system changes in time, depending on the
stimulus and network parameters, with strong consequences on the spectrum of
${\mathcal{L}}$, and thereby, on the power spectrum briefly discussed in the
supplementary section 6.3.4.

##### Local non linearities vs network effects.

Additional non linearities take place in retinal dynamics. Ion channels have a
non linear behaviour inducing phenomena such as bifurcation and bursting,
essential, for example, in the development of the retina were bursting
Starburst Amacrine Cells generate retinal waves (Hennig et al., , 2009;
Karvouniari et al., , 2019; Cessac and Matzakou-Karvouniari, , 2022). In
addition, gain control plays also a central role in the response to spatio-
temporal stimuli inducing, for example, retinal anticipation (Berry et al., ,
1999; Chen et al., 2013a, ; Souihel and Cessac, , 2021).

Although our study was limited to responses to full-field flashes we would
like to extend the consequence of our analysis to more complex stimuli. First,
the presence of peaks in the power spectrum implies the existence of
resonances, that is preferred frequencies for the RGCs. Exciting a cell with a
resonance frequency will produce a maximal response. When applying a stimulus
like the Chirp stimulus (Hilgen et al., , 2022; Baden et al., , 2016) there is
a phase where periodic flashes, with constant contrast but increasing
frequency, are applied. One observes frequently a bump in the experimental
RGCs response that might correspond to such a resonance.

In addition to preferred time frequencies, our analysis also shows that the
response of RGCs, induced by the network, may also involve specific space
scales. Mathematically, these scales appear in the eigenvectors of the
transport operator (see eq. (49) in the Supplementary section 6.3.5). The
practical implication would be that, presenting a localised, time periodic,
stimulus at a resonant frequency and with small radius, and increasing slowly
this radius, one may observe scales where the response is maximal. One of
these scales may correspond to the size of the RF but we conjecture that there
should be other, larger, scales where this phenomenon appears. This would
actually be a way to disentangle intrinsic responses of cells, from network
induced responses, by blocking the ACs synapses (e.g. strychnine for
glycinergic cells). More generally, the existence of time resonances and
preferred space scales would also induce resonant response to moving objects
with the appropriate speed. Such effects could also be related to mechanisms
giving rise to anticipatory waves (Menz et al., , 2020).

To conclude, this research provides mathematical insights to explore the
potential role of the amacrine cells network in vision processing. First, it
brings in the field of retinal modelling methods and concepts from dynamical
systems theory, especially, the idea that concerted population activities can
be understood with a few reduced parameters, still reproducing experimental
observations. Second, it makes a step forward to formalizing the concept that
the response of RGCs to stimuli is not only the result of intrinsic cells
characteristic (e.g. their morphology), but it also depends of their
interaction with other cells, and on the stimulus itself.

This paper was a first step to confront this approach to experiments, as a
proof of concept opening up several research tracks. First, it could be used
to infer the potential structure of spike correlations generated by a RGCs
population when responding to a moving object. While it is easy to show, in
this model, that lateral inhibition decorrelates the response to full field
stimuli, the situation is dramatically different when an object is moving in
its visual field. The resulting wave of activity generates non trivial,
transient, spatio-temporal correlations which may contain fundamental
information deciphered by the cortex. The corresponding non stationary spike
train distributions can be mathematically (Cessac, , 2011; Cessac et al., ,
2009, 2021) and numerically (Cessac et al., , 2017) studied. Pushing forward
this formalism could therefore be a key toward understanding better how the
retina encodes the visual information in a visual world made of motion.

Second, with the necessity of improving their biological plausibility (see
discussion above), this class of model could be used to fit online
experimental retinal responses to a simulator, by adapting the stimulus to the
model prediction, in a closed loop process (Benda et al., , 2007). This is in
the spirit of current trends of research attempting to construct data driven
models of the retina (Schröder et al., , 2020), although our approach targets
to reduce as much as possible the number of parameters and maximize the
mathematical control so as to avoid the "black box" effect.

Finally, our work could be used in future studies to explore the role of other
RGC subclasses or other retinal neurons and their interactions. In addition,
it could be used to disassemble the components of other retinal circuits, by
manipulating the activity of specific neurons. It could also potentially
benefit research in other parts of the nervous system, as fundamental
properties of the inner retina are shared with other parts of the brain.

## Acknowledgments

This project was funded by the Leverhulme Trust (RPG-2016-315 to ES and BC),
by Newcastle University (Faculty of Medical Sciences Graduate School).

## Supplementary

##  6 Derivation of the model equations and properties

###  6.1 Generic equations for voltages

We start from the fundamental equation of neuronal dynamics (charge
conservation) and make several simplifying assumptions.

First, we assume that neurons are reduced to points. They all have the same
leak conductance, $g_{L}$, the same leak reversal potential,
${\mathcal{E}}_{L}$ and the same membrane capacity, $C$. They are passive (no
active ion channels). This last assumption is based on the fact that most of
the retinal neurons considered here (BCs and ACs) are not spiking. From these
assumptions we get rid of equations for activation-inactivation variables
present in most spiking models (Ermentrout and Terman, , 2010). For RGCs,
spiking is ruled by a LN type model (see main text).

As exposed in the main text we model the CNO effect by a current of the form
$-g_{CNO_{T}}\,\left(\,V_{q}-{\mathcal{E}}_{CNO_{T}}\,\right)$ where $T$ is
the cell type (e.g. ACs or RGCs) so that the CNO conductance is constant and
depends only on the cell type. ${\mathcal{E}}_{CNO_{T}}$ is the corresponding
reversal potential. We insist, however, that this is only a modelling short-
cut allowing to mathematically treat inhibitory and excitatory effects on an
equal footing, with only one parameter, the CNO conductance $g_{CNO_{T}}$.

The membrane potential $V_{q}$ of neuron $q$ then obeys the equation:

$C\frac{dV_{q}}{dt}=-g_{L}\left(\,V_{q}-{\mathcal{E}}_{L}\,\right)-g_{CNO_{T}}\left(\,V_{q}-{\mathcal{E}}_{CNO_{T}}\,\right)-\sum_{Y,syn}\sum_{p\in
Y}g^{(Y,p)}_{q}\,\left(\,V_{q}-{\mathcal{E}}_{Y}\,\right)\,+\,i_{q}(t),$ (22)

where $\sum_{Y,syn}$ is the sum over all possible types of synaptic currents
shaping the membrane potential of neuron $q$ (for example the current
resulting from glutamate transmitter with NMDA receptors). The sum $\sum_{p\in
Y}$ stands for the contribution of presynaptic neurons $p$ connecting to post-
synaptic neuron $q$ via a synapse of type $Y$. The term $i_{q}$ is an external
current (here, the photoreceptors current that input BCs). The synaptic
conductance $g^{(Y,p)}_{q}$ depends, in general, on the pre-synaptic voltage
$V_{p}$ and additional activation or inactivation variables. Here, we assume
that it depends only on the pre-synaptic voltage and take the form
$g^{(Y,p)}_{q}=\lambda^{(Y,p)}_{q}{\mathcal{N}}^{(Y)}(V_{p})$ where
$\lambda^{(Y,p)}_{q}$ is a constant, positive factor, and
${\mathcal{N}}^{(Y)}$ is the piecewise linear function introduced in (4).

We shift the voltages by ${\mathcal{E}}_{L}$ (i.e. replacing $V_{q}$ by
$V_{q}-{\mathcal{E}}_{L}$, ${\mathcal{E}}_{X}$ by
${\mathcal{E}}_{X}-{\mathcal{E}}_{L}$ and so on). We introduce the
characteristic time scale of the membrane response for neuron $q$:

$\tau_{q}=\frac{C}{g_{L}\,+\,g_{CNO_{T}}\,+\,\sum_{Y,syn}\sum_{p\in
Y}\lambda^{(Y,p)}_{q}{\mathcal{N}}^{(Y)}(V_{p})}$

$\tau_{q}$ depends on the pre-synaptic neurons states (via the conductances
$g^{(Y,p)}_{q}$). We can write the time dependence of voltages as
$V_{p}(t)=V_{p}^{\ast}+\delta V_{p}(t)$ where $V_{p}^{\ast}$ is the rest
state, i.e. the voltage reached by neurons, in the absence of stimulation,
given by:

$V_{q}^{\ast}=\frac{g_{CNO_{T}}\,{\mathcal{E}}_{CNO_{T}}\,+\,\sum_{Y}\sum_{p\in
Y}\lambda^{(Y,p)}_{q}\,{\mathcal{N}}^{(Y)}(V_{p}^{\ast})\,{\mathcal{E}}_{Y}}{g_{L}\,+\,g_{CNO_{T}}\,+\,\sum_{Y,syn}\sum_{p\in
Y}\lambda^{(Y,p)}_{q}{\mathcal{N}}^{(Y)}(V_{p}^{\ast})}.$ (23)

This is a complex, non linear, system of equations where the rest voltages of
all neurons are dependent.

The term $\delta V_{p}(t)$ corresponds to the time variations of the voltage.
What makes the dynamics complex are precisely these time variations inducing
important non linear effects. Here, however, we are investigating effects on
the retinal response when applying CNO, considering that CNO acts on a time
scale quite larger than the characteristic time $\tau_{q}$. That is, we
neglect the fluctuations $\delta V_{p}(t)$ and assume that the characteristic
time $\tau_{q}$ is only controlled by the rest part of the voltages
$V_{p}^{\ast}$, depending on the CNO conductance. Thus, we can replace
${\mathcal{N}}^{(Y)}(V_{p})$ by ${\mathcal{N}}^{(Y)}(V_{p}^{\ast})$ in the
expression of $\tau_{q}$:

$\tau_{q}=\frac{\tau_{L}}{1+\tau_{L}\left(\,\frac{\zeta_{T,q}}{{\mathcal{E}}_{CNO_{T}}}+\sum_{Y,syn}\frac{1}{{\mathcal{E}}_{Y}}\sum_{p\in
Y}W^{(Y,p)}_{q}{\mathcal{N}}^{Y}(V_{p}^{\ast})\,\right)},$ (24)

where $\tau_{L}=\frac{C}{g_{L}}$ is the leak characteristic time.

Finally, we introduce the synaptic weights:

$W^{(Y,p)}_{q}=\frac{\lambda^{(Y,p)}_{q}\,{\mathcal{E}}_{Y}}{C},$ (25)

the CNO polarization:

$\zeta_{T,q}=\frac{g_{CNO_{T}}\,{\mathcal{E}}_{CNO_{T}}}{C},$ (26)

and the forcing term $F_{q}(t)=\,\frac{i_{q}(t)}{C}$, and we write:

$\frac{dV_{q}}{dt}=-\frac{1}{\tau_{q}}V_{q}\,+\,\sum_{Y}\sum_{p\in
Y}W^{(Y,p)}_{q}\,{\mathcal{N}}^{(Y)}(V_{p})\,+\,\zeta_{T,q}\,+\,F_{q}(t),\quad
q=1\dots N,$ (27)

where we dropped $"syn"$ in $\sum_{Y,syn}$ for legibility. This is the general
form of cells voltage equations used in the paper, leading to eq. (3).

###  6.2 Rest states, characteristic times and CNO dependence

We study now the dependence of the rest states and characteristic times in
CNO, for the model (3). The details of the computations can be found on the
web page https://team.inria.fr/biovision/mathematicalderivationscno/. From
(24) the characteristic times take the form:

$\left\\{\begin{array}[]{lll}\tau_{B_{i}}&=&\frac{\tau_{L}}{1+\frac{\tau_{L}}{{\mathcal{E}}_{A}}\sum_{j=1}^{N}W^{{A}_{j}}_{{B}_{i}}{\mathcal{N}}^{A}(V_{A_{j}}^{\ast})},\\\
\tau_{A_{j}}&=&\frac{\tau_{L}}{1+\tau_{L}\left(\,\frac{\zeta_{A}}{{\mathcal{E}}_{CNO_{A}}}+\frac{1}{{\mathcal{E}}_{B}}\sum_{i=1}^{N}W^{{B}_{i}}_{{A}_{j}}{\mathcal{N}}^{B}(V_{B_{i}}^{\ast})\,\right)},\\\
\tau_{G_{k}}&=&\frac{\tau_{L}}{1+\tau_{L}\left(\,\frac{\zeta_{G}}{{\mathcal{E}}_{CNO_{G}}}+\frac{1}{{\mathcal{E}}_{B}}\sum_{i=1}^{N}W^{{B}_{i}}_{{G}_{k}}{\mathcal{N}}^{B}(V_{B_{i}}^{\ast})+\frac{1}{{\mathcal{E}}_{A}}\sum_{j=1}^{N}W^{{A}_{j}}_{{G}_{k}}{\mathcal{N}}^{A}(V_{A_{j}}^{\ast})\,\right)}.\end{array}\right.$
(28)

Here, we see an evident and expected dependence on CNO. Increasing $\zeta_{A}$
decreases the characteristic times of the membrane. This is, however, not
complete, as the rest states depend themselves on CNO. Indeed, assume, for
simplicity, that the system is uniform in space, so that the characteristic
times do not depend on the cell index (only on the cell type). Then, the rest
states are given by (in vector form):

$\left\\{\begin{array}[]{lll}\vec{V}^{\ast}_{B}&=&\tau_{B}\,\left(\,{\mathcal{I}}_{N}\,-\,\tau_{A}\tau_{B}W^{{A}}_{{B}}W^{{B}}_{{A}}\,\right)^{-1}.\left(\,\tau_{A}W^{{A}}_{{B}}\zeta_{A}\,-\,\tau_{A}\theta_{B}\,W^{{A}}_{{B}}W^{{B}}_{{A}}\,-\,\theta_{A}\,W^{{A}}_{{B}}\,\right)\vec{1}_{N},\\\
\vec{V}^{\ast}_{A}&=&\tau_{A}\,\left(\,{\mathcal{I}}_{N}\,-\,\tau_{A}\tau_{B}W^{{B}}_{{A}}W^{{A}}_{{B}}\,\right)^{-1}.\left(\,\zeta_{A}\,-\,\tau_{B}\theta_{A}W^{{B}}_{{A}}W^{{A}}_{{B}}\,-\,\theta_{B}\,W^{{B}}_{{A}}\,\right)\vec{1}_{N},\\\
\vec{V}^{\ast}_{G}&=&\tau_{G}W^{{A}}_{{G}}\vec{V}^{\ast}_{A}\,+\,\tau_{G}W^{{B}}_{{G}}\,\vec{V}^{\ast}_{B}\,+\,\tau_{G}\left(\,\zeta_{G}\,-\,\theta_{A}W^{{A}}_{{G}}\,-\,\theta_{B}W^{{B}}_{{G}}\,\right).\vec{1}_{N},\end{array}.\right.$
(29)

where ${\mathcal{I}}_{N}$ is the $N\times N$ identity matrix and $\vec{1}_{N}$
is the unit vector in $N$ dimension. Here, we have assumed that the rest
states are above threshold and that
${\mathcal{I}}\,-\,\tau_{A}\tau_{B}W^{{A}}_{{B}}W^{{B}}_{{A}}$ and
${\mathcal{I}}\,-\,\tau_{A}\tau_{B}W^{{B}}_{{A}}W^{{A}}_{{B}}$ are invertible.
See the web page https://team.inria.fr/biovision/mathematicalderivationscno/
for the general case.

The system of equations (28) and (29) are thus non linearly entangled.
Especially, there is a nice feedback effect appearing in feedback "loop" terms
like $W^{{A}}_{{B}}W^{{B}}_{{A}}$. Actually, e.g. the term
$\left(\,{\mathcal{I}}_{N}\,-\,\tau_{A}\tau_{B}W^{{A}}_{{B}}W^{{B}}_{{A}}\,\right)^{-1}$
involves loops $\left(\,W^{{A}}_{{B}}W^{{B}}_{{A}}\,\right)^{n}$ of order
$n>0$, with an amplitude decaying with $n$. Consider now the role of the
parameter $\zeta_{A}$. It impacts directly ACs and indirectly, through the
network, the BCs and RGCs activity. One observes, first, a direct effect on
the ACs rest state. Looking at the numerator of (29), increasing $\zeta_{A}$
would have the effect of depolarizing the cell (if
${\mathcal{E}}_{CNO_{A}}>0$) or hyperpolarizing it (if
${\mathcal{E}}_{CNO_{A}}<0$), if $\tau_{A}$ and $\tau_{B}$ were independent of
CNO. However, CNO acts as well on these times. The net effect of $\zeta_{A}$
is thus non linear and cannot be anticipated by simple arguments.

In the case of section 3.1.4 with a reduced set of parameters the equations
for characteristic times reduce to (8) while the rest states are given by (9).

###  6.3 Mathematical analysis of network dynamics

In the following, the notation
$\mathop{\mathrm{diag}}\left(\,x_{n}\,\right)_{n=1}^{N}$ denotes a diagonal
$N\times N$ matrix with diagonal entries $x_{n}$.

####  6.3.1 Joint dynamics

The joint dynamics of all cells type is given by the dynamical system (3). We
use Greek indices $\alpha,\beta,\gamma=1\dots 3N$ and define the state vector
$\vec{{\mathcal{X}}}$ with entries:

${\mathcal{X}}_{\alpha}=\left\\{\begin{array}[]{llll}&V_{{B}_{i}},&\alpha=i,&i=1\dots
N;\\\ &V_{{A}_{j}},&\alpha=N+j,&j=1\dots N;\\\
&V_{{G}_{k}},&\alpha=2N+k,&k=1\dots N.\end{array}\right.$

We introduce $\vec{{\mathcal{F}}}$ with entries:

${\mathcal{F}}_{\alpha}=\left\\{\begin{array}[]{llll}&F_{{B}_{i}},&\alpha=i,&i=1\dots
N;\\\ &\zeta_{A},&\alpha=N+j,&j=1\dots N;\\\ &\zeta_{G},&\alpha=2N+k,&k=1\dots
N;\end{array}\right.$

and the rectification vector $\vec{{\mathcal{R}}}(\vec{{\mathcal{X}}})$ with
entries:

${\mathcal{R}}_{\alpha}(\vec{{\mathcal{X}}})=\left\\{\begin{array}[]{llll}&{\mathcal{N}}^{(B)}\left(\,V_{{B}_{i}}\,\right),&\alpha=i,&i=1\dots
N;\\\ &{\mathcal{N}}^{(A)}\left(\,V_{{A}_{j}}\,\right),&\alpha=N+j,&j=1\dots
N;\\\ &0,&\alpha=2N+k,&k=1\dots N;\end{array}\right.$

We introduce the $3N\times 3N$ matrices:

${\mathcal{T}}=\left(\,\begin{array}[]{cccccc}&\mathop{\mathrm{diag}}\left[\,\tau_{B_{i}}\,\right]_{i=1\dots
N}&&0_{NN}&&0_{NN}\\\
&0_{NN}&&\mathop{\mathrm{diag}}\left[\,\tau_{A_{j}}\,\right]_{j=1\dots
N}&&0_{NN}\\\
&0_{NN}&&0_{NN}&&\mathop{\mathrm{diag}}\left[\,\tau_{G_{k}}\,\right]_{k=1\dots
N}\end{array}\,\right),$ (30)

characterizing the characteristic integration times of cells,

${\mathcal{W}}=\left(\,\begin{array}[]{cccccc}&0_{NN}&&W^{{A}}_{{B}}&&0_{NN}\\\
&W^{{B}}_{{A}}&&0_{NN}&&0_{NN}\\\
&W^{{B}}_{{G}}&&W^{{A}}_{{G}}&&0_{NN}\end{array}\,\right),$ (31)

summarizing chemical synapses interactions. Then, the dynamical system (3)
reads, in vector form:

$\frac{d\vec{{\mathcal{X}}}}{dt}\,=\,-\,{\mathcal{T}}^{-1}.\vec{{\mathcal{X}}}+{\mathcal{W}}.\vec{{\mathcal{R}}}(\vec{{\mathcal{X}}})+\vec{{\mathcal{F}}}(t).$
(32)

We remark that eq. (32) has a specific product structure: the dynamics of RGCs
is driven by BCs and ACs with no feedback. This means that one can study first
the coupled dynamics of BCs and ACs and then the effect on RGCs.

####  6.3.2 Linear evolution

We consider the evolution of eq. (32) from an initial time $t_{0}$. Typically,
$t_{0}$ is a reference time where the network is at rest, before the stimulus
is applied. The dynamical system has almost the form of a non-autonomous
linear system driven by the term $\vec{{\mathcal{F}}}(t)$. There is however a
weak non linearity, due to the piecewise linear rectification appearing in the
term $\vec{{\mathcal{R}}}(\vec{{\mathcal{X}}})$. Therefore, when the voltages
of all cells are large enough the system is linear. Mathematically, there is a
domain of $\mathbbm{R}^{3N}$:

${\mathcal{D}}=\left\\{\,V_{{B}_{i}}\geq\theta_{B},V_{{A}_{j}}\geq\theta_{A},i,j=1\dots
N\,\right\\},$ (33)

where $R\left(\,\vec{{\mathcal{X}}}\,\right)$ is linear so that eq. (32) is
linear too (check (Cessac, , 2022) for more details).

From now on we consider this linear case. We write
${\mathcal{L}}\,=\,-\,{\mathcal{T}}^{-1}\,+\,{\mathcal{W}}$ so that:

${\mathcal{L}}=\left(\,\begin{array}[]{cccccc}&-\mathop{\mathrm{diag}}\left[\,\frac{1}{\tau_{B_{i}}}\,\right]_{i=1\dots
N}&&W^{{A}}_{{B}}&&0_{NN}\\\
&W^{{B}}_{{A}}&&-\mathop{\mathrm{diag}}\left[\,\frac{1}{\tau_{A_{j}}}\,\right]_{j=1\dots
N}&&0_{NN}\\\
&W^{{B}}_{{G}}&&W^{{A}}_{{G}}&&-\mathop{\mathrm{diag}}\left[\,\frac{1}{\tau_{G_{k}}}\,\right]_{k=1\dots
N}\end{array}\,\right),$ (34)

We introduce the $N$ dimensional vector $\vec{1}_{N}=\vec{(}1)_{i=1}^{N}$, and
the $3\,N$ dimensional vector
$\vec{{\mathcal{C}}}=\left(\,\begin{array}[]{ccc}-\theta_{A}\,W^{{A}}_{{B}}.\vec{1}_{N}\\\
-\theta_{B}\,W^{{B}}_{{A}}.\vec{1}_{N}\\\
-\left(\,\theta_{B}\,W^{{B}}_{{G}}.\vec{1}_{N}\,+\,\theta_{A}\,W^{{A}}_{{G}}\,\right).\vec{1}_{N}\end{array}\,\right)$
and (32) reads
$\frac{d\vec{{\mathcal{X}}}}{dt}={\mathcal{L}}.\vec{{\mathcal{X}}}\,+\,\vec{{\mathcal{F}}}(t)\,+\,\vec{{\mathcal{C}}}$.

We assume that ${\mathcal{L}}$ is invertible. This assumption, and more
generally, the spectrum of ${\mathcal{L}}$ is further discussed in section
6.3.4. The general solution of eq. (32) is:

$\vec{{\mathcal{X}}}(t)=e^{{\mathcal{L}}(t-t_{0})}.\vec{{\mathcal{X}}}(t_{0})\,+\,\int_{t_{0}}^{t}e^{{\mathcal{L}}(t-s)}\,\vec{{\mathcal{F}}}(s)\,ds\,-\,{\mathcal{L}}^{-1}.\left[\,{\mathcal{I}}_{3N,3N}\,-\,e^{{\mathcal{L}}(t-t_{0})}\,\right].\vec{{\mathcal{C}}}.$
(35)

where ${\mathcal{I}}_{3N,3N}$ is the $3N$ dimensional identity matrix.

Although this equation is general, it actually stands when one can define a
notion of asymptotic regime. That is, when ${\mathcal{L}}$ has stable
eigenvalues (eigenvalues with a strictly negative real part). The spectrum of
${\mathcal{L}}$ is studied below and conditions ensuring the stability of
eigenvalues are given. Here, we are going to assume that eigenvalues are all
stable and that $t-t_{0}$ is large so that we can remove the transient term
$e^{{\mathcal{L}}(t-t_{0})}.\vec{{\mathcal{X}}}(t_{0})$ depending on the
initial condition $\vec{{\mathcal{X}}}(t_{0})$. In addition, the last term
converges to:

$\vec{{\mathcal{X}}}^{\ast}=-{\mathcal{L}}^{-1}.\vec{{\mathcal{C}}},$ (36)

the rest state of the linear system, which vanishes whenever the thresholds
$\theta_{A},\theta_{B}$ are set to $0$.

####  6.3.3 Derivation of eq. (10)

We note the eigenvalues of ${\mathcal{L}}$, $\lambda_{\beta},\beta=1\dots 3N$
and its eigenvectors, ${\mathcal{P}}_{\beta}$ (the columns of the matrix
${\mathcal{P}}$ transforming ${\mathcal{L}}$ in diagonal form). We consider
first the case $\vec{{\mathcal{C}}}=\vec{0}$. We have then, from (35):

${\mathcal{X}}_{\alpha}(t)=\sum_{\beta=1}^{3N}\sum_{\gamma=1}^{3N}{\mathcal{P}}_{\alpha\beta}{\mathcal{P}}^{-1}_{\beta\gamma}\int_{-\infty}^{t}e^{\lambda_{\beta}(t-s)}\,{\mathcal{F}}_{\gamma}(s)\,ds,$

where ${\mathcal{F}}_{\gamma}=F_{{B}_{i}}$, $\gamma=i=1\dots N$ (BCs).

We recall that, from (5),
$F_{{B}_{\gamma}}(t)=\frac{{\mathcal{V}}^{(drive)}_{{\gamma}}}{\tau_{B_{\gamma}}}\,+\,\frac{d{\mathcal{V}}^{(drive)}_{{\gamma}}}{dt}$,
so that:

$\int_{-\infty}^{t}e^{\lambda_{\beta}(t-s)}\,{\mathcal{F}}_{\gamma}(s)\,ds={\mathcal{V}}^{(drive)}_{{\gamma}}(t)\,+\,\varpi_{\beta\gamma}\,\int_{-\infty}^{t}e^{\lambda_{\beta}(t-s)}\,{\mathcal{V}}^{(drive)}_{{\gamma}}(s)\,ds,\quad\gamma=1\dots
N,$

for B cells, with
$\varpi_{\beta\gamma}=\frac{1}{\tau_{B_{\gamma}}}+\lambda_{\beta}$ and using
${\mathcal{V}}^{(drive)}_{{\gamma}}(-\infty)=0$.

For $\gamma=N+1\dots 2N$, ${\mathcal{F}}_{\gamma}=\zeta_{A}$ (ACs) we have:

$\int_{-\infty}^{t}e^{\lambda_{\beta}(t-s)}\zeta_{A}\,ds=-\frac{\zeta_{A}}{\lambda_{\beta}}.$

Finally, for $\gamma=2N+1\dots 3N$, ${\mathcal{F}}_{\gamma}=\zeta_{G}$ (RGCs):

$\int_{-\infty}^{t}e^{\lambda_{\beta}(t-s)}\zeta_{G}\,ds=-\frac{\zeta_{G}}{\lambda_{\beta}}.$

We remark that:

$\begin{split}\sum_{\beta=1}^{3N}\sum_{\gamma=1}^{3N}{\mathcal{P}}_{\alpha\beta}{\mathcal{P}}^{-1}_{\beta\gamma}{\mathcal{V}}^{(drive)}_{{\gamma}}(t)&=\sum_{\gamma=1}^{3N}{\mathcal{V}}^{(drive)}_{{\gamma}}(t)\left(\,\sum_{\beta=1}^{3N}{\mathcal{P}}_{\alpha\beta}{\mathcal{P}}^{-1}_{\beta\gamma}\,\right)\\\
&=\sum_{\gamma=1}^{3N}{\mathcal{V}}^{(drive)}_{{\gamma}}(t)\delta_{\alpha\gamma}\\\
&={\mathcal{V}}^{(drive)}_{{\alpha}}(t).\end{split}$

It follows that:

${\mathcal{X}}_{\alpha}(t)={\mathcal{V}}^{(drive)}_{{\alpha}}(t)\,+{\mathcal{E}}^{(drive)}_{{\alpha}}(t)\,+\,{\mathcal{E}}^{(CNO_{A})}_{{\alpha}}\,+\,{\mathcal{E}}^{(CNO_{G})}_{{\alpha}},\quad\alpha=1\dots
3N.$

with:

${\mathcal{E}}^{(drive)}_{{\alpha}}(t)=\sum_{\beta=1}^{3N}\sum_{\gamma=1}^{N}{\mathcal{P}}_{\alpha\beta}\,{\mathcal{P}}^{-1}_{\beta\gamma}\varpi_{\beta\gamma}\,\int_{-\infty}^{t}e^{\lambda_{\beta}(t-s)}\,{\mathcal{V}}^{(drive)}_{{\gamma}}(s)\,ds,\quad\alpha=1\dots
3N,$
${\mathcal{E}}^{(CNO_{A})}_{{\alpha}}=\,-\,\zeta_{A}\,\sum_{\beta=1}^{3N}\sum_{\gamma=N+1}^{2N}\frac{{\mathcal{P}}_{\alpha\beta}\,{\mathcal{P}}^{-1}_{\beta\gamma}}{\lambda_{\beta}},\alpha=N+1\dots
2N;$
${\mathcal{E}}^{(CNO_{G})}_{{\alpha}}=\zeta_{G}\,\tau_{G},\quad\alpha=2N+1\dots
3N,$

where the last result comes from $\lambda_{\beta}\,=\,-\,\frac{1}{\tau_{G}}$,
for $\beta=2N+1\dots 3N$ (see next section). The expression of
${\mathcal{X}}_{\alpha}(t)$ corresponds to (10).

When ${\mathcal{C}}\neq\vec{0}$, there is an additional term corresponding to
the rest state (36).

####  6.3.4 Eigenvalues and eigenvectors of ${\mathcal{L}}$

##### Linear case.

We start from the eq. (34) of the linear operator ruling the dynamics in the
set ${\mathcal{D}}$ defined by (33). We consider, as in the main text, the
case where all characteristic times $\tau_{B_{i}}$ are equal to $\tau_{B}$,
all characteristic times $\tau_{A_{j}}$ are equal to $\tau_{A}$ and all
characteristic times $\tau_{G_{k}}$ are equal to $\tau_{G}$. Using the same
notations as the main text we have:

${\mathcal{L}}=\left(\,\begin{array}[]{cccccc}&-\frac{{\mathcal{I}}_{NN}}{\tau_{B}}&&-w^{-}\,\Gamma^{{A}}_{{B}}&&0_{NN}\\\
&w^{+}\,\Gamma^{{B}}_{{A}}&&-\frac{{\mathcal{I}}_{N}}{\tau_{A}}&&0_{NN}\\\
&w^{B}_{G}\,\Gamma^{{B}}_{{G}}&&w^{A}_{G}\,\Gamma^{{A}}_{{G}}&&-\frac{{\mathcal{I}}_{NN}}{\tau_{G}}\end{array},\,\right)$

where $0_{NN}$ is the $N\times N$ $0$ matrix and $I_{NN}$ the $N\times N$ $0$
identity matrix.

##### Eigenvalues and eigenvectors.

We consider the case where $\Gamma^{{A}}_{{B}}=\Gamma^{{B}}_{{A}}$. We note
$\kappa_{n},n=1\dots N$, the eigenvalues of $\Gamma^{{A}}_{{B}}$ ordered as
$\left|\,\kappa_{1}\,\right|\leq\left|\,\kappa_{2}\,\right|\leq\dots\leq\left|\,\kappa_{n}\,\right|$
and note the normalized eigenvectors $\vec{\psi}_{n}$, $n=1\dots N$.

We seek the eigenvalues, $\lambda_{\beta}$, and eigenvectors,
$\vec{{\mathcal{P}}}_{\beta}$, $\beta=1\dots 3N$, of ${\mathcal{L}}$. It is
evident, from the form of ${\mathcal{L}}$, that there are $N$ eigenvalues
$\lambda_{\beta}=-\frac{1}{\tau_{G}},\vec{{\mathcal{P}}}_{\beta}=\vec{e}_{\beta}$
where $\vec{e}_{\beta}$ is the canonical basis vector in direction $\beta$. We
attribute them the indices $\beta=2N+1\dots 3N$ as this indexing corresponds
to the form of ${\mathcal{L}}$ when $w^{-}=w^{+}=0$. We seek the $2N$
remaining eigenvalues-eigenvectors assuming that
$\vec{{\mathcal{P}}}_{\beta}$s has the form:

$\vec{{\mathcal{P}}}_{\beta}=\left(\,\begin{array}[]{ccc}\vec{\psi}_{n}\\\
\rho_{n}\vec{\psi}_{n}\\\ \vec{\phi}_{n}\end{array}\,\right),\quad n=1\dots
N,$ (37)

where $\rho_{n}$ is an unknown parameter and $\vec{\phi}_{n}$ a $N$
dimensional vector, to be determined, from the characteristic equation:

${\mathcal{L}}.\vec{{\mathcal{P}}}_{\beta}=\lambda_{\beta}\vec{{\mathcal{P}}}_{\beta}.$

This leads to the system of equations:

$\left\\{\begin{array}[]{llll}\lambda_{\beta}&=&-\frac{1}{\tau_{B}}-w^{-}\,\rho_{n}\,\kappa_{n};\\\
w^{+}\kappa_{n}-\frac{\rho_{n}}{\tau_{A}}&=&\rho_{n}\,\lambda_{\beta};\\\
\lambda_{\beta}\,\vec{\phi}_{n}&=&\left(\,w^{B}_{G}\,\Gamma^{{B}}_{{G}}+w^{A}_{G}\,\rho_{n}\,\Gamma^{{A}}_{{G}}\,\right)\,\vec{\psi}_{n}-\frac{1}{\tau_{G}}\,\vec{\phi}_{n},\end{array},\right.$
(38)

We first assume that $w^{-},w^{+}>0$ and later discuss the limit when these
quantities tend to zero. Combining the two first equations leads to a second-
order polynomial in the variable $\rho_{n}$,

$w^{-}\,\kappa_{n}\,\rho_{n}^{2}\,-\,\frac{1}{\tau}\,\rho_{n}\,+\,w^{+}\,\kappa_{n}=0,$
(39)

giving $2$ solutions for each $n$:

$\rho_{n}^{\pm}=\left\\{\begin{array}[]{lll}\frac{1}{2\,\tau\,w^{-}\,\kappa_{n}}\left(\,1\,\pm\,\sqrt{1-4\,\mu\,\kappa_{n}^{2}}\,\right),&\quad\kappa_{n}\neq
0,\frac{1}{\tau}\neq 0;\\\
w^{+}\,\tau\,\kappa_{n},&\quad\kappa_{n}=0,\frac{1}{\tau}\neq 0;\\\
\pm\sqrt{-\frac{w^{+}}{w^{-}}},&\quad\frac{1}{\tau}=0.\end{array}\right.$ (40)

where:

$\frac{1}{\tau}=\frac{1}{\tau_{A}}-\frac{1}{\tau_{B}}.$ (41)

and:

$\mu=w^{-}\,w^{+}\,\tau^{2}\geq 0.$ (42)

The $2N$ first eigenvalues of ${\mathcal{L}}$ are therefore given by:

$\lambda_{n}^{\pm}=\left\\{\begin{array}[]{llll}-\frac{1}{2\,\tau_{AB}}\mp\frac{1}{2\,\tau}\,\sqrt{1-4\,\mu\,\kappa_{n}^{2}},&\quad\frac{1}{\tau}\neq
0;\\\ &&\\\
-\frac{1}{\tau_{A}}\,\mp\,\sqrt{-w^{-}\,w^{+}\,\kappa_{n}^{2}},&\quad\frac{1}{\tau}=0.\end{array}\right.$
(43)

where:

$\frac{1}{\tau_{AB}}=\frac{1}{\tau_{A}}+\frac{1}{\tau_{B}}.$ (44)

We finally obtain $2N$ vectors $\vec{\phi}_{n}$:

$\vec{\phi}_{n}^{\pm}=\frac{1}{\lambda_{n}^{\pm}+\frac{1}{\tau_{G}}}\,\left(\,w^{B}_{G}\,\Gamma^{{B}}_{{G}}+w^{A}_{G}\,\rho_{n}^{\pm}\,\Gamma^{{A}}_{{G}}\,\right)\,\vec{\psi}_{n}.$
(45)

Let us now discuss the limit when $w^{-}$ or $w^{+}$ or both tend to $0$. If
$w^{-}=0$, $\rho_{n}=w^{+}\kappa_{n}\tau$ from (39). If $w^{+}=0$ there are
two solutions of $\eqref{eq:Polyrho}$, $\rho_{n}=0$ or $\rho_{n}=\frac{1}{\tau
w^{-}\kappa_{n}}$. Finally, when $w^{-}=w^{+}=0$, $\rho_{n}=0$ and the ansatz
(37) does not apply. Actually, in this case, ${\mathcal{L}}$ is diagonal, the
$N$ first eigenvalues are $-\frac{1}{\tau_{B}}$, the $N$ next eigenvalues are
$-\frac{1}{\tau_{A}}$. We have, in this case:
$\lambda_{n}^{+}=-\frac{1}{\tau_{B}}$ and
$\lambda_{n}^{-}=-\frac{1}{\tau_{A}}$. Therefore, we order eigenvalues and
eigenvectors of ${\mathcal{M}}$ such that the $N$ first eigenvalues are
$\lambda_{n}^{+},n=1\dots N$, and the $N$ next are $\lambda_{n}^{-},n=N+1\dots
2N.$

We finally end up with the following form for the eigenvalues and eigenvectors
of ${\mathcal{L}}$:

$\begin{array}[]{lllll}\lambda_{\beta}=\lambda_{n}^{+},\vec{{\mathcal{P}}}_{\beta}=\left(\,\begin{array}[]{ccc}\vec{\psi}_{n}\\\
\rho^{+}_{n}\vec{\psi}_{n}\\\
\frac{1}{\lambda_{n}^{+}+\frac{1}{\tau_{G}}}\,\left(\,w^{B}_{G}\,\Gamma^{{B}}_{{G}}+w^{A}_{G}\,\rho_{n}^{+}\,\Gamma^{{A}}_{{G}}\,\right)\,\vec{\psi}_{n}\end{array}\,\right),&\beta=n=1\dots
N,&\\\
\lambda_{\beta}=\lambda_{n}^{-},\vec{{\mathcal{P}}}_{\beta}=\left(\,\begin{array}[]{ccc}\vec{\psi}_{n}\\\
\rho^{-}_{n}\vec{\psi}_{n}\\\
\frac{1}{\lambda_{n}^{-}+\frac{1}{\tau_{G}}}\,\left(\,w^{B}_{G}\,\Gamma^{{B}}_{{G}}+w^{A}_{G}\,\rho_{n}^{-}\,\Gamma^{{A}}_{{G}}\,\right)\,\vec{\psi}_{n}\end{array}\,\right),&\beta=N+1\dots
2N,&n=1\dots N,\\\
\lambda_{\beta}=-\frac{1}{\tau_{G}},\vec{{\mathcal{P}}}_{\beta}=\vec{e}_{\beta},&\beta=2N+1\dots
3N.\end{array}$

(46)

##### Skeleton.

The eigenvalues $\lambda_{n}^{\pm}$ in 43 can be real or complex conjugated.
By increasing $\mu$, they become complex when:

$\mu>\frac{1}{4\,\kappa_{n}}\equiv\mu_{n,c}.$ (47)

In this case the real part is $-\frac{1}{2\,\tau_{AB}}$, the imaginary part is
$\pm\frac{1}{2\,\tau}\,\sqrt{1-4\,\mu\,\kappa_{n}^{2}}$. If
$\mu\leq\mu_{n,c}$, eigenvalues $\lambda_{\beta}$ are real with a negative
real part. This implies that the linear dynamical system (32) is stable.

The N equations (47) define what we call the "skeleton of the RFs map". In the
main text, we introduced the quantities
$r=\frac{\tau_{A}}{\tau_{B}},s=\frac{w^{-}}{w^{+}}$. In these variables, the
critical condition (47) reads:

$s_{n_{c}}=\frac{1}{w^{+^{2}}\,\tau_{B}^{2}}\,\frac{1}{4\,\kappa_{n}^{2}}\,\frac{\left(\,1-r\,\right)^{2}}{r^{2}}.$
(48)

This defines two critical lines symmetric with respect to $r=1$. These lines
are invariant by the variable change
$r^{\prime}=r,s^{\prime}=\left(\,\frac{{w^{\prime}}^{+}\,{\tau^{\prime}}_{B}}{w^{+}\,\tau_{B}}\,\right)^{2}\,s$.
This allows to map the skeleton obtained from a set of values
${\tau^{\prime}}_{B},{w^{\prime}}^{+}$ to the skeleton obtained with
references value $\tau_{B},w^{+}$.

##### Rectification.

In this paper we have essentially considered a situation where cells are not
rectified, whereas the full model (eq. (3)) considers rectification, in
agreement with realistic biological systems. The mathematical effect of
rectification of a cell ${A}_{j}$, is to set to zero the corresponding row in
the matrix $W^{{A}}_{{B}}$. This has several consequences. First, we cannot
apply anymore the useful Ansatz used in the section, that is
$\Gamma^{{A}}_{{B}}=\Gamma^{{B}}_{{A}}$. In addition, the vanishing of only
one row in ${\mathcal{L}}$ completely modifies its spectrum. However,
thresholding in rectification corresponds to partition the phase space of the
model, a compact subset of $\mathbbm{R}^{3N}$, into convex subdomains
delimited by hyperplanes. In each of these domains the matrix
${\mathcal{W}}.\vec{{\mathcal{R}}}(\vec{{\mathcal{X}}})$ appearing in eq. (32)
is linear with a number of zero eigenvalues corresponding to the number of
rectified cells. This matrix acts as a projector on the complementary subspace
of its kernel. In each of these subdomains eq. (10) applies. One can actually
compute, for a given stimulus, the time of entrance and exit in a new
subdomain with the effect of modifying the eigenvalues and eigenvectors
appearing in eq. (10). The resulting equation is quite complex though and will
require further investigations. See (Cessac, , 2020) for more details.

####  6.3.5 Nearest neighbours connectivity

We consider the case where the connectivity matrices
$\Gamma^{{B}}_{{A}}=\Gamma^{{A}}_{{B}}$ have nearest neighbours symmetric
connections. We also assume that the dynamics hold on a square lattice with
null boundary conditions. We define $\alpha=i_{x}\in\left\\{\,1\dots
L=N\,\right\\}$ in one dimension and
$\alpha=i_{x}+(i_{y}-1)\,L\in\left\\{\,1\dots L^{2}=N\,\right\\}$ in two
dimensions. We also set $n=n_{x}\in\left\\{\,1\dots L=N\,\right\\}$ in one
dimension and $n=n_{x}+(n_{y}-1)\,L\in\left\\{\,1\dots L^{2}=N\,\right\\}$ in
two dimensions. Then, the eigenvectors and eigenvalues of these matrices have
the form:

$\begin{array}[]{lll}\psi_{\alpha,n}=\left(\,\frac{2}{L+1}\,\right)^{\frac{d}{2}}\prod_{l}\,\sin\left(\,\frac{n_{l}\pi}{L+1}\,i_{l}\,\right),\\\
\kappa_{n}=2\,\sum_{l}\left[\,\cos\left(\,\frac{n_{l}\pi}{L+1}\,\right)\,\right];\end{array}$
(49)

with $l=x$ for $d=1$ and $l=x,y$ for $d=2$. Especially, in one dimension:

$\psi_{\alpha,n}=\sqrt{\frac{2}{L+1}}\sin\left(\,\frac{\alpha\pi}{L+1}\,n\,\right),$

The quantum numbers $\left(\,n_{x},n_{y}\,\right)$ define a wave vector
$\vec{k}_{n}=\left(\,\frac{n_{x}\pi}{L+1},\frac{n_{y}\pi}{L+1}\,\right)$
corresponding to wave lengths
$\left(\,\frac{L+1}{n_{x}},\frac{L+1}{n_{y}}\,\right)$. 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,1)}=2\left(\,\cos\left(\,\frac{\pi}{L+1}\,\right)+\cos\left(\,\frac{\pi}{L+1}\,\right)-2\,\right)$.
Each of these eigenmodes is related to a characteristic time
$\tau_{n}=\frac{1}{\lambda_{n}}$.

Eigenvalues $\kappa_{n}$ can be positive or negative. However, from eq. (43)
this has no impact on the eigenvalues as what matters is $\kappa_{n}^{2}$.
This induces however a symmetry $\kappa_{n}\to-\kappa_{n}$ that can be seen in
the skeleton Figure 4, not forgetting that this figure is in log scale.

## References

  * Asari and Meister, (2012) Asari, H. and Meister, M. (2012).  Divergence of visual channels in the inner retina.  Nature Neuroscience.
  * Baccus et al., (2008) Baccus, S. A., Ölveczky, B. P., Manu, M., and Meister, M. (2008).  A retinal circuit that computes object motion.  Journal of Neuroscience, 28(27):6807–6817.
  * Baden et al., (2016) Baden, T., Berens, P., Franke, K., Román Rosón, M., Bethge, M., and Euler, T. (2016).  The functional diversity of retinal ganglion cells in the mouse.  Nature, 529(7586):345–350.
  * Benda et al., (2007) Benda, J., Gollisch, T., Machens, C. K., and Herz, A. V. (2007).  From response to stimulus: adaptive sampling in sensory physiology.  Current Opinion in Neurobiology, 17(4):430–436.  Sensory systems.
  * Berry et al., (1999) Berry, M., Brivanlou, I., Jordan, T., and Meister, M. (1999).  Anticipation of moving stimuli by the retina.  Nature, 398(6725):334—338.
  * Besharse and Bok, (2011) Besharse, J. and Bok, D. (2011).  The Retina and its Disorders.  Elsevier Science.
  * Cessac, (2011) Cessac, B. (2011).  A discrete time neural network model with spiking neurons ii. dynamics with noise.  Journal of Mathematical Biology, 62(6):863–900.
  * Cessac, (2020) Cessac, B. (2020).  Recent Trends in Chaotic, Nonlinear and Complex Dynamics, in honour of Prof. Miguel A.F. Sanjuán on his 60th Birthday., chapter The retina as a dynamical system.  World Scientific.
  * Cessac, (2022) Cessac, B. (2022).  Retinal processing: Insights from mathematical modelling.  Journal of Imaging, 8(1).
  * Cessac et al., (2021) Cessac, B., Ampuero, I., and Cofré, R. (2021).  Linear response of general observables in spiking neuronal network models.  Entropy, 23(2).
  * Cessac et al., (2017) Cessac, B., Kornprobst, P., Kraria, S., Nasser, H., Pamplona, D., Portelli, G., and Viéville, T. (2017).  PRANAS: A new platform for retinal analysis and simulation.  Frontiers in Neuroinformatics, 11:49.
  * Cessac and Matzakou-Karvouniari, (2022) Cessac, B. and Matzakou-Karvouniari, D. (2022).  The non linear dynamics of retinal waves.  Physica D: Nonlinear Phenomena, 439:133436.
  * Cessac et al., (2009) Cessac, B., Rostro-Gonzalez, H., Vasquez, J.-C., and Viéville, T. (2009).  How gibbs distribution may naturally arise from synaptic adaptation mechanisms: a model based argumentation.  J. Stat. Phys,, 136(3):565–602.
  * Chalupa and Werner, (2004) Chalupa, L. M. and Werner, J., editors (2004).  The visual neurosciences.  MIT Press.  Two volumes.
  * (15) Chen, E. Y., Marre, O., Fisher, C., Schwartz, G., Levy, J., da Silviera, R. A., and Berry, M. (2013a).  Alert response to motion onset in the retina.  Journal of Neuroscience, 33(1):120–132.
  * (16) Chen, E. Y., Marre, O., Fisher, C., Schwartz, G., Levy, J., da Silviera, R. A., and Berry, M. J. (2013b).  Alert response to motion onset in the retina.  Journal of Neuroscience, 33(1):120–132.
  * Chichilnisky, (2001) Chichilnisky, E. J. (2001).  A simple white noise analysis of neuronal light responses.  Network (Bristol, England), 12(2):199–213.
  * Daw, (2012) Daw, N. (2012).  How Vision Works: The Physiological Mechanisms Behind What We See.  Oxford University Press, USA.
  * de Vries et al., (2011) de Vries, S. E. J., Baccus, S. A., and Meister, M. (2011).  The projective field of a retinal amacrine cell.  Journal of Neuroscience, 31(23):8595–8604.
  * Demb and Singer, (2015) Demb, J. B. and Singer, J. H. (2015).  Functional Circuitry of the Retina, volume 1.
  * Diamond, (2017) Diamond, J. S. (2017).  Inhibitory Interneurons in the Retina: Types, Circuitry, and Function.  Annual Review of Vision Science, 3:1–24.
  * Ermentrout and Terman, (2010) Ermentrout, G. B. and Terman, D. (2010).  Foundations of Mathematical Neuroscience.  Interdisciplinary Applied Mathematics. Springer.
  * Euler et al., (2014) Euler, T., Haverkamp, S., Schubert, T., and Baden, T. (2014).  Retinal bipolar cells: Elementary building blocks of vision.
  * FitzHugh, (1969) FitzHugh, R. (1969).  Mathematical models of excitation and propagation in nerve, chapter 1.  McGraw-Hill Companies.
  * Franke and Baden, (2017) Franke, K. and Baden, T. (2017).  General features of inhibition in the inner retina.  Journal of Physiology, 595(16):5507–5515.
  * Gollisch and Meister, (2010) Gollisch, T. and Meister, M. (2010).  Eye Smarter than Scientists Believed: Neural Computations in Circuits of the Retina.  Neuron, 65(2):150–164.
  * Hennig et al., (2009) Hennig, M., Adams, C., Willshaw, D., and Sernagor, E. (2009).  Early-stage waves in the retinal network emerge close to a critical state transition between local and global functional connectivity.  J Neurosci.
  * Hilgen et al., (2022) Hilgen, G., Kartsaki, E., Kartysh, V., Cessac, B., and Sernagor, E. (2022).  A novel approach to the functional classification of retinal ganglion cells.  Open Biology, 12(3).
  * Hilgen et al., (2017) Hilgen, G., Pirmoradian, S., Pamplona, D., Kornprobst, P., Cessac, B., Hennig, M. H., and Sernagor, E. (2017).  Pan-retinal characterisation of Light Responses from Ganglion Cells in the Developing Mouse Retina.  Scientific Reports, 7.
  * Hodgkin and Huxley, (1952) Hodgkin, A. L. and Huxley, A. F. (1952).  A quantitative description of membrane current and its application to conduction and excitation in nerve.  The Journal of Physiology.
  * Ichinose et al., (2014) Ichinose, T., Fyk-Kolodziej, B., and Cohn, J. (2014).  Roles of ON cone bipolar cell subtypes in temporal coding in the mouse retina.  Journal of Neuroscience, 34(26):8761–8771.
  * Kartsaki, (2022) Kartsaki, E. (2022).  How specific classes of retinal cells contribute to vision: a computational model.  Theses, Université Côte d’Azur, Inria, France ; Newcastle University.
  * Karvouniari et al., (2019) Karvouniari, D., Gil, L., Marre, O., Picaud, S., and Cessac, B. (2019).  A biophysical model explains the oscillatory behaviour of immature starburst amacrine cells.  Scientific Reports, 9:1859.
  * Manu et al., (2022) Manu, M., McIntosh, L. T., Kastner, D. B., Naecker, B. N., and Baccus, S. A. (2022).  Synchronous inhibitory pathways create both efficiency and diversity in the retina.  Proceedings of the National Academy of Sciences, 119(4):e2116589119.
  * Marr, (1982) Marr, D. (1982).  Vision.  W.H. Freeman and Co.
  * (36) Masland, R. H. (2012a).  The Neuronal Organization of the Retina.
  * (37) Masland, R. H. (2012b).  The tasks of amacrine cells.  Visual neuroscience.
  * McIntosh et al., (2017) McIntosh, L. T., Maheswaranathan, N., Nayebi, A., Ganguli, S., and Baccus, S. A. (2017).  Deep learning models of the retinal response to natural scenes.  Advances in neural information processing systems, 29.
  * Menz et al., (2020) Menz, M. D., Lee, D., and Baccus, S. A. (2020).  Representations of the amacrine cell population underlying retinal motion anticipation.  bioRxiv.
  * Morris and Lecar, (1981) Morris, C. and Lecar, H. (1981).  Voltage oscillations in the barnacle giant muscle fiber.  Biophys J, 35(1):193–213.
  * Nagumo et al., (1962) Nagumo, J., Arimoto, S., and Yoshizawa, S. (1962).  An active pulse transmission line simulating nerve axon.  Proceedings of the IRE, 50:2061–2070.
  * Pamplona et al., (2021) Pamplona, D., Hilgen, G., Hennig, M., Cessac, B., Sernagor, E., and Kornprobst, P. (2021).  Large visual neuron assemblies receptive fields estimation using a super-resolution approach.
  * Paninski, (2003) Paninski, L. (2003).  Convergence properties of three spike-triggered analysis techniques.  Network: Computation in Neural Systems, 14(3):437–464.
  * Pillow et al., (2008) Pillow, J. W., Shlens, J., Paninski, L., Sher, A., Litke, A. M., Chichilnisky, E. J., and Simoncelli, E. P. (2008).  Spatio-temporal correlations and visual signalling in a complete neuronal population.  Nature, 454(7207):995–999.
  * Protti et al., (2014) Protti, D., Di Marco, S., Huang, J., Vonhoff, C., Nguyen, V., and Solomon, S. (2014).  Inner retinal inhibition shapes the receptive field of retinal ganglion cells in primate.  J Physiol., 592(1):49–65.
  * Roth, (2016) Roth, B. L. (2016).  DREADDs for Neuroscientists.  Neuron, 89(4):683–694.
  * Sanes and Masland, (2015) Sanes, J. R. and Masland, R. H. (2015).  The Types of Retinal Ganglion Cells: Current Status and Implications for Neuronal Classification.  Annual Review of Neuroscience, 38(1):221–246.
  * Schreyer, (2018) Schreyer, H. M. (2018).  Nonlinearities in bipolar cells and their role for encoding visual signals.  Doctoral thesis.
  * Schröder et al., (2020) Schröder, C., Klindt, D., Strauss, S., Franke, K., Bethge, M., Euler, T., and Berens, P. (2020).  System identification with biophysical constraints: A circuit model of the inner retina.  bioRxiv, 2020.06.16.154203.
  * Schwartz et al., (2006) Schwartz, O., Pillow, J. W., Rust, N. C., and Simoncelli, E. P. (2006).  Spike-triggered neural characterization.  Journal of Vision, 6(4):484–507.
  * Simoncelli et al., (2004) Simoncelli, E. P., Paninski, L., Pillow, J., and Schwartz, O. (2004).  Characterization of Neural Responses with Stochastic Stimuli.  The cognitive neurosciences, page 20.
  * Souihel and Cessac, (2021) Souihel, S. and Cessac, B. (2021).  On the potential role of lateral connectivity in retinal anticipation.  J. Math. Neurosc., 11(3).
  * Urban and Roth, (2015) Urban, D. J. and Roth, B. L. (2015).  DREADDs (Designer Receptors Exclusively Activated by Designer Drugs): Chemogenetic Tools with Therapeutic Utility.  Annual Review of Pharmacology and Toxicology, 55(1):399–417.
  * Wässle, (2004) Wässle, H. (2004).  Parallel processing in the mammalian retina.